Intrinsic optical imaging has revealed a representation of eye position smoothly mapped across the surface of the inferior parietal lobule in behaving monkeys. We demonstrate here that blood vessels imaged along with the cortex have large signals tuned sometimes, but not always, to match the surrounding tissue. The relationship between the vessels and surrounding tissue in both space and time was explored using independent component analysis (ICA). Working only with single-trial data, ICA discovered a sequence of regions corresponding to the vascular propagation of activated signals from remote loci into the blood vessels. The vascular signals form a novel map of cortical function—the functional angioarchitecture—superimposed upon the cortical functional architecture. Furthermore, the incorporation of temporal aspects in optical data permitted the tuning of the inferior parietal lobule to be tracked in time through the task, demonstrating the expression of unusual tuning properties that might be exploited for higher cognitive functions.
Areas 7a and dorsal prelunate cortex (DP) of the inferior parietal lobe in a behaving monkey are crucial for generating spatial percepts (Siegel and Read 1997). The responses of single neurons in 7a and DP are modulated by both eye position and the position of the stimulus on the retina (Andersen and others 1985; Read and Siegel 1997). The topographic organization of gain fields (Siegel and others 2003), retinotopy (Heider and others 2005), and attention (Raffi and Siegel 2005) across the cortical surface has been established using intrinsic optical imaging of cortex in the behaving monkeys.
Intrinsic optical imaging uses small changes in reflected light from the surface of the illuminated brain to measure neuronal metabolism and, by extension, neuronal activity (Malonek and Grinvald 1996). Application of this technique with appropriate chronic recording techniques permit repeated measurements of functional architectures while visual, motor, and cognitive behaviors are performed.
Intrinsic imaging studies often examine regions devoid of substantial blood vessels or exclude blood vessels from study. In recent parietal studies (Siegel and others 2003; Heider and others 2004; Raffi and Siegel 2005), the vasculature was anecdotally observed to be tuned to the same variables as the neurons. If the vascular signal arose from local signal sources as generally accepted (Malonek and Grinvald 1996; Woolsey and others 1996; Malonek and others 1997; Logothetis 2003), then neighboring cortex and vessels should always have the same dependence on eye position. This was not always the case based on visual inspection. Presumably, the hemoglobin-containing blood cells that comprise a portion of the intrinsic signal constantly move, following tortuous paths defined by the vasculature and draining beds of venules and veins. These paths may lead the deoxygenated blood astray from its cortical origins.
The key issue examined here is whether the cortical signals and the closely opposed vascular fields have the same, or different, selectivity to stimuli. Although it is relatively straightforward to trace blood vessels by hand, it is not as easy to determine which segments have similar tuning. Similarly, when considering the cortical surface, segmenting the cortex into regions that have similar tuning is not a trivial task. One could draw boundaries based on some tuning measure for some time slice on a pixel-by-pixel basis, and perhaps this process could be automated; however, there can be a considerable arbitrariness in the decision process. Thus, the question of the similarity between cortical signals and the draining vascular fields needs to be approached in 2 stages. First, the relevant regional cortical fields need to be computationally defined. Second, the tuning between fields requires comparison.
In the first step of the data analysis, intrinsic cortical from hundreds of millions of data points in space and time were segmented into statistically independent components to explore in detail the relationships between the tissue and vascular beds in the association cortex of behaving monkey. With this approach, called “independent component analysis” (ICA) (Bell and Sejnowski 1995), the data were grouped into collections of pixels with similar informational content (Duann and others 2002). ICA is blind to any spatial or temporal continuity so that the structures that are discovered are not constrained by any assumptions other than that the resulting components are spatially independent. The spatial contiguity and temporal course of the gain fields that are extracted from the data were not preordained by a priori expectations. The second step was to analyze the components for physiological properties. As expected, the cortical patches had coherent temporal signals and were tuned to eye position over time. Segments of blood vessels also were tuned to eye position in time. Comparison with nearby regions revealed that the cortical patches did not always match the dependence on time or on eye position of the nearby blood vessel signals, raising questions about the interpretation of functional imaging data collected at lower spatial resolutions.
Monkeys were prepared for chronic optical imaging using established methods (Siegel and others 2003; Heider and others 2004; Raffi and Siegel 2005). Briefly, substantial head holders were used to provide stability between the camera and the animal's skull to the 1-μm level, even while the monkey was performing a reaction time task. Eye position was monitored with an infrared eye tracker to within 1°. The monkey performed a gain field mapping task (Read and Siegel 1997) for which he fixated a 0.5° spot and attended to an expansion flow pattern. Release of a lever when the motion became unstructured led to a juice reward. All studies were carried out according to National Institute of Health Guidelines for Animal Research and approved by Rutgers University Animal Care and Facilities Committee.
Optical signals were collected while the cortex was illuminated with a 100-W halogen light source powered by a stabilized direct coupled source bandpass filtered at 605 nm. The wavelength 605 nm emphasizes deoxygenated hemoglobin changes (Malonek and others 1997; Vanzetta and others 2004). An Optical Imaging Company (Rehovot, IL) Imager 2001 system was used to collect the images. Every 16 trials, a reference image was collected over 256 frames at 30 Hz while the monkey was not under behavioral control; this reference is specific to the Imager 2001 system and is only used to increase the number of bits of the analog-to-digital conversions. Data frames relative to the reference were collected at 2 Hz synchronized to the fixation onset. Off-line, the differences from the reference image were then added to the reference image to yield a signal with a final resolution of ∼16 bits. These images in arbitrary luminance units were then stored for further analysis (Siegel and others 2003). Only images from behaviorally correct trials were used. Spatial resolution was 30-μm spatial per pixel over an ∼8 × 11–mm domain.
The order of stimulus presentation was randomly selected in a fixed block design by the computer controlling the behavior (Read and Siegel 1997). Thus, the monkey never knew what the next stimulus would be until the trial began.
Independent Component Analysis
Images with 360 × 240–pixel resolution were collected. Earlier studies used a baseline normalization approach to eliminate time by subtracting a baseline response from a visually evoked response (Siegel and others 2003). Here, ICA was used in order to exploit the spatial and temporal information in the optical signals. ICA segregates data based on the information content in the temporal stream (Bell and Sejnowski 1995; Duann and others 2002).
The spatial independence assumption made in applying ICA to the optical imaging data is consistent with the principle of brain modularity, namely, that different brain regions perform different functions, with different time courses of activity (though not necessarily independent, particularly when only a few hundred or fewer time points are available). Spatial modularity, plus the high spatial resolution of optical images, allows the use of ICA to identify maximally “spatially” independent regions with distinguishable time courses. Specifically, the input matrix to ICA, xj,k, is the optical image data used in ICA training, where j = 1, …, T is number of time frames and k = 1, …, N is the number of pixels in each frame of optical image. k is computed as k = Iη + J, where (I, J) are the indices of the pixel and η is the number of pixels in a row. Note that the ICA algorithm does not have any knowledge of the pixel location parameters (I, J, η).
ICA finds an “unmixing matrix,” “W”, to perform component separation and recover the underlying independent sources, ui,k = Wi,j × xj,k, where i = 1, …, M is number of components and M ≤ T. The brain activity of a region of interest can be obtained by projecting selected ICA component back onto to the original data space, x′j,k = Wj,i−1 × ui,k. Components are ordered according to the mean energy of back-projections: ξi = meank(∑jx′j,k2).
For each session, ICA training data consisted of ∼300–900 concatenated 8-point, 360 × 240–pixel epochs (∼30–100 trials for each condition). The mean of all images (2400–7200 images) was subtracted from each image. In the analysis, the average image was computed for all frames (as more clearly indicated in the text) and subtracted from all frames to remove the DC offset Ak = ∑j = 1Txj,k. Otherwise, the first principal component analysis (PCA) component would be this DC offset (Ak ∼ A(I, J)). The resulting analysis was then performed on signals considered as a percentage change from the average of all frames. The use of the percentage change permits comparison of the amplitude signal with other published studies.
The percentage change spatial–temporal data was then processed with PCA to reduce the temporal dimensionality of the data from 2400–7200 temporal points to 200 principal components that account for 96.3 ± 0.7% (n = 17 experiments) of the variance of the data. This reduction was necessary to guarantee the convergence of ICA training. The dimension of the unmixing matrix W should, based on our empirical results, be smaller than where Ninput is the number of voxels. Using a typical set of original image data, this unmixing matrix is approximately 16 000 × 16 000 elements. Using a PCA reduction of the time values to 200 components, the unmixing matrix becomes a more manageable ∼4000 × 4000 elements.
The PCA components typically selected blood vessels but never patches of cortex, as described in Results. The 200 principal components were then back-projected to the original space to be analyzed by ICA, yielding a 200-element vector for each of 360 × 40 pixels. The 200 element per pixel vectors were randomly presented to the ICA algorithm (i.e., ICA could not use any spatial information). ICA is an iterative procedure and converged when the 200 independent components were maximally spatially independent of each other. Of these, usually 100 components had spatially coherent and clearly defined regions indicating a putative biological source, with the remainder having scattered “noisy” signals.
Comparisons of Regions of Activity
Regions of activity (ROAs) were computed from the “mixing matrix.” The mean and standard deviation of the values in the mixing matrix were computed, and a Z-score was assigned to each value. Pixels were included in a region of interest for a component if the absolute Z-scores for the values in the mixing matrix were greater than 2 (|z| > 2). These Z-scores are color coded in figures of ROA.
Similarity Measures of ROAs
The similarity of the ROA between experimental runs was computed using an overlap ratio measure (ORM). ORM, as defined in equation (1), measures the overlap between the ROAs from 2 component maps. Corresponding component maps from different data sets are coregistered by shifting one of the maps to spatially fit the other according to landmarks, which can be easily found in the image. In our data set, the large intraparietal vein lying across the entire image horizontally and some surrounding arterioles were selected as landmarks to coregister 2 component maps. After component maps were well coregistered to each other, the significant pixels (|z| > 2.0) repeatedly found in both component maps at the same locations were counted as indicated in numerator of equation (1). The number of commonly found pixels was then divided by the square root of the multiplication of the total pixel number of the 2 ROAs to yield the ORM expressed as a percent.
Selection of Independent Components Based on Contiguity Measures
ROAs were collections of pixels that exceeded a certain threshold. These pixels could be contiguous or scattered across the image. To assess the contiguity for each pixel, a defining independent component that maximally contributes to the activity (time course) of the pixel was selected. This was determined by looking at the matrix ui,k = Wj,i × x′j,k, where W and x are defined as above. Each element ui,k defines the contribution of the ith independent component to the kth pixel. So, throughout the kth column of u, the element i* is defined such that ui*, k is maximum. Hence, i* is the component that maximally contributes to the kth pixel. Because each pixel can only have a single maximum across its independent components, each pixel is labeled with only one independent component. This provides a segmentation of the image by the independent components. Then a 3 × 3 mask was superimposed upon each pixel. If more than 4 pixels had the same components, this pixel was considered to be part of a contiguous component. If fewer than 4 pixels were the same, the pixel was treated as noise and dropped from the analysis. The remaining pixels were then used to create a histogram of the number of contiguous pixels for each component.
Determination of Gain Field Tuning
In order to determine how each component was tuned with respect to the varied eye position, linear regression was performed upon the trial-by-trial components data using a standard general linear model with PROC GLM (SAS Co., Durham, NC).
Before plotting, the signs of these regression terms were multiplied by −1, so the direction of the amplitude of the components matched those expected from electrophysiological measurements (Siegel and others 2003). The term ϵ(I, t, i) is the residual error.
As noted elsewhere (Siegel and others 2003; Heider and others 2004; Raffi and Siegel 2005), the low signal-to-noise in intrinsic optical data results in low amounts of the variance accounted for by selected models. Two questions arise: First, what is the appropriate model for the data? Second, does the data deviate significantly from a random distribution?
The first question was addressed by comparing different models. A criterion for comparing models based on the sum of squares of error is a poor criterion because increasing the number of parameters necessarily decreases the sum square error (Siegel and Birks 1988; Heider and others 2004). The Akaike Information Criteria (AIC) penalizes for the increase in the number of parameters and has been used to select between different classes of general linear models (Siegel and Birks 1988; Heider and others 2004). The AIC was computed for the results of regression using equation (2) as well as a second order with interaction regression equation:
This was performed for every time slice for each component. It was found that for the substantial majority of components and the majority of time points, the linear model was a better representation of the data using the AIC and was hence used throughout this study.
Significance of Regression Model
The correlation coefficient (or R2) values are low when using optical data on a pixel-by-pixel basis as noted in our earlier work. The correlation coefficient was also low when evaluated for the independent components. To determine whether the gain fields that are computed from the regression coefficients were significantly different from zero, a Monte Carlo analysis was used (Siegel and others 2003). In short, the regression coefficients and hence the vector (βx, βy) were computed using equation (2) for half the data selected at random (without replacement) from one experimental data set. This was repeated for 500 random selections. The distribution of gain field vectors was computed and compared with that expected for a uniform distribution using a circular bivariate statistic, Hotelling's 1-sample t-test (Batschelet 1981). This statistic uses the direction as well as the amplitude of the response. This same analysis was performed after randomizing the trial-by-trial relationship between the stimuli and the measured components, using half of the original data set each time. For the latter case, the circular bivariate statistic indicated that the gain fields obtained with the shuffled data were not significantly different from a uniform distribution with a mean gain field amplitude of zero.
It should be noted that these analyses were performed on a time slice by time slice basis across all components for all experiments. The voluminous data were reduced and represented as needed in the accompanying figures.
The gain field in inferior parietal lobule was studied in 2 monkeys performing a reaction time task (Siegel and others 2003) during which they fixated a small red light and detected changes in motion stimuli (Fig. 1). In this task, the monkey pulled back a lever at the onset of a 0.1° red point and simultaneously fixated the point. The eye position was monitored with an infrared system, and eye movements terminate the trial. After 2000 ms, the monkey was presented with a 10° expanding optic flow stimulus over the red fixation point. At a random time 3000–4500 ms after the optic flow onset, the stimulus became unstructured (Read and Siegel 1997), and the monkey had to release the key within an 800-ms reaction time window for a juice reward. There was a 5000- to 7000-ms delay between trials that permitted the optical signal to decay. Further, the order of the stimulus was randomized (Methods). The vertical and horizontal position of the fixation points were varied across trials, whereas the visual stimulus was always presented over the fixation point, to assess the effect of eye position on the visual response (i.e., the gain fields). As the animal fixated 1 of 9 positions, expanding optic flow was presented over the fovea. In 2 monkeys, 3 cortical areas were simultaneously imaged: caudal parietal area (PEc) (Raffi and others 2002), area 7a, and the dorsal prelunate cortical gyrus (DP) (Andersen and Siegel 1990).
Gain field maps, as described in an earlier study, suggest that there is a gradual modulation of the gain eye signal with position on the cortex. ICA was used to computationally segment the maps derived from the optical data using both temporal and spatial information. The ICA identifies groups of pixels that behave similarly in time; these groups are termed “independent components.” The independent components are weighted sums of the original data that are uncorrelated from each other in an information theoretic sense (Bell and Sejnowski 1995; McKeown and others 1998; Duann and others 2002). The particular implementation used here treated each pixel autonomously, and the algorithm had no knowledge of the location or time for the input data. As a result, groups of pixels in each independent component were matched in their properties solely as a result of their maximal mutual information and minimal entropy.
Description of Independent Components
An examination of the independent components indicated that the cortical surface was segmented into 3 types of regions that covered the imaged region: 1) irregularly shaped contiguous patches that overlaid the cortex (Fig. 2A1,B1), 2) thin branching components that overlaid blood vessels (Fig. 2C1,D1), and 3) scattered isolated pixels. (Forty of the 200 components are illustrated in Fig. 3.) The “component number” was assigned based on ordering by mean energy, ξi (see Methods).
In order to select independent components, 2 approaches were used. First, independent components were identified “by eye,” eliminating those that clearly consisted of scattered pixels. Second, the contiguity of pixels was calculated as described in Methods (Fig. 4). Each pixel was evaluated to determine which component maximally contributed (Fig. 4A). This effectively segmented the image into many regions. The noisy components are actually in this image but cannot be easily observed. Each pixel was thus labeled with an independent component. If a pixel for a particular component had at least 3 neighboring identical components (out of the 8 neighbors), it was added to the component's bin; this resulted in a distribution of the number of pixels with contiguous pixels (Fig. 4B). Components with more than 50 pixels per bin were examined further. These matched well those selected by eye.
The shapes of “patch” and “vessel” components were consistent across 6 weeks of experimentation; the components measured twice within a day were also internally consistent as shown by subdividing data sets into 2 and repeating the analysis (Fig. 5). Although the component number (ordered by ξi) for each particular patch of cortex could vary across days, the match of the shape and location of the patch or blood vessel across the cortex was remarkable. The reproducibility of the patches across days and within a days' experiment argues against these patches arising from an artificial methodological segregation of the continuous gain field maps.
To quantify the similarity of the independent components and their spatial locations for the brain patches and blood vessels obtained in ICA component maps between experiments, an ORM was devised (Duann and others 2002); matched patches could easily be found between experiments, and there was a high degree of spatial correlation between them (ORM = 82.7 ± 4.7%, N = 10). The reproducibility of the patches across experiments with different noise suggests that the patches were not critically dependent on the numerical analysis, consistent with a biological origin.
Temporal Ordering of Cortical and Vascular Activation
Clues to the origin of these independent components were found in their time courses, by examining the temporal order by which nearby independent components (ICs) and branches were activated. For many pairs of patches and vessels, the vessel was activated with a delay relative to a patch that was spatially close to a vessel. This sequence of the responses of the patch and the vessels was consistent with the transfer of deoxygenated hemoglobin from cortex to nearby vessels (second columns of Figs 2 and 6). If this was indeed so, then the eye position tuning of the patch ought to be transmitted through the deoxygenated hemoglobin to the putative draining blood vessel (Chaigneau and others 2003). There should be a match between the tuning of the patch and the nearby draining vessel. In order to determine if each component was tuned with respect to the varied eye position, components were explicitly reordered using the stimulus condition after linear and quadratic regressions (eqs 2 and 3) determined the dependence of each component on the vertical and horizontal eye position. The regressions defined the gain field in time (see Methods).
Time Course of Gain Field
The characteristics of the time course of the ICs provided insights into the relationship between the patches and the vessels. The gain fields for patches were visualized in time by plotting (βx, βy) as a vector in time (second columns of Figs 2 and 6). The tuning observed in the optical signal prior to stimulus onset unexpectedly contained information about eye fixation and corresponds to the so-called “pure” fixation signal observed with single-unit recordings in area 7a (Andersen and Siegel 1990; Read and Siegel 1997). Over time, the strength of the gain field response increased, and its direction evolved. The final tuning of the patch matched that obtained using temporal averaging in an earlier study (Siegel and others 2003).
Comparison and Significance of Models
The linear regression model was a better representation of the data than the quadratic regression model, based on the AIC (see Methods). The AIC was compared for each time slice of the linear and the quadratic models (eqs 2 and 3). Because there were 8 time slices and 200 components, 1600 measurement regressions were made in the experiment as shown in Figure 2. Of these 1600, 1350 (84.4%) were better modeled by a linear component than by a quadratic. When the linear coefficients were compared for the 2 models, they were within 2% of each other and highly correlated (Fig. 7A). Similar findings were obtained with other experiments; therefore, the linear model was used to represent the gain fields.
The correlation coefficients for these linear regressions, like those of earlier studies (Siegel and others 2003; Heider and others 2004), were quite low (R2 = 0.013 ± 0.012, mean and standard deviation of 8 time slices by 200 components for experiment 03-19-2000/run 4), as expected given the low signal-to-noise of optical data. Nonetheless, these gain fields were shown to be significant at each time point by the use of a Monte Carlo approach (see Methods, Siegel and others 2003). When the trial-by-trial components were randomized with respect to the stimulus conditions, the amplitude of the gain fields statistically vanished. When 500 shufflings were performed, the mean amplitude of the gain field across time was 0.05 ± 0.03 ADC per degree, n = 8, where ADC are raw units of the analog-to-digital converter. For comparison the mean amplitude of the unshuffled data was 1.77 ± 0.57 ADC per degree, n = 8. The shuffled data appeared as a cluster of points close to the origin (Fig. 7B). Such an analysis was performed on every time slice of every component. For the experiment of Figure 2, 1598/1600 gain fields were significantly different from a uniform distribution of gain fields at a level of P < 0.001 (Hotelling's 1-sample test, Batschelet 1981). (Note that due to the computationally intensive nature of this shuffling analysis, the raw data were used rather than normalizing to a percentage change signal.)
Time Dependence of the Gain Field
Because the linear model was the best representation of the data, it was used to examine the time course. The gain fields in many patches were found to develop in 2 stages: In Figures 2 and 7, there was initially a development of tuning in the vertical direction followed by a later tuning to the contralateral field. The curving trajectory of the gain field with time suggests that the vertical and horizontal modulation of the vascular signal occur with different time courses. This was a common finding for both the cortical patches and the vascular patches, suggesting multiple streams of eye position information arriving to the upper layers of the imaged area 7a and DP or ongoing modulation of eye position by putative cognitive or sensory signals.
That this curving trajectory was not simply a result of noise in the signals can be addressed in 2 ways from the biological data. First, within a day, nearby components that drained into apparently the same vessel had the same time course (third column of Figs 2 and 6). Second, a particular component had the same tuning across days. In Figure 8, 2 components with a high similarity in shape measured from 2 different days were assessed for their gain field tuning. The minimum and maximum of the range of the intercepts and the slopes were scaled to be (0, 1). The gain field tuning for both components was remarkably similar, switching from favoring upper eye positions to lower eye position and back again, with a contralateral shift.
Comparison between Vascular and Cortical Signal
For many of the vessels, there was a reasonable match between the draining vessel and the surrounding cortex, suggesting that the vascular signal is a good indicator of the cortical deoxygenation signal and by inference the neuronal activity. A fortuitously close anatomical arrangement of 4 blood vessels near the intraparietal sulcus illustrates a severe violation of these assumptions (Fig. 9). The upper vessel appears to collect blood from the superior parietal gyrus (Fig. 9B), presumably PEc (Raffi and others 2002), whereas the lower vessel likely collects blood from the inferior parietal gyral area 7a (Fig. 9C). The source of the blood for the intermediate vessels is not apparent (Fig. 9D,E). Each of these vessels has different delays in the time course for the amplitude of the components (cf., Fig. 9Bα,Cα). The gain field tuning for eye position also differed (cf., Fig. 9Bβ,Cβ). A functional magnetic resonance imaging (fMRI) measurement would mix and produce some weighted average of these signals, coming from at least 2 different nearby cortical areas being combined. In Figure 9F, a hypothetical signal was used to represent a coarse spatial average as might be found with fMRI. It was constructed by combining the trial data for all 4 components and then computing the intercept (Fig. 9Fα) and the gain field tuning (Fig. 9Bβ). These hypothetical measurements do not match any of the individual components. Any imaging approach with lower resolution than these vessels would likely report a combined signal or signal a gain field gradient in space.
PCA of Optical Data
Blind separation is a stringent test for locating blood vessels and separating them from cortical regions. It was possible that the pixels representing the blood vessels contained signals that were so dominated by the light scatter of the nearby tissue that they would be segregated by ICA, but this did not occur. Could another method like PCA suffice? When we performed a principal component decomposition of the optical data, only the blood vessels were found in the principal components; this linear analysis failed to find the patches of cortex (Fig. 10). This occurred because PCA is dominated by the largest signals and tends to mix signals that have low power, such as the cortical patches. The differences between ICA and PCA appear to be crucial in providing the segregation needed to explore the relationship between the vascular and tissue signals we describe here.
The goal of this study was to examine the relationship between metabolic signals in cortical tissue and nearby vasculature in the inferior parietal lobule of the behaving monkey. This was accomplished first by segmenting the hundreds of millions of pixels in space and time from an experimental session into regions using an unsupervised method and second by comparing the tuning of nearby regions. ICA was used to identify the spatial–temporal sets of pixels with similar informational content. The analysis revealed short segments of blood vessels and multiple spatial patches within the apparently continuously mapped cortex.
The tuning properties of the cortical patches and vascular segments were determined using regression analysis. The tuning between the cortex and nearby vessels and cortex did not always match. Independent patches were found proximal to blood vessels that were not apparent in the original gain field studies (Siegel and others 2003). The patches and vessels identified by the ICA were remarkably similar across days and within a days' experiment. This was clear both from visual inspection and from the ORM. This makes it highly likely that the patch size and shape revealed by ICA reflects the underlying biological processes.
Comparison with Other Methods
PCA has been successfully used to remove blood vessels and other artifacts (e.g., from respiration) from imaging data. Varimax (Kaiser 1958) and Promax (Hendrickson and White 1964) further segregate the data by rotating principal components toward a “simpler structure” to concentrate the variance into relatively few pixels or few time points. However, these methods are computation intensive when the number of principal components is large, and convergence is difficult to achieve with 200 principal components (Duann and others 2002). Varimax is further limited to orthogonal rotations in the principal component subspace (Mocks and Verleger 1986) and cannot account for activity from nonorthogonal brain sources (Donchin and others 1986).
Another state-of-the-art approach to blind separation is extended spatial decorrelation (ESD), which has been used to extract orientation patterns in V1 (Schiessl and others 2000; Stetter and others 2000). ESD requires that the data be averaged over repeated presentations grouped by the stimulus condition. By using averaged data, the size of the data set and the amount of computation are reduced but at the expense of altering variance estimates and parameter significance. In contrast, ICA uses all the raw data, without spatial filtering or averaging, which is a substantial advantage for low signal-to-noise measurements, such as those found with optical imaging and fMRI; fewer repeats allow a broader sampling of stimulus space, which makes it possible to investigate more subtle interrelationships between the temporal and spatial aspects of cortical processing.
Finally, and perhaps most importantly, the present approach does not assume a priori knowledge of the stimulus conditions so that the segregation is indeed “blind.” Spatial or temporal filters are useful when the stimulus or cognitive dimensions are known a priori (e.g., orientation tuning in primary visual), but this approach may miss novel and unexpected features. The disadvantage of using ICA is that the underlying sources of the spatial and temporal components may not be known and may require additional biological evidence to interpret them (Brown and others 2001; Jung and others 2001; Lee and others 2002; Makeig and others 2002; Tang and others 2002).
Source of the Blood Flow Signal
Sequences of activation in tissue and vasculature have been demonstrated in rats using optical methods (Woolsey and others 1996; Erinjeri and Woolsey 2002; Sheth and others 2004). Malonek and Grinvald (1996) used spectral analysis to separate the contribution of the hemoglobin and volume derivative signals, demonstrating the oxygenation overshoot. Others have examined changes in blood flow at the capillary level, providing evidence that blood can be shunted specifically to microregions of cortex.
Our results extend these findings to association cortex in primates and uncovers evidence for mesodomains (“meso” here refers to the 500- to 4000-μm scale) of vascular draining, which are unrelated to a columnar organization as in striate cortex (Vanzetta and others 2004, 2005). The mesodomains do not necessarily match the draining vessels and can vary over time. Indeed, the spatial distribution of tuning can abruptly change where blood vessels are in close opposition.
Temporal Modulation of Gain Fields
The temporal tuning of the parietal cortical patches identified by ICA was unexpected. In early visual cortex, the time dependence for the orientation tuning of patch of cortex develops smoothly, without any shift from the preferred orientation (Malonek and Grinvald 1996). In contrast, the vertical and horizontal contributions to the gain field varied, as might be expected if caused by activation and inactivation of multiple processes.
There are several possible sources for the signals found in the inferior parietal lobule. The upper layers of area 7a and DP cortex receive feedback from multiple areas (e.g., area 7a receives feedback from prefrontal areas 46 and 8a, as well as the anterior superior temporal polysensory area, STPa) (Siegel and Read 1997). The optical signals at 605 nm are measurements of the amount of deoxygenated hemoglobin. The primary source of the deoxygenated hemoglobin intrinsic signal is from the neural tissue that has the maximum metabolic impact, namely, the smallest afferent fibers and dendrites with their associated large surface area to volume (Malonek and Grinvald 1996; Woolsey and others 1996; Malonek and others 1997; Erinjeri and Woolsey 2002; Logothetis 2003; Sheth and others 2004). One possible explanation for the 2 stages of the time course is different temporal components from the 2 projective areas with afferents to layer II/III. For example, STPa might contain faster vertical signals, whereas area 46 has slower horizontal eye position signals. Differences in timing of laminar contributions could also have a role (Woolsey and others 1996). Alternatively, the early fixation signal and the later gain field responses could have different temporal tuning; electrophysiological data to date suggest no relationships between the two (Read and Siegel 1997). Lastly, there may be shifts in blood flow between different capillary beds as a function of the neural activity (Chaigneau and others 2003) or indeed between different layers of cortex (Woolsey and others 1996).
A Hidden Functional Architecture?
ICA also revealed that the cortical patches have a dimension of as much as 2 × 4 mm, which is larger than the 1-mm columns of striate cortex and other areas but commensurate to stripes in area V2. There are 2 possibilities for the origin of the segmentation of the map. First, there could be an underlying functional neural architecture that the particular stimulus and behavioral paradigm used here fails to reveal. There might be a different behavioral task or visual stimuli that would lead to the segregation of function along the patch's borders. In optical studies of spatial attention, repeating patches of 860 μm are found embedded within the eye position gain field map (Raffi and Siegel 2005); however, these patches are spatially very different in character from the mesodomains described here. Hence, at least the attentional architecture cannot account for the mesodomains discovered by ICA.
Second, these maps could be segmented based on angioarchitectonic principles, that is, the organization of the vasculature constrains the neuronal maps. The latter seems more plausible because many of the patches could be matched up with a nearby blood vessel defined by the ICA (e.g., Figs 2 and 6) with 500- to 1000-ms lags between them. Further, when the gain field tuning of these spatially neighboring components were compared, they often matched well to its nearby cortex. The temporal lag in the time course and the match in the gain field tuning suggested a 2-step process: first, an increase in neural activity, deoxygenating the surrounding microvasculature (Malonek and Grinvald 1996); second, this deoxygenated blood in the surrounding extracellular space and microcapillaries is collected into the draining veins.
Such a straightforward transfer of deoxygenated blood between the cortex and vessels would be expected to have no real effect on the signals assessed in the 1- to 5-mm range. Thus, imaging studies with a lower resolution system such as fMRI for which the vascular signals might dominate would still arrive at conclusions that would be valid for the nearby drained cortex. However, not all patches and nearby blood vessels were matched in their tuning. For some vessels, there appeared to be absolutely no relationship to the nearby cortex or nearby vessels, and the temporal relationship could be disrupted. This suggests that the draining field for that blood vessel is elsewhere, for example, either deep in the cortex or spatially remote.
In primary visual cortex of awake monkey, a different relationship between the vasculature and the nearby cortex has been reported (Vanzetta and others 2004, 2005). Imaging at the same 605-nm wavelength does indeed show that blood vessels appear tuned for orientation tuning; however, the relationship between the cortex and the vessels is less clear. Borders of draining fields, such as those seen in Figure 2 or 6 of the current work, were not observed in Vanzetta's work.
These differences may be a result of differences in the analysis. In the prior V1 study (Vanzetta and others 2004), the cortex and vasculature were segregated according to a PCA; here the ICA selected regions in the images with similar content defined in informational terms. In both studies, the cortical signals generally preceded the vascular components. However, in the V1 study, the vascular time course was the same for every blood vessel in the images. In the current study, different blood vessels are independently identified. Indeed multiple capillaries joining together to form a larger vessel were segmented by our approach. Another difference was that the cortical signals in the V1 study were for the entire imaged cortex; in the present study, the cortical region was subdivided by ICA into multiple patches and multiple segments of blood vessels. As a consequence, many pairings of patches intimately associated with draining vascular fields can be identified, providing richer information about the hemodynamics.
Alternatively, differences between the 2 studies may arise from the biological characteristics of the cortical regions studied. In monkey V1, the functional architecture is strongly determined by the underlying neural circuitry, such as orientation or ocular dominance. In contrast, in the imaged portions of the inferior parietal lobule, the linear nature of the gain, the very large receptive fields, and the gain field functional architecture suggest that almost all of its cortex will be altered by a stimulus (Read and Siegel 1997; Siegel and others 2003; Heider and others 2004). Hence, the more distributed nature of the gain field representation may provide a better platform to reveal the vascular effects.
Implications for Functional Brain Imaging
This lack of a match between spatially close vessels (Fig. 9) has important implications for fMRI studies, in which the pixel size is minimally 200- to 300-fold greater than the 30 μm in the optical studies. In the fMRI studies, a basic assumption is that the signal directly represents the neural activity at that location (Vanzetta and Grinvald 1999; Logothetis 2002). Further, it assumes no mixing between adjacent regions. The lack of spatial resolution will lead to incorrect interpretations of different signals across adjacent vessels as a cortical topography, a presumed gradient of function in an array of neurons when, in fact, there is none.
Although the substantial contribution of the vasculature to the fMRI signal might seem at first to create difficulties for the functional human imaging experimentation, these effects can be exploited with a combination of local measurements of variability and angiography. ROAs arising from multiple blood vessels ought to be identifiable by independent components with high variance or broad tuning. An additional measure needs to be incorporated to distinguish noise from vascular signals in functional imaging experiments. These highly variable voxels could be correlated with noninvasive high-resolution angiography down to the submillimeter range (Bernstein and others 2001; Reichenbach and Haacke 2001; Schad 2001; Hall and others 2002). Simulation of the hemodynamic flows through the 3-dimensional venous drainage could be tempered by local functional estimates of time-dependent tuning. The progression of cortical tuning from the tips of the smallest venous vessels to the larger vessels could be thus accurately rendered. Two maps should emerge from such an analysis: one of function mapped onto vasculature, perhaps at the submillimeter scale, and the second would be of the cortical tissue itself.
This functional angioarchitecture map (fAAM) would be constrained by a combination of the vascular regions segregated by independent component and detailed physical angiography. The second functional corticoarchitecture map (fCAM) would depict the function mapped onto cortex, constructed by excluding reconstructed vasculature signals of the functional angioarchitecture. The fAAM could be used to assess risk for stroke by the presence of singularities or discontinuities. The impact of intracranial aneurysms and carcinomas on function could be similarly imaged by exceptional increases or decreases in the fAAM. Therefore, the detail afforded by the ICA of the vascular-based signals should provide a powerful new means to map higher cognitive function in human and nonhuman primates.
The eye position gain field of the inferior parietal lobule was imaged with 605-nm light to measure deoxygenated hemoglobin signals, and the spatiotemporal data were segmented into components that had maximal information independence. The components corresponded to segments of blood vessel and cortical mesodomains that were tuned to the eye position. Matches in the tuning of cortex and nearby vessels were often, but not always, found with an appropriate delay, indicating draining of cortical deoxygenated blood by vessels. Lower resolution fMRI signals in the vicinity of blood vessels may be biased by vessels that are not matched in their tuning properties to the neighboring cortex, necessitating care in functional interpretation. The vascular signals comprise a fAAM as well as a segmented fCAM of the cortex and could be exploited with high-resolution imaging techniques to reveal neurological dysfunction.
Early discussions of blood flow issues with Gabor Jando as well as his contributions to the published data used in this study are acknowledged. The use of these published data also collected by Milena Raffi, Raymond Phinney, and Jessica Turner is acknowledged. Careful reading and discussion of this study by Larry Cohen and Scott Makeig is gratefully acknowledged. This work was supported by the National Institutes of Health Grants EY-09223 (RMS) and 1S10RR-12873 (RMS), The Whitehall Foundation (RMS), National Science Foundation Grant NPACI RUT223 (RMS), Howard Hughes Medical Institute (TS), and Swartz Foundation (J-RD, T-PJ, and TS). Computational resources of the Center for Computational Neuroscience, Rutgers, Newark, are acknowledged as are the massive file transfer and storage services provided by the Storage Resource Broker team headed by George Kremenek at the San Diego Supercomputer Center. Conflict of Interest: None declared.