Independent component analysis of high-resolution imaging data identifies distinct functional domains

www.elsevier.com/locate/ynimg NeuroImage 34 (2007) 94 – 108

Independent component analysis of high-resolution imaging data identifies distinct functional domains Jürgen Reidl,a,b,1 Jens Starke,a,b,e David B. Omer,c Amiram Grinvald,c and Hartwig Spors a,d,⁎,1 a

Win Group of Olfactory Dynamics, Heidelberger Akademie der Wissenschaften, Germany Interdisciplinary Center for Scientific Computing and Institute of Applied Mathematics, University of Heidelberg, Im Neuenheimer Feld 294, D-69120 Heidelberg, Germany c Department of Neurobiology, Weizmann Institute of Science, Rehovot 76100, Israel d Max-Planck-Institute for medical research, Jahnstrasse 29, D-69120 Heidelberg, Germany e Department of Mathematics, Technical University of Denmark, Matematiktorvet, Building 303, DK-2800 Kgs. Lyngby, Denmark b

Received 7 October 2005; revised 10 August 2006; accepted 13 August 2006 Available online 25 October 2006 In the vertebrate brain external stimuli are often represented in distinct functional domains distributed across the cortical surface. Fast imaging techniques used to measure patterns of population activity record movies with many pixels and many frames, i.e., data sets with high dimensionality. Here we demonstrate that principal component analysis (PCA) followed by spatial independent component analysis (sICA), can be exploited to reduce the dimensionality of data sets recorded in the olfactory bulb and the somatosensory cortex of mice as well as the visual cortex of monkeys, without loosing the stimulusspecific responses. Different neuronal populations are separated based on their stimulus-specific spatiotemporal activation. Both, spatial and temporal response characteristics can be objectively obtained, simultaneously. In the olfactory bulb, groups of glomeruli with different response latencies can be identified. This is shown for recordings of olfactory receptor neuron input measured with a calcium-sensitive axon tracer and for network dynamics measured with the voltagesensitive dye RH 1838. In the somatosensory cortex, barrels responding to the stimulation of single whiskers can be automatically detected. In the visual cortex orientation columns can be extracted. In all cases artifacts due to movement, heartbeat or respiration were separated from the functional signal by sICA and could be removed from the data set. sICA following PCA is therefore a powerful technique for data compression, unbiased analysis and dissection of imaging data of population activity, collected with high spatial and temporal resolution. © 2006 Elsevier Inc. All rights reserved.

⁎ Corresponding author. WIN Group of Olfactory Dynamics, Max-PlanckInstitute for Medical Research, Jahnstrasse 29, D-69120 Heidelberg, Germany. Fax: +49 6221 486 459. E-mail address: [email protected] (H. Spors). 1 Equal contribution. Available online on ScienceDirect (www.sciencedirect.com). 1053-8119/$ - see front matter © 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2006.08.031

Introduction One of the most important goals in systems neuroscience is to understand how brain regions respond to the sensory inputs. Generally a sensory stimulus is presented and the neuronal activity (action potentials, changes in membrane potential or ion concentration) is measured. This approach is complicated by several problems. Firstly, measurements of the neuronal activity need high spatial and high temporal resolution to examine the distributed and dynamic representation of the stimulus. Data sets resulting from these measurements have often around 10,000 to 1,000,000 pixels per frame and several hundred to several thousand frames (time points). This corresponds to the size of several 100 MB up to several gigabyte. The high dimensionality makes data analysis time consuming and complicated. This leads to difficulties analyzing and interpreting the experimental results and generating mathematical models. Spatial or temporal binning reduces the dimensionality to some extend, however, fine spatial or temporal structure is lost by such averaging. It is important to keep in mind that different neuronal populations with different spatial location can have different response time courses to the same stimulus and that this can carry information about the stimulus (e.g., Spors and Grinvald, 2002; Spors et al., 2006). Secondly, spontaneous neuronal activity and stimulus-evoked activity are mixed. Thirdly the evoked neuronal activity depends on the amplitude of the spontaneous neuronal activity just prior to the stimulation (Arieli et al., 1996). Finally especially when pushing the spatial and temporal resolution to its limits, recorded data are contaminated by noise sources like heartbeat and respiration artifacts, camera noise (e.g., fixed pattern noise or dark noise) and shot noise. Signals can therefore be seen as a mixture of a varying number of signal sources including stimulus-evoked activity, spontaneous activity and the different noise sources. The ultimate data analysis approach should be able to separate these different sources contributing to the measured complex patterns. It should maintain both the original spatial and temporal resolution and not need any averaging. In this respect one of the most important principles is the decomposition of a complex spatiotemporal data set into simple spatial patterns (modes) and

J. Reidl et al. / NeuroImage 34 (2007) 94–108

corresponding time courses. This principle is widely used to reduce the dimensionality of data sets in order to visualize and analyze complicated data sets and to separate the different signal components and noise sources (Stetter et al., 2000; Vanzetta et al., 2004). Examples of this decomposition in neuroscience are the singular value decomposition or principal component analysis (PCA) and its derivates (Karhunen, 1946; Kantz and Schreiber, 1997; Kantz et al., 1998) of MEGs (Fuchs et al., 1992) and optical imaging data (Mitra and Pesaran, 1999; Stetter et al., 2001; Vanzetta et al., 2004). Although PCA optimally reduces dimensionality of a data set, the obtained modes do not directly correspond to single (biological) signal sources (Brown et al., 2001; Roberts and Everson, 2001) but define an orthogonal basis of modes with maximal contributions to the total signal. Algorithms for blind source separation like the independent component analysis (ICA) focus on the extraction of the original sources by investigation of their statistical independency. These algorithms have already successfully been used for analyzing data of EEG (Makeig et al., 1996), MEG (Kelso et al., 1998; Barros et al., 2000), fMRI (McKeown et al., 1998a,b; Formisano et al., 2004; Aragri et al., 2006), imaging of intrinsic optical signals (Schiessl et al., 2000), and optical single cell measurements as reviewed in Brown et al. (2001) or voltage-sensitive dye measurements with low spatial resolution (Maeda et al., 2001; Inagaki et al., 2003). The data obtained with the abovementioned population recording methods either lack high spatial or high temporal resolution. Thus, so far ICA has not been tested for imaging data of population activity with both high spatial and temporal resolution. Here we examine the use of these modern statistical data analysis tools using population imaging data obtained in three different brain regions. We demonstrate for these three different brain areas that spatial ICA (sICA) analysis is able to identify several functional groups as different sources based on their different response characteristics. Furthermore, spontaneous neuronal activity, systematic artifacts (e.g., blood vessel artifacts) and noise sources can be separated by the sICA. This separation can be performed without a priori knowledge of the response characteristic, or subjective choice of filtering parameters. Data averaging and filtering are not necessary however they can be used for data preprocessing. The performance of the sICA analysis is demonstrated for data obtained with calciumsensitive and voltage-sensitive dyes in the mouse olfactory bulb, mouse somatosensory cortex, and the monkey visual cortex. Materials and methods Animal preparation, staining, and data collection Population responses to odor stimuli were recorded in the dorsal olfactory bulb of 10 C57Bl6 mice using calcium-sensitive dye imaging of the receptor neuron input (Wachowiak and Cohen, 2001) or using voltage-sensitive dye (RH 1838) imaging of the electrical network activity (Spors and Grinvald, 2002). Responses to whisker deflection were measured in the primary somatosensory cortex of 2 C57Bl6 mice using voltage-sensitive dye imaging (RH 1691, Derdikman et al., 2003). Responses to horizontal and vertical moving gratings were recorded in 2 Macaca fascicularis monkeys using voltage-sensitive dye imaging (RH 1838 or RH 1916). The surgical procedure has been reported in detail previously (Slovin et al., 2002). All animal care and procedures were in accordance with the animal ethics guidelines of the Max Planck Society and the Weizmann Institute of Science. Procedures for imaging olfactory

95

bulb activity based on calcium-sensitive dyes (Wachowiak and Cohen, 2001) and on voltage-sensitive dyes (Spors and Grinvald, 2002) have been recently described in detail elsewhere (Shoham et al., 1999; Grinvald et al., 1999). Here only aspects relevant to the present study are discussed. Surgery Mice were anesthetized using urethane (1.5 g/kg i.p.) or pentobartial (Narcuren, 50 mg/kg, i.p.). Anesthetic was supplemented throughout the experiments. All incision lines and pressure points were injected with a local anesthetic. The body temperature was kept between 36.5°C and 37.5°C using a heating pad and a rectal probe. Animals were breathing freely. The respiration of the animals was recorded using a piezzo electric strip (WPI, Sarasota). The ECG was monitored throughout the experiment and recorded during data acquisition. Calcium-sensitive dye imaging was performed through the thinned bone. For voltage-sensitive dye imaging a craniotomy over one olfactory bulb (2 × 4 mm), separately over both olfactory bulbs, or over the barrel cortex (4 × 4 mm) was made. The dura was removed to facilitate staining with the voltage-sensitive dye (RH 1838 for the olfactory bulb and RH 1691 for the barrel cortex). The tissue was stained by bath application for 60–90 min. After staining and extensive rinsing the imaging chamber was filled with agar (1.5%) and covered with a glass coverslip. Opening the cisterna magna in combination with the closed cranial window reduced pulsation of the brain due to heartbeat and respiration. Stimulation Odors were presented using a custom build olfactometer (Spors and Grinvald, 2002). Odors were purchased from Sigma-Aldrich or Fluka Chemie (Steinheim, Germany). For the somatosensory experiments single whiskers were moved by a piezzo electric stimulator (amplitude > 10°). In the monkeys experiments the visual stimulus used was a patch of square drifting gratings in two orientations 0° and 90° (contrast 100%, size 6 × 6° visual angle (DVA), spatial frequency– 2 cycles per DVA, temporal frequency 2 DVA/second, i.e., optimal stimulus parameters to determine orientation columns). The stimuli were presented on a CRT monitor (mean luminance ~50 cd/m2). The trials started with a tight fixation of the monkey on a fixation spot within a 2 × 2 DVA window. The monkeys were trained to keep a tight fixation during the trial. 200 ms after the starting of the trial the stimulus appeared on the screen for duration of 200 ms. Imaging Imaging data were acquired at 14 to 200 Hz frame rate (128 × 128 pixels) using a modified Fuji Deltaron HR 1700 camera with the Dyedaq software package (Shoham et al., 1999) and a tandem photo lens system (n.a. = 0.46). The tissue was illuminated with a halogen lamp powered by a stabilized power source (JQE 15-25 M, Kepco). Filter sets were HQ480/40 (exciter), LP505 (dichroic), and HQ 535/50 (emitter) for calcium imaging and HQ 630/20 (exciter), LP650 (dichroic), and HQ665LP (emitter) for voltage-sensitive dye imaging. This camera system provided excellent signal to noise ratio as long as the light levels allowed to nearly saturate the camera. At higher frame rates (200 Hz) the fixed pattern noise of the cmos-chip and the offsets between individual (horizontal) amplifiers increased noticeably. Data preprocessing Data were analyzed using custom made software written in the Matlab (www.mathworks.com, versions 6.5 and 7) environment.

96

J. Reidl et al. / NeuroImage 34 (2007) 94–108

Before submitting data to the sICA analysis standard routines for optical imaging data were applied: Briefly, data were divided by the average of 5 to 20 frames before stimulus onset to convert raw data into dF/F values and to minimize common noise. Regions with very low resting fluorescence or underlying bone were excluded from the analysis. After selecting the region of interest the data set had about 10,000 pixels. To try if data preprocessing improves the source separation, in some instances, bleaching and heartbeat were removed prior to the sICA. To remove the bleaching, for each pixel the average of the non-stimulated traces was low-pass filtered in time. The resulting bleaching traces were subtracted from all data that included responses to stimuli. Data analysis using PCA and sICA Data properties and data format The spatial–temporal imaging data show the measured intensity at each camera pixel rk with k = {1,…, Np}. Np is the number of recorded pixels of the region of interest at each time frame tm with m = {1,… Nf}. Nf is the number of frames recorded. This can be described in a more compact manner as matrix elements Skm = S(rk, tm). The data can be represented as a superposition of so-called modes Vi(rk) and their corresponding timedependent amplitudes ui(tm) X

ui ðtm ÞVi ðrk Þ ¼ Sðrk ;tm Þ

i

Principal component analysis (PCA, Karhunen–Loève expansion) To approximate the total spatiotemporal signal with as few modes as possible, only those with the highest contributions to the overall signal are taken into account. This is done by subtracting the spatial mean of each frame S V¼ Sðrk ;tm Þ−

X

Sðrk ;tm Þ

and minimizing the sum over the squared deviations for all time frames and pixels: XX k

S Vðrk ;tm Þ−

m

N X

Spatial independent component analysis (sICA) Based on the assumption that biological functional units do not move on the surface of the brain (i.e., they are spatially stationary structures), sICA determines modes with maximal spatial structure by using entropy measures. This assumption is fulfilled for our data sets since it has been shown using single unit (Hubel and Wiesel, 1968; Woolsey and Van der Loos, 1970) and population recordings (Orbach et al., 1985; Bonhoeffer and Grinvald, 1991; Kleinfeld and Delaney, 1996; Derdikman et al., 2003) that neurons with similar functional properties and receptive fields are clustered in space. In this context, the spatial structure is the difference of the intensity distribution of the spatial pixels of the actual modes compared to a Gaussian distribution (Roberts and Everson, 2001). In sICA the entropy of the spatial structure is minimized to extract single spatial structures by finding suitable linear combinations Aj = Bji V˜ i of the modes V˜ i previously calculated by PCA. In our case the number of PCA modes (upper limit for index i) and the number of ICA modes (upper limit for index j) are chosen to be equal. The various ICA approaches differ in the estimators for the information of the spatial structures. FastICA (Hyvarinen and Oja, 1997, available at http://www.cis.hut.fi/projects/ica/fastica/) maximizes the kurtosis, RADICAL (Learned-Miller and Fisher, 2003) has implemented direct estimators, like spacing estimates of the entropy. Further approaches make use of kernel generalized canonical correlation analysis (Kernel ICAs, Bach and Jordan, 2002). The ICA algorithm implemented in our analysis to further process the whitened spatial modes V˜ i is JADER (Cardoso, 1999; available at http://sig.enst.fr/~cardoso/stuff.html). This results in the decomposing matrix B to calculate the sICA modes A with f A ¼ BdV

rk

E¼

(see, e.g., Roberts and Everson, 2001) and is a requirement for the JADER ICA algorithm.

!2 ui ðtm ÞVi ðrk Þ

i¼1

This leads to an eigenvalue problem typically solved numerically by singular value decomposition (see, e.g., Numerical Recipes, Press et al., 1992). In most cases N can be chosen much smaller than the number of pixels Np or time frames Nf to approximate the given data. For our data sets only 5 to 20 modes are necessary to capture more than 80% of the original variance (Fig. 2). Subsequently to PCA for all numbers N of approximating modes the spatial modes Vi(rk) are normalized by dividing by the square-root of their original variances Vi ðrk Þ f Vi ðrk Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P 2 rk Vi ðrk Þ This whole procedure – subtracting the mean, orthogonal decomposition as done, e.g., by PCA, and normalization of the variances to the value of one – is often referred to as ‘whitening’

The above described approach was chosen instead of the more common inspection of the independency of the time courses ui(tm), i.e., the structure in the spatial dimension was exploited keeping the original number of pixels. An advantage of the investigation of the spatial modes instead of their temporal properties is the up to 100fold larger number of pixels compared to the number of time frames in the measurements. This allows for better estimates of cumulants or other estimates of their inherent information. Another feature which ensures good separability using sICA is the non-Gaussian distribution of the pixel intensities for groups of glomeruli (e.g., kurtosis is higher in the spatial dimension compared to the temporal. Dimension, e.g., for calcium imaging data: 8.92 ± 0.43 vs. 2.46 ± 0.02, mean ± SEM, spatial vs. temporal kurtosis). Creation of the artificial data set In order to examine the performance of the sICA in separating known sources and to test the robustness against noise we created artificial data sets. Two sources with glomerulus-like round structures and two sources with blood vessel-like elongated spatial structures were mixed. The ‘glomerulus sources’ were activated at ‘respiration rate’ with shifted phase and occupied only a part of the movie. The ‘blood vessel sources’ were mixed in at two different frequencies, one with a frequency corresponding to the heart rate and one with a frequency corresponding to the respiration rate (Fig.

J. Reidl et al. / NeuroImage 34 (2007) 94–108

2). After mixing the sources normal distributed noise was added with varying standard deviations in order to compare the performance of sICA for different signal-to-noise ratios. Scaling and sorting of the modes and time series sICA analysis can separate a high number of modes that contain stimulus responses, spontaneous activity, and artefacts to varying extend. The amplitudes of the calculated time courses and the corresponding spatial modes are only meaningful in their combination since the total amplitude in space and time is given as their product. A common way to normalize them is to set all the variances of the spatial modes to one and to take the total variances (the total spatiotemporal intensity) for the time courses, but within such a scaling the overall magnitudes of the single glomerular activities are difficult to compare. To facilitate the interpretation of the results and to obtain a meaningful scale for the comparison of different glomerular responses and contributions of artefacts we chose a different normalization procedure: Spatial modes were normalized by the absolute maximum value of the corresponding time course and vice versa the time courses were normalized by the absolute maximum of the spatial modes. Furthermore, modes and time courses were sorted according to their stimulus dependency. Stimulus dependency was calculated for each mode separately. The time series were cut into two parts: a first part before the stimulus application (‘baseline’) and a second part containing the time window of the stimulus application (‘stimulus response’). For both parts the variance was calculated separately. The ratio of the ‘stimulus response’ variance and the ‘baseline’ variance was taken as a measure of stimulus dependency. Results Quality of separation of different sources Individual functional units, i.e., different glomeruli in the olfactory system can respond with different time courses to the same stimulus and the same glomerulus can respond with different time courses to different stimuli (Spors and Grinvald, 2002; Spors et al., 2006). A single glomerulus responds as one unit as has recently been shown by 2 photon laser scanning microscopy (Wachowiak et al., 2004). Individual glomeruli can thus be seen as different sources contributing to a complex spatiotemporal pattern. To test if blind source separation can separate two groups of glomeruli and additional sources (e.g., blood vessel and respiration artifacts), we created an artificial data set (see Materials and methods) consisting of two sets of glomeruli. These glomeruli are modeled to respond with different phase relationships relative to an assumed breathing trace during the time interval of stimulus presentation. Artifacts at heart rate and respiration frequency were added to the entire series of movie frames. The spatial distribution of the sources and their time courses are shown in Fig. 1A. After mixing of the sources, normal distributed noise was added (see movie frames in Figs. 1B1 and B2). PCA analysis already results in a good separation of the time series and the spatial modes, that either contain only glomerular or only ‘blood vessel’ patterns (Fig. 1C). However, the same spatial structures (e.g., the same two groups of glomeruli) appear twice as varying superposition in different spatial modes and the original time courses cannot be recovered (Fig. 1C). Additional consideration of higher statistical moments using sICA allows for a better source separation. sICA analysis results in spatial modes that are more similar to the

97

original patterns (compare Figs. 1A and D) than the modes recovered with PCA. To determine the influence of noise on the identification of the initial sources we plotted the spatiotemporal correlation of the sources recovered by PCA and sICA with the original sources as a function of noise level (Fig. 1E1). The noise level is taken as the ratio of the standard deviations (SD) of noise and signal. Fig. 1E1 shows that even at high noise levels up to six times of the ‘signal’ SD the original sources and their time courses can be recovered with high fidelity using sICA. Conversely, PCA modes have only a small correlation coefficient with the original sources almost independent of the noise level. These findings suggest that PCA is able to separate odor responses from noise in certain conditions and that sICA in addition can separate the different glomerular groups. In the displayed case the spatially overlapping blood vessel source (source 4) and a glomerulus (from source 1) could be separated by PCA and by sICA. Using additional examples of simulated data we tested if overlapping sets of glomeruli can be separated by sICA. This is indeed the case as long as the sources of glomeruli are activated with different time courses (data not shown). At the same time the robustness against random noise is reduced by the overlap of the glomeruli. As dimensionality reduction using PCA and sICA does not preserve the entire variance of a given data set we calculated how much of the original variance is maintained after PCA and sICA analysis as a function of the noise level for the 4 leading modes (Fig. 1E2). Dimensionality reduction Principal component analysis (PCA) can substantially reduce the dimensionality of complex data sets (Karhunen, 1946; see, e.g., Kantz and Schreiber, 1997; Kantz et al., 1998) while preserving most of the non-random fine details in the measured data (Fig. 2A). To test how much of the spatiotemporal variance is retained from our imaging data we submitted 38 representative data sets from 12 animals (6 for calcium-sensitive dye imaging (CaSDI), 6 for voltage-sensitive dye imaging (VSDI)) to a standard PCA algorithm (see Materials and methods). 22 data sets (calciumsensitive dye imaging in the olfactory bulb) consisted of 4 × 56 frames and ~ 10000 pixels each. 16 data sets (voltage-sensitive dye imaging in the olfactory bulb) consist of 4 × 156 frames and ~ 10000 pixels each. 16 data sets (voltage-sensitive dye imaging in the somatosensory cortex) consist of 4 × 74 frames and ~ 10000 pixels each. Fig. 2A shows one representative example of fluorescence changes in response to odor presentation measured with calcium imaging. The approximation with 6 modes (2.7% of original dimensionality) is very similar to the original data. The difference between the original data set and the data after dimensionality reduction contains only negligible stimulus-dependent glomerular responses (lower row in Fig. 2A). In this example 80.3% of the original variance was kept by the 6 modes and their time courses. In the other data sets 80% of the total variance was kept with 6 modes for CaSDI and 20 modes for VSDI, corresponding to a data compression to 2.6% and 3.2% respectively (Fig. 2B, blue traces). The lost variance could be attributed to noise and artifacts or to fine details of the stimulus response. We therefore asked next how much of the stimulusspecific response and how much of the non-stimulus-related signals are retained after dimensionality reduction. The amount of original variance that was still contained in the data set after dimensionality reduction using PCA was calculated separately for epochs with stimulus and epochs without stimulus. Dimensionality

98

J. Reidl et al. / NeuroImage 34 (2007) 94–108

Fig. 1. sICA separates signal components better than PCA—artificial data. (A) Artificial spatial patterns with corresponding time courses (sources). The first two sources simulate two groups of ‘glomeruli’ with different phase relationship to the respiration (phase shift: see inset). The third and fourth sources simulate ‘respiration’ and ‘heartbeat’ artifact. They are present during the entire movie. Normal distributed noise was added to all movie frames. (B) Artificial spatiotemporal data produced by linear combination of the sources shown in panel (A). (B1) The movie frames without noise at the time when the sources 1 and 2 are introduced (see arrow in A) show the glomeruli like signal and the artifacts. Panel (B2) same as panel (B1) after adding the noise (arrow in E1 for noise level). (B3) Combination of the two sICA modes corresponding best to the glomerulus like sources. This results in a cleaned signal. (B4) Difference between spatiotemporal pattern with noise (B2) and the cleaned signal in panel (B3). (C) Modes and time courses obtained by PCA analysis. The modes are mixtures of the sources. (D) Modes and time courses obtained by sICA analysis after PCA. The patterns are highly correlated with the sources. The correlation is better for the intermittently active sources 1 and 2. (E) Dependence of the source recovery on the noise level. The noise level is given as the ratio of standard deviation of the noise and the standard deviation of the signal. (E1) Spatiotemporal correlation of the first and the second source and the corresponding mode recovered by PCA (gray) and sICA (black) as a function of the noise level. (E2) Proportion of the total variance maintained in the leading 5 PCA modes depending on the noise level.

reduction kept most of the original variance for the stimulus epochs (> 80% for reduction to 2,6% of dimensionality, Fig. 2B, green traces) while more of the original variance was lost in traces without stimulus (> 40% for reduction to 2,6% of dimensionality, Fig. 2B, red traces). Epochs with stimulus also contain the contributions of noise and artifacts. Assuming independency of these different contributions the variance for stimulus epochs is the

sum of the variances of stimulus-related signal and control signal (containing only spontaneous activity, artifacts and noise). The small decrease of the variance of stimulus epochs with reduction of dimensionality is due to the reduction of the noise contribution, while the stimulus-dependent signal is almost completely retained in the first 6 PCA modes. This demonstrates that dimensionality reduction with PCA is possible for data sets obtained with voltage-

J. Reidl et al. / NeuroImage 34 (2007) 94–108

99

Fig. 2. Dimensionality of measured odor representations can be reduced. (A) PCA and sICA maintain fine spatiotemporal features of odor responses. Top, input to the olfactory bulb measured with the calcium-sensitive axon tracer calcium green dextran. Left, resting fluorescence. Each movie frame is the integral of 77 ms. Stimulation at t = 0 with ethyl butyrate (1% liquid dilution). Middle, approximation by sum of the 6 leading modes (of 224 total modes) calculated using PCA and sICA (Jade algorithm). 80.3% of the original variance is kept. Bottom, difference between the original signal and the sum of the 6 leading modes. Scale bar 0.2 mm. (B) Proportion of the original variance depending on the number of leading modes recombined. (B1) Data were obtained with the calcium-sensitive axon tracer Calcium Green Dextran. Calculated block wise and subsequently averaged for 22 blocks (from 6 mice). Each block consists of 2 epochs (2 × 56 frames) with stimulus (ethyl butyrate 1%) and 2 epochs without stimulus. (B2) Data obtained with the voltage-sensitive dye RH1838. Calculated for 16 blocks (from 4 mice). Each block consists of 2 epochs (2 × 156 frames) with stimulation (ethyl butyrate 1%) and 2 control traces (2 × 156 frames). (B3) Data obtained with the voltage-sensitive dye RH1691. Calculated for 16 blocks (from 1 mouse). Each block consists of 2 epochs (2 × 71 frames) with whisker stimulation and 2 control traces (2 × 71 frames).

or calcium-sensitive dyes in the olfactory bulb and somatosensory cortex of mice. Separation of functional groups and vascular patterns in single trials As a first test of PCA and sICA analysis on experimental spatiotemporal data sets we used a single trial (56 frames) recorded with a calcium-sensitive axon tracer (Calcium Green dextran; Spors et al., 2006). In the main olfactory bulb axon terminals of olfactory

receptor neurons were selectively stained. Fluorescence changes in freely breathing mice were recorded at 14 Hz frame rate during stimulation with an odor (Ethyl butyrate 1%). Individual movie frames of this data set are displayed in Fig. 3A. PCA analysis of the response to a single stimulus presentation results in a few leading spatial modes with high variances containing most of the glomerular responses (arrows) and blood vessel artifacts (arrow heads, Fig. 3B1). The same spatial structures (i.e., groups of glomeruli) appear in several of the spatial PCA modes. The corresponding time courses (Fig. 3C1) are modulated with the respiration as expected for an

100

J. Reidl et al. / NeuroImage 34 (2007) 94–108

Fig. 3. Separation of glomerular groups and blood vessel artifact by sICA analysis of odor-evoked spatiotemporal patterns. (A) Single movie frames (frame integral 77 ms) at time points indicated in ms relative to the beginning of stimulation with ethyl butyrate (1% liquid dilution) recorded using the calcium-sensitive tracer calcium green dextran. These complex spatiotemporal patterns represent the input to the olfactory bulb as the dye has been selectively loaded into the terminals of the olfactory receptor neurons. Left panel: resting fluorescence and blood vessel pattern. Scale bar 0.2 mm, Color scale ΔF/F − 0.5% to 3%. (B1) Leading four PCA modes contain mixtures of the caudo-lateral and anterior medial group of glomeruli (arrows) and the blood vessel pattern (arrow heads). (B2) sICA separates groups of glomeruli (arrows) in two different spatial modes. A third mode contains parts of the blood vessel pattern (arrow heads, compare (A) left panel). The forth sICA mode shown here contains a weakly activated glomerulus (lower right) and an artifact (top left). sICA results in significantly better separation of glomerular groups. 4 out of 20 calculated spatial modes are displayed. (C) The time series of the spatial modes. (C1) The time courses of the principal components (PC) 1 and 3 follow the respiration. The time courses of PC1 and PC2 are both modulated at the same frequency as the recorded heartbeat. (C2) The time courses of the independent components (IC) 1 and 3 follow the respiration with different phase relationship and varying degree of modulation. Only IC 2 is modulated at the same frequency as the recorded heartbeat. Single heartbeats were missed because of the limited sampling rate of 14 Hz.

olfactory response. The degree of modulation varies in the different time courses (compare red and blue time courses). The time course of the strongest principal component appears to be the mixture of at least two time courses. The first one is odor-specific and modulated with the respiration frequency. The second one is a fast modulation at the heartbeat frequency. The blood vessel pattern itself is not visible in the corresponding spatial mode ‘PC1’ due to the overlapping glomerular pattern.

Using sICA (Jade algorithm, see Materials and methods) glomeruli and the blood vessel pattern can be further separated into different spatial modes (Fig. 3B2). Here glomeruli or groups of glomeruli with different temporal response patterns are separated into different spatial modes. Their time courses still show the prominent respiration but no heartbeat modulation (IC 1 and 3 in Fig. 3C2). The heartbeat artifact is now mainly concentrated in one mode with a time course that is not stimulus-dependent (Figs. 3B2 and C2:

J. Reidl et al. / NeuroImage 34 (2007) 94–108

IC 2). The response of the two groups of glomeruli (IC 1 and 3) is shifted in time and the degree of their modulation with respiration varies. These different temporal responses can now be clearly attributed to the groups of glomeruli shown in the spatial modes IC1 and IC3 in Fig. 3B2, as the same glomeruli no longer appear in multiple modes. Thus, manual selection and comparison of the time courses of many individual glomeruli (standard data analysis), are not any more necessary in order to find groups of glomeruli with common response characteristics. Reproducibility of the sICA analysis To test the reproducibility of the separation into functional groups we repeated the sICA analysis on individual trials with the presentation of the same odor stimulus. This approach depends on the high reproducibility of the response to the same stimulus. Both the spatial maps (Fig. 4A) and the phase shifts between the time courses of activation (Fig. 4B) are very similar across trials. Small difference between the responses, here the ratio between the response in the first and the second inspiration, is kept by sICA analysis. We test the reproducibility of the spatial glomerular patterns in 6 animals with 8 to 10 repetitions of the identical stimulus (Fig. 4C). The two modes with the largest ratio of stimulus-related variance divided by non-stimulus-related variance were selected and compared across repetitions. Fig. 4C shows the reproducibility of the two main modes across 4 out of 8 repetitions of the same stimulus (ethyl butyrate 1%, different animal than Fig. 4A). For each mode the average across the 8 repetitions was calculated. The difference between the modes calculated from the individual repetitions and the average mode is shown in the 2nd row of Fig. 4C. This analysis was repeated on data from 6 mice. The average correlations between corresponding spatial modes were 0.85 ± 0.01 (mean ± SEM) and 0.77 ± 0.02, for the first and second mode respectively. Spatial lowpass filtering (Gauss filter, σ = 13.3 μm) increased the correlation coefficients only by 0.037 ± 0.001 (mean ± SEM) indicating that most high-frequency spatial noise was already removed by the sICA analysis. Thus separation into functional groups is highly reproducible using sICA. For the validity of the sICA analysis it is important to show that functional neuronal groups – here glomeruli – identified in the sICA modes are indeed activated by the stimulus and that glomeruli that are clustered in the same modes have very similar time courses. We compared the time courses of the sICA modes with the raw fluorescence time courses after conversion to dF/F and bleaching correction. The averaged fluorescence change of manually chosen glomeruli was very similar to the one of the mode, that contained them (Fig. 4D). The average correlation to the same mode (r = 0.74 ± 0.04, mean ± SEM, n = 6 experiments with 1 to 4 repetitions each) was significantly higher than the correlation to other modes (r = 0.54 ± 0.04, p ≪ 0.001, Wilcoxon test). Furthermore, time courses of glomeruli that are contained in the same mode are almost identical (r = 0.91 ± 0.01, 3 glomeruli per mode, n = 6 experiments with 1 to 4 repetitions each) while time courses of glomeruli contained in different modes are still similar but not as highly correlated (r = 0.76 ± 01, p ≪ 0.001, Wilcoxon test). sICA identifies corresponding glomerular groups in VSDI and CaSDI Identification of functional subgroups with blind source separation methods should ideally work for different imaging methods. After establishing the effectiveness of sICA in combina-

101

tion with calcium imaging data, we tested data recorded with the voltage-sensitive dye RH 1838. This dye mainly reports the electrical activity in the network of the olfactory bulb (Spors and Grinvald, 2002). The same strain of mice and the same odorant were used as with the calcium imaging. Analysis of voltagesensitive dye imaging data was expected to provide additional difficulty for the analysis, because of the smaller fractional change, and because of the less localized response (Fig. 5A, upper left panel) (Spors and Grinvald, 2002). Nonetheless glomeruli, or activated regions, and blood vessel artifacts could be separated (Fig. 5A, IC1–IC3). The time courses of the two modes (IC1 and IC2) contained a reproducible response after stimulus application. One mode responded tonically, i.e., with weak modulation at respiration frequency, the second mode contained a full response with each inspiration cycle and had a shorter latency than the time course of the first mode (IC1). Thus sICA analysis works also well for voltage-sensitive dye imaging data sets. The quantitative comparison with the calcium imaging data (see above, Figs. 3 and 4) is complicated because of the difference in anesthetic used (urethane for voltage-sensitive dye imaging, Narcuren for calcium imaging) and the resulting difference in respiration rate. However the coarse localization of the glomeruli with fast response kinetics and strong respiration modulation was comparable. This is consistent with the hypothesis that a substantial part of the electrical network dynamics measured with the voltage-sensitive dye are determined by the spatiotemporal dynamics to the olfactory bulb which can be measured with the calcium-sensitive dye. sICA analysis for other brain regions The ideal blind source separation method for fast high-resolution imaging should be applicable to different brain areas. We tested sICA analysis on data sets of two other brain regions. In the following we illustrate sICA results from the analysis of whisker responses in the somatosensory cortex and responses to orientated gratings in the visual cortex of awake behaving monkeys. In both cases sICA separates biological signals, biological noise and camera noise. Somatosensory cortex The whisker barrel field of the somatosensory cortex (Woolsey and Van der Loos, 1970) is not only a well suited structure to test new imaging techniques (Shoham et al. 1999, Murakami et al., 2004) but also for testing data analysis. The anatomical structure of the socalled barrel cortex is well known and the responses of individual neurons (Brecht and Sakmann, 2002; Brecht et al., 2003; Manns et al., 2004), as well as populations of neurons (Kleinfeld and Delaney, 1996; Derdikman et al., 2003) are well characterized. We therefore tested sICA analysis on single whisker responses measured with a voltage-sensitive dye in the mouse. The electrical response to stimulation of a single whisker is fast and spreads rapidly to the neighboring columns (Fig. 6A). When stimulating different whiskers, the centers of the electrically activated areas differ. However the total activated area shows considerable overlap (Fig. 6A, time point 25 ms). For sICA analysis we combined two data sets each containing 16 stimulations of one whisker (stimulated whiskers were C1 and D2). Four out of 30 calculated modes are shown in Fig. 6B. Among these are two spatial modes whose areas of activation hardly overlap and whose time courses are stimulus-specific. Furthermore, a mode containing a blood vessels pattern could be identified. One example of a mode containing a regular spatial

102

J. Reidl et al. / NeuroImage 34 (2007) 94–108

J. Reidl et al. / NeuroImage 34 (2007) 94–108

103

Fig. 5. sICA separates neuronal subgroups in the network of the mouse olfactory bulb. (A) Electrical activity of olfactory bulb neurons measured with the voltage-sensitive dye RH 1838 in a mouse stimulated with the odor ethyl butyrate, urethane anesthesia. Left, map averaged over the first respiration cycle (clipping 0 to 0.1%) and green image. Middle, three out of 20 calculated modes are displayed. (B) The anterior and more medial activated region (IC1) responds with longer latency and modulates less with the respiration than the more posterior lateral one (IC2). During the odor response the time courses of the first two modes are shifted by almost half a respiration cycle (compare respiration trace). A third mode correlates well with the blood vessel pattern and is modulated at the heart rate. The slow increase during the odor stimulation resembles the intrinsic blood vessel response. Clipping − 0.05% to 0.17%, scale bar 0.2 mm.

pattern – so-called fixed pattern – is shown in Fig. 6B under the figure heading “pattern noise”. This pattern is specific to cameras with a Cmos detector like the Fuji Deltaron HR 1700, used here. The time courses of the two modes 1 and 2 that show a barrel-like outline follow the stimulation (Fig. 6C). The time courses of mode 1 (blue traces) and mode 2 (red trace) have different amplitudes and different time courses when stimulating whiskers C1 and D2. Mode 1 shows the strongest response of all sICA modes to stimulation of the whisker C1. We conclude that the whisker C1 is the principal whisker of the barrel lighting up in mode 1. In analogy whisker D2 is the principal whisker of the barrel lighting up in mode 2. As known from electrophysiological recordings (Brecht et al., 2003; Manns et

al., 2004) the responses to the surround whiskers show smaller amplitudes, have a longer latencies and delayed peaks (Fig. 6C). Visual cortex of the awake behaving monkey One of the most interesting and challenging preparations for fast in vivo population imaging is to record from trained awake behaving animals. Electrical responses to horizontal and vertical moving gratings were recorded in two awake monkeys in area 17 and 18 of the visual cortex. The standard differential map obtained by division of responses to the two gratings was used as a reference pattern (Fig. 7A). Before entering the data into the sICA 14 responses to a horizontal grating and 14 responses to a vertical

Fig. 4. Separation into functional groups of glomeruli by sICA is reproducible. (A and B) Display of the spatial modes and the corresponding time courses analogous to Fig. 3. Here sICA analysis was repeated on two individual responses to stimulation with ethyl butyrate 2s (1% liquid dilution) in the same animal. (A) The selected 4 spatial modes (out of 20 calculated spatial modes) contain three groups of glomeruli and a blood vessel artifact. The grouping of the activated glomeruli is highly reproducible across the repetitions. Correlation of these patterns with the average of all corresponding modes across all repetitions is 0.90 ± 0.02 (mode 1), 0.92 ± 0.01 (mode 2), 0.62 ± 0.18 (mode 3), 0.8 ± 0.05 (mode 4, mean ± SD, n = 10 repetitions). Scale bar 0.2 mm. (B) The relative latency, phase shift, and degree of modulation of the glomerular groups are preserved across the repetitions of the stimulation and analysis. (C) sICA was calculated separately on 8 responses to the same odor. (C1) Four repetitions (out of 8) of the first sICA mode. The difference between the sICA modes calculated from a single repetition to the average of all repetitions is shown below. Average correlation of the single mode with the average of all modes is 0.86 ± 0.01 (mean ± SD, n = 8 repetitions). (C2) Comparison of the second sICA mode. Average correlation of the single mode with the average of all modes is 0.83 ± 0.03 (mean ± SD, n = 8 repetitions). Color bar as in panel (A). (D) Comparison of sICA time courses and time course of glomeruli of bleaching corrected raw data. (D1) The time course of mode1 (black trace, also see (A) and (B)) was compared to the time course of the fluorescence change of two caudal glomeruli that appear in mode1 (see (A) and inset). (D2) The time course of mode2 (black trace) was compared to the time course of two anterior medial glomeruli that appear in mode2 (see (A) and inset).

104

J. Reidl et al. / NeuroImage 34 (2007) 94–108

Fig. 6. sICA separates functional domains in the somatosensory cortex. (A) Movie frames of evoked electrical responses in the barrel field of the somatosensory cortex after stimulation of two different single whiskers. The depolarization starts localized in the barrel corresponding to the moved whisker and spreads over a wide region within the next 10 ms (following two movie frames). Time in ms relative to stimulus onset. The top and bottom rows show the response to the stimulation of the C1 and D2 whisker, respectively. Scale bar 1 mm. (B) sICA analysis of data containing both responses (see Materials and methods) recovers modes (2 out of 30) in which the two activated whisker barrels show up without spread (Mode 1 and Mode 2). Examples of further modes containing the blood vessel artifact and the camera pattern noise are given below. The heartbeat mode contains major blood vessels (arrows), compare blood vessel pattern on left side. (C) The time courses of the two modes in response to the stimulus are shown for the average of 8 repetitions of each stimulus. Both modes (whisker barrels) respond to both stimuli, however the amplitude is bigger and the latency is shorter if the principal whisker is stimulated.

J. Reidl et al. / NeuroImage 34 (2007) 94–108

105

Fig. 7. sICA recovers voltage-sensitive dye mapping signal during imaging in the awake monkey. Responses to horizontal and vertical moving gratings were recorded in the visual cortex of awake monkeys with the voltage-sensitive dye. 14 repetitions of horizontal and vertical grating were averaged separately. 200 ms before to 200 ms after the response onset were combined for both responses. (A) The differential map calculated by dividing responses to the vertical grating by the responses to the horizontal grating during the first 150 ms of the response. (B) The blood vessel pattern recorded with green illumination. (C) 4 out of 10 sICA modes are shown. In the mode containing the mapping signals dark patches (neuronal groups mainly responding to the horizontal grating) are labeled with green stars to allow the comparison to the differential map. Green stars in the remaining three modes have no clear relationship to the dark or white patches. Blood vessels and the electrode are marked by arrows. Gauss band-pass filter σlow-pass = 31.3 μm, σhigh-pass = 13.3 μm.

grating were averaged separately. sICA of the combined data containing responses to both stimuli results in several spatial modes. One of these modes has a high similarity to the standard differential map (Fig. 7C, mapping signal). Other modes contain camera pattern noise, blood vessel and electrode artifacts (Fig. 7C). A further example from a second animal is not shown here. Thus even in the most difficult preparation for voltage-sensitive dye imaging, the awake monkey, artifacts and the functional architecture can be extracted using sICA. Noise removal using sICA Measurements of responses to weak stimuli and measurements at high frame rates are often limited by the signal to noise ratio. We therefore examined if sICA can be used to remove different types of noise and artifacts from population responses measured with high spatial and temporal resolution. sICA reproducibly identifies blood vessel artifact patterns and camera pattern noise (Figs. 3–7). As a first test of the performance of sICA for noise and artifact removal we identified modes with time courses that did not contain dynamics correlated to the known stimulus. These modes were recombined to form a noise/artifact signal, which was subtracted from the original data. In Fig. 8 the response to 16 stimulations of the whisker C1 measured with the voltagesensitive dye RH 1691 in the somatosensory cortex is displayed.

The images in the top row were only normalized to the resting fluorescence intensity. A prominent blood vessel artifact (Fig. 8, top row, arrows) has about the same amplitude as the evoked signal. Traditionally these artifacts are removed by averaging several non-stimulated traces and subtracting them. This analysis removes, e.g., the blood vessel artifact (Fig. 8, middle row). However camera pattern noise in form of diagonal stripes is still prominent. Subtracting manually selected ‘noise modes’ yields superior signal to noise ratio (Fig. 8, bottom row). This indicates that sICA analysis cannot only be used to identify neuronal groups based on their different responses to the same stimulus, but also can be exploited to remove biological and camera-related noise types. Discussion An optimal data analysis tool for in vivo imaging should be able to objectively extract biologically meaningful information out of data sets with high spatial and temporal resolution (i.e., very high dimensionality). Here we show that sICA can successfully separate several functional neuronal groups and artifacts contained in data sets from three different brain regions using either calcium- or voltage-sensitive indicators. Neuronal groups responding to the presented stimuli can be reproducibly identified using sICA and the dimensionality of the data set is reduced effectively.

106

J. Reidl et al. / NeuroImage 34 (2007) 94–108

Fig. 8. Removal of blood vessel artifacts and pattern noise by sICA. Top, response to stimulation of the C1 whisker measured with the voltage-sensitive dye RH1691 in the mouse somatosensory cortex. Time is given in ms relative to the stimulus onset. Images were normalized to baseline frames (that are not displayed here). The prominent blood vessel artifact (arrows) has a different time course than the neuronal response. Average of 8 repetitions. Middle, subtraction of the average of the none-stimulated movies. The blood vessel artifact was removed; however the pattern noise of the camera became more prominent. Bottom, after sICA without preprocessing, independent components that don’t respond to the stimulus are subtracted.

Why does sICA work well with fast imaging of population activity? At the level of neuronal populations it is possible to identify functional groups or clusters of neurons since neurons or dendritic processes of neurons with similar response characteristics are often clustered (Mountcastle, 1957; Hubel and Wiesel, 1965). In the olfactory bulb this clustering is given by the organization of the receptor neuron afferents into spherical neuropil structures called glomeruli (Ressler et al., 1994; Vassar et al., 1994). In the somatosensory cortex so-called barrels receive specific input from single whiskers (Woolsey and Van der Loos, 1970). The spatial structures respond as units that do not move in space, whereas their temporal response profiles vary with stimulus identity and stimulus strength. An essential requirement for sICA algorithms is that the underlying patterns have significant non-Gaussian distributions, e.g., high kurtosis or high skewness. This is given for our experimental data because the individual functional units with high local intensities (here the glomeruli or initial response of barrels) are small relative to the imaged background area. Due to the small size of the glomeruli the variance of the corresponding modes is small compared to the kurtosis. The synchronous responses of one or several small glomeruli result in asymmetric spatial intensity distributions with high skewness. Comparing sICA to other analysis methods Conventional data analysis requires choosing or outlining regions of interest in order to obtain time courses of activation. This is typically done manually based on the time averaged response maps. In many cases these response maps have to be filtered in space using filter parameters which are to some extend arbitrary. sICA analysis automatically yields time courses that are not based on subjective outlining of activated areas and any assumptions about the size of the underlying functional neuronal groups. Functional units that share specific stimuli or that have similar time courses are

automatically clustered without the need for laborious outlining of the individual regions of interest. Comparison of a large number of time courses can be avoided using sICA analysis. Differential imaging works very well if stimuli are orthogonal in stimulus space (e.g., comparing horizontal and vertical gratings). In the olfactory system, however, different stimuli activate overlapping sets of glomeruli based on the broad tuning spectrum of the olfactory receptor neurons—therefore, differential imaging is only of very limited use here. PCA has been successfully used to separate functional responses and artifacts for intrinsic signal imaging data (Stetter et al., 2000; Vanzetta et al., 2004). However, it is common knowledge that often the resulting components are not related to the biological functional domains. While PCA often yields good separation of stimulus responses and artifacts in some preparations (e.g., intrinsic imaging of the visual cortex of the monkey) it was not providing satisfactory identification of artifacts and biological responses in the data sets obtained from the olfactory bulb and the somatosensory cortex. Here we demonstrate that sICA performs this job better (Figs. 3–7) and can further improve the separation of noise patterns (e.g., camera noise), of various artifacts and several different neuronal subpopulations with distinct activity patterns. Temporal correlation between physiological parameters (e.g., electrocardiogram (ECG) or respiration) and stimulus timing on one side and biological noise or biological responses on the other side has been successfully used in a number of studies (Rector et al., 1999; Kalatsky and Stryker, 2003). This approach is particularly useful if the physiological ‘clock signal’ is rhythmic and regular, or if the time courses of the biological response can be easily determined and used as a model function (Rector et al., 2001). In the case of the olfactory bulb different glomeruli respond with different time courses in a stimulus-specific manner and the modulation with the physiological ‘clock signal’ the respiration varies between glomeruli (Spors et al., 2006). A blind source separation technique that works without assumptions about the temporal profile of biological artifacts and biological responses is therefore advantageous.

J. Reidl et al. / NeuroImage 34 (2007) 94–108

Limitations of the sICA analysis Conventional ICA algorithms like Fast ICA, JADE, Maximum Likelihood can only obtain a number of sources that is lower than or equal to the number of samples that were recorded (in our case the number of frames, Brown et al., 2001), however some algorithms (Cardoso, 1991; Lee et al., 1999) can recover more sources than samples. This limitation often does not come into effect since 20 modes are sufficient to capture almost completely the stimulusrelated variance and because computation time increases nonlinearly with the number of sources that the algorithm (e.g., JADE, Cardoso, 1999) searches for. The part of measured signal that is not contained in the sICA approximation (Fig. 1A) does not contain obvious responses to a stimulus. Stimulus-specific variance was reduced less compared to the baseline variance. That is, more of the noise artifacts and high-frequency noise are ignored than of the biological signal. A limitation of the current algorithms is the need to subtract the mean of each time point or the mean of each pixel in time in order to whiten the data as a first step before the sICA using the JADER algorithm. As the signals in our measurements are very small and distributed around 0 one might think that whitening of the data is not necessary. Indeed the violation of the whitening requirement results in almost identical spatial modes (r = 0.933 ± 0.008, mean ± SEM, calcium imaging, n = 10 movies∗20 modes). However due to theoretical reasons we prefer to whiten the data before entering them into the sICA algorithm even if this makes it difficult to identify nethyperpolarization as it would haven been possible if the mean over the entire frame was not subtracted. In cases where net-hyperpolarization is searched for, the sICA patterns can be used to determine putative regions of inhibition. The time courses of the averaged pixels values of these regions can subsequently be calculated from the raw data to judge if net-hyperpolarization was measured. Conclusions In summary sICA analysis allows for an objective separation of different functional neuronal groups and the simultaneous recovery of activation time courses of these functional groups. Artifacts caused by heartbeat, respiration, or camera pattern noise can be identified and subtracted from the original data. Data compression and reduction of random high-frequency noise are useful byproducts of the combined PCA and sICA analysis. Acknowledgments We thank Rina Hildesheim for synthesizing the voltage-sensitive dyes. We thank Andreas Schaefer, Thomas Kuner, Matt Wachowiak, and Rainer Friedrich for their technical help and discussions. We thank Guo Qingzhu and Thomas Künsting for their help with programming. The work was funded by the Heidelberger Akademie der Wissenschaften (WIN Program), the MPG, the BMBF and ISF. References Aragri, A., Scarabino, T., Seifritz, E., Comani, S., Cirillo, S., Tedeschi, G., Esposito, F., Di Salle, F., 2006. How does spatial extent of fMRI datasets affect independent component analysis decomposition? Hum. Brain Mapp. 27, 736–746. Arieli, A., Sterkin, A., Grinvald, A., Aertsen, A., 1996. Dynamics of ongoing activity: explanation of the large variability in evoked cortical responses. Science 273, 1868–1871.

107

Bach, F.R., Jordan, M.I., 2002. Kernel independent component analysis. J. Mach. Learn. Res. 3, 1–48. Barros, A.K., Vigario, R., Jousmaki, V., Ohnishi, N., 2000. Extraction of event-related signals from multichannel bioelectrical measurements. IEEE Trans. Biomed. Eng. 47, 583–588. Bonhoeffer, T., Grinvald, A., 1991. Iso-orientation domains in cat visual cortex are arranged in pinwheel-like patterns. Nature 353, 429–431. Brecht, M., Sakmann, B., 2002. Dynamic representation of whisker deflection by synaptic potentials in spiny stellate and pyramidal cells in the barrels and septa of layer 4 rat somatosensory cortex. J. Physiol. 543, 49–70. Brecht, M., Roth, A., Sakmann, B., 2003. Dynamic receptive fields of reconstructed pyramidal cells in layers 3 and 2 of rat somatosensory barrel cortex. J. Physiol. 553, 243–265. Brown, G.D., Yamada, S., Sejnowski, T.J., 2001. Independent component analysis at the neural cocktail party. Trends Neurosci. 24, 54–63. Cardoso, J.F., 1991. Super-symmetric decomposition of the fourth-order cumulant tensor. Blind identification of more sources than sensors. Toronto 3109–3112. Cardoso, J.F., 1999. High-order contrasts for independent component analysis. Neural Comput. 11, 157–192. Derdikman, D., Hildesheim, R., Ahissar, E., Arieli, A., Grinvald, A., 2003. Imaging spatiotemporal dynamics of surround inhibition in the barrels somatosensory cortex. J. Neurosci. 23, 3100–3105. Formisano, E., Esposito, F., Di Salle, F., Goebel, R., 2004. Cortex-based independent component analysis of fMRI time series. Magn. Reson. Imaging 22, 1493–1504. Fuchs, A., Kelso, J.A.S., Haken, H., 1992. Phase transition in the human brain: spatial mode dynamics. Int. J. Bifurcation Chaos 2, 917–939. Grinvald, A., Shoham, D., Shmuel, A., Glaser, D.E., Vanzetta, I., Shtoyerman, E., Slovin, H., Wijnbergen, C., Hildesheim, R., Sterkin, A., Arieli, A., 1999. In-vivo optical imaging of cortical architecture and dynamics. In: Windhorst, U., Johansson, H. (Eds.), Modern Techniques in Neuroscience Research. Springer, Berlin, pp. 893–969. Hubel, D.H., Wiesel, T.N., 1965. Receptive fields and functional architecture in 2 nonstriate visual areas (18 and 19) of cat. J. Neurophysiol. 28, 229–289. Hubel, D.H., Wiesel, T.N., 1968. Receptive fields and functional architecture of monkey striate cortex. J. Physiol. (Lond.) 195, 215. Hyvarinen, A., Oja, E., 1997. A fast fixed-point algorithm for independent component analysis. Neural Comput. 9, 1483–1492. Inagaki, S., Katura, T., Kawaguchi, H., Song, W.J., 2003. Isolation of neural activities from respiratory and heartbeat noises for in vivo optical recording in guinea pigs using independent component analysis. Neurosci. Lett. 352, 9–12. Kalatsky, V.A., Stryker, M.P., 2003. New paradigm for optical imaging: temporally encoded maps of intrinsic signal. Neuron 38, 529–545. Kantz, H., Schreiber, T., 1997. Nonlinear Time Series Analysis. Camebridge University Press. Kantz, H., Kurths, J., Mayer-Kress, G., 1998. Nonlinear Analysis of Physiological Data. Springer-Verlag. Karhunen, K., 1946. Zur spektralen Theorie stochastischer Prozesse. Ann. Acad. Sci. Fenn., A 1. Math. - Phys. 37. Kelso, J.A.S., Fuchs, A., Lancaster, R., Holroyd, T., Cheyne, D., Weinberg, H., 1998. Dynamic cortical activity in the human brain reveals motor equivalence. Nature 392, 814–818. Kleinfeld, D., Delaney, K.R., 1996. Distributed representation of vibrissa movement in the upper layers of somatosensory cortex revealed with voltage-sensitive dyes. J. Comp. Neurol. 375, 89–108. Learned-Miller, E., Fisher, J., 2003. ICA using spacings estimates of the entropy. J. Mach. Learn. Res. 4, 1271–1295. Lee, T.-W., Lewicki, M.S., Girolami, M., Sejnowsk, T.J., 1999. Blind source separation of more sources than mixtures using overcomplete representations. Signal Processing Letters. IEEE 6, 87–90. Maeda, S., Inagaki, S., Kawaguchi, H., Song, W.J., 2001. Separation of signal and noise from in vivo optical recording in Guinea pigs using independent component analysis. Neurosci. Lett. 302, 137–140.

108

J. Reidl et al. / NeuroImage 34 (2007) 94–108

Makeig, S., Bell, A.J., Jung, T.-P., Sejnowski, T.J., 1996. Independent component analysis of electroencephalographic data. In: Touretzky, D., Mozer, M., Hasselmo, M. (Eds.), Advances in Neural Information Processing Systems, vol. 8. MIT Press, pp. 145–151. Manns, I.D., Sakmann, B., Brecht, M., 2004. Sub- and suprathreshold receptive field properties of pyramidal neurones in layers 5A and 5B of rat somatosensory barrel cortex. J. Physiol. Online 556, 601–622. McKeown, M.J., Makeig, S., Brown, G.G., Jung, T.P., Kindermann, S.S., Bell, A.J., Sejnowski, T.J., 1998a. Analysis of fMRI data by blind separation into independent spatial components. Hum. Brain Mapp. 6, 160–188. McKeown, M.J., Jung, T.P., Makeig, S., Brown, G., Kindermann, S.S., Lee, T.W., Sejnowski, T.J., 1998b. Spatially independent activity patterns in functional MRI data during the Stroop color-naming task. Proc. Natl. Acad. Sci. 95, 803–810. Mitra, P.P., Pesaran, B., 1999. Analysis of dynamic brain imaging data. Biophys. J. 76, 691–708. Mountcastle, V.B., 1957. Modality and topographic properties of single neurons of cats somatic sensory cortex. J. Neurophysiol. 20, 408–434. Murakami, H., Kamatani, D., Hishida, R., Takao, T., Kudoh, M., Kawaguchi, T., Tanaka, R., Shibuki, K., 2004. Short-term plasticity visualized with flavoprotein autofluorescence in the somatosensory cortex of anaesthetized rats. Eur. J. Neurosci. 19, 1352–1360. Orbach, H.S., Cohen, L.B., Grinvald, A., 1985. Optical mapping of electrical activity in rat somatosensory and visual cortex. J. Neurosci. 5, 1886–1895. Press, W., Teukolsky, S., Vetterling, W., Flannery, B., 1992. Numerical Recipes in C. Cambridge University Press, Cambridge, New York. Rector, D.M., Rogers, R.F., George, J.S., 1999. A focusing image probe for assessing neural activity in vivo. J. Neurosci. Methods 91, 135–145. Rector, D.M., Rogers, R.F., Schwaber, J.S., Harper, R.M., George, J.S., 2001. Scattered-light imaging in vivo tracks fast and slow processes of neurophysiological activation. NeuroImage 14, 977–994. Ressler, K.J., Sullivan, S.L., Buck, L.B., 1994. Information coding in the olfactory system—Evidence for a stereotyped and highly organized epitope map in the olfactory-bulb. Cell 79, 1245–1255. Roberts, S., Everson, R., 2001. Principal Component Analysis—Principles and Practice. Cambridge University Press, Cambridge.

Schiessl, I., Stetter, M., Mayhew, J.E., McLoughlin, N., Lund, J.S., Obermayer, K., 2000. Blind signal separation from optical imaging recordings with extended spatial decorrelation. IEEE Trans. Biomed. Eng. 47, 573–577. Shoham, D., Glaser, D.E., Arieli, A., Kenet, T., Wijnbergen, C., Toledo, Y., Hildesheim, R., Grinvald, A., 1999. Imaging cortical dynamics at high spatial and temporal resolution with novel blue voltage-sensitive dyes. Neuron 24, 791–802. Slovin, H., Arieli, A., Hildesheim, R., Grinvald, A., 2002. Long-term voltage-sensitive dye imaging reveals cortical dynamics in behaving monkeys. J. Neurophysiol. 88, 3421–3438. Spors, H., Grinvald, A., 2002. Spatio-temporal dynamics of odor representations in the Mammalian olfactory bulb. Neuron 34, 301–315. Spors, H., Wachowiak, M., Cohen, L.B., Friedrich, R.W., 2006. Temporal dynamics and latency patterns of receptor neuron input to the olfactory bulb. J. Neurosci. 26, 1247–1259. Stetter, M., Schiessl, I., Otto, T., Sengpiel, F., Hubener, M., Bonhoeffer, T., Obermayer, K., 2000. Principal component analysis and blind separation of sources for optical imaging of intrinsic signals. NeuroImage 11, 482–490. Stetter, M., Greve, H., Galizia, C.G., Obermayer, K., 2001. Analysis of calcium imaging signals from the honeybee brain by nonlinear models. NeuroImage 13, 119–128. Vanzetta, I., Slovin, H., Omer, D.B., Grinvald, A., 2004. Columnar resolution of blood volume and oximetry functional maps in the behaving monkey: implications for fMRI. Neuron 42, 843–854. Vassar, R., Chao, S., Sitcher, Nunez, J.M., Vosshall, L.B., Axel, R., 1994. Topographic organization of sensory projections to the olfactory-bulb. Cell 79, 981–991. Wachowiak, M., Cohen, L.B., 2001. Representation of odorants by receptor neuron input to the mouse olfactory bulb. Neuron 32, 723–735. Wachowiak, M., Denk, W., Friedrich, R.W., 2004. Functional organization of sensory input to the olfactory bulb glomerulus analyzed by twophoton calcium imaging. Proc. Natl. Acad. Sci. 101, 9097–9102. Woolsey, T.A., Van der Loos, H., 1970. The structural organization of layer IV in the somatosensory region (S I) of mouse cerebral cortex: the description of a cortical field composed of discrete cytoarchitectonic units. Brain Res. 17, 205–242.