Receptive field characterization by spike-triggered independent component analysis

Abstract

The spikes generated by a neuron in response to stimuli provide information about the nature of the stimuli and also about the functional organization of the circuit in which the neuron is embedded. Spike-triggered analysis techniques such as spike-triggered covariance (STC) have been proposed to characterize the receptive field properties of neurons. So far, they have been able to provide only limited information about the functional organization of neural circuitry; in particular, STC tends to generate subfields that are mixed observations of independent processes. We address this problem by adding a criterion that sources are independent, resulting in an approach we call spike-triggered independent component analysis (ST-ICA). The method exploits the central limit theorem to find the directions in the high-dimensional stimulus space of spike-triggered data that are most independent. We demonstrate the improvement of the ST-ICA method over STC analysis using simulated neurons. When tested on data obtained from the H1 neuron in the fly visual system, it predicts a spatial arrangement of functional subunits with adjacent receptive fields. The properties of these subunits strongly resemble the known properties of elementary movement detector inputs to the H1 neuron. Using the ST-ICA method, we derive a model that captures functional and physiological properties of fly motion vision.

Introduction

One of the goals of computational neuroscience is to understand how sensory information is represented and processed in the nervous system. A core strategy has been to develop quantitative stimulus–response models of sensory systems. To gain substantial insight from this approach, we need to be able to relate the structure of the resulting models to the functional or physiological properties of the brain. 

The classical approach to the quantification of sensory responses is to count the number of spikes fired by a neuron as individual stimulus dimensions are varied sequentially (Hubel & Wiesel, 1962, 1998). A typical stimulus of this kind in the visual domain is a drifting sine-wave grating. The direction of motion that maximizes the neuron's firing rate is found by systematically changing the orientation. The next step is often to study the neuron's spatial tuning by changing the spatial frequency of gratings moving along the preferred direction and so on. This approach assumes that the stimulus dimensions are separable, and it also assumes substantial prior knowledge of the relevant stimuli with which to probe the cell. An alternative approach is spike-triggered analysis, which relaxes some of these assumptions, while in some cases imposing others. This approach has its origins in the Wiener–Volterra theory of nonlinear systems identification (Bialek & de Ruyter van Steveninck, 2005; Brenner, Bialek, & de Ruyter van Steveninck, 2000; de Ruyter van Steveninck & Bialek, 1988; Marmarelis & Marmarelis, 1978; Schwartz, Pillow, Rust, & Simoncelli, 2006; Theunissen, Sen, & Doupe, 2000; Wu, David, & Gallant, 2006). The simplest type of spike-triggered analysis is the spike-triggered average (STA), in which the first moment of the spike-triggered stimulus distribution indicates the average feature in stimulus space that causes the neuron to fire (Chichilnisky, 2001; Citron, Emerson, & Levick, 1988). When the system is linear, or approximately linear followed by a static nonlinearity, then this is the most common stimulus feature eliciting a spike. In this case, the linear filter obtained is the “receptive field” and together with the static nonlinearity, completely defines the system. The STA approach can be readily applied to simple cells in the mammalian primary visual cortex (Ringach, 2004). However, even in this cortical area, phase averaging (complex) cells are found, whose characterization requires the estimation of higher order terms. The H1 tangential cell in the fly lobula plate is also a phase averaging cell and thus cannot be characterized fully by the STA method which does not take higher order terms into account. 

The second order Wiener–Volterra term, or covariance matrix, can be computed relatively easily and we can calculate the dimensions that modulate the neuronal firing. The spike-triggered covariance (STC) approach applies eigenvalue decomposition to the covariance matrix and determines the subspace of orthogonal stimulus dimensions along which the variance of the neural response is maximal or minimal (Bialek & de Ruyter van Steveninck, 2005; Brenner et al., 2000; Rust, Schwartz, Movshon, & Simoncelli, 2004; Schwartz et al., 2006; Simoncelli, Paninski, Pillow, & Schwartz, 2004). Studies using STC methods have typically recovered a series of parallel filters for a neuron (Bialek & de Ruyter van Steveninck, 2005; Rust, Schwartz, Movshon, & Simoncelli, 2005; Touryan, Felson, & Dan, 2005). Although the models derived using STC give good predictions of the stimulus–response relationship, their relationship to functional circuit properties is unclear. 

Our goal is to relate the spike-triggered stimulus distribution to the structure of the functional circuit, or subunits carrying out components of processing in the system. STC recovers stimulus dimensions along which the variance is maximal. These dimensions cannot be attributed to individual functional units but reflect the combined contribution of all units. Independent component analysis (ICA) is a method that has been used to identify sources from observations that are mixtures of the sources (Simoncelli & Olhausen, 2001). Using ICA in combination with a white noise stimulus, we can find functional structure of the input elements to a neuron that can be described by orthogonal linear filters and symmetric nonlinearities. 

We describe the ICA algorithm as it applies to the spike-triggered framework and use it to analyze data recorded from the fly H1 neuron (e.g., Hausen, 1984). H1 is a heterolateral interneuron in the motion-processing system that responds predominantly to horizontal back-to-front motion within an entire visual hemisphere (e.g., Krapp, Hengstenberg, & Egelhaaf, 2001). It receives retinotopically organized input from small field elements which analyze the direction of local motion and whose properties are thought to be equivalent to those of elementary movement detectors (EMDs) (Reichardt, 1957; rev: Borst & Egelhaaf, 1993; Buchner, 1984; Reichardt, 1987). Here, we show that the filters derived from the spike-triggered independent component analysis (ST-ICA) are consistent with the functional substructure of EMDs. 

Experimental methods

Stimuli

We used a stimulus comprised of a number (B) of parallel adjacent bars, each of which had in any given frame an intensity value of either −1.0 (black) or 1.0 (white). The overall average luminance indicated by intensity 0.0 was maintained at 20 cd/m ², and intensities −1.0 and +1.0 reflect 0.5 and 40.0 cd/m ², respectively. The stimulus was generated using Expo (developed by Prof Peter Lennie, University of Rochester) running on a dual processor Apple Macintosh computer and projected onto a rear projection screen measuring 100 × 70 cm with a DepthQ 3D Projector (InFocus, Oregon, USA). 

The ideal unbiased stimulus for spike-triggered analysis is a spherically symmetric, Gaussian white noise stimulus (Paninski, 2003; Simoncelli et al., 2004). We found that such a stimulus elicits only a weak response in the H1 neuron and requires exceedingly long experiments to effectively recover all the filters. Instead, we used the cubic symmetric stimulus described above that elicits a stronger response. This stimulus is adequate for the spike-triggered analysis, as it has zero average and equal variance across all stimulus dimensions (Rust et al., 2004). To further verify, this we used the same stimulus on simulated linear–nonlinear–linear Poisson (LNLP) models and the filters recovered are not biased in any direction. Stimuli with unequal variance along the dimensions have also been used (Touryan et al., 2005) but are not guaranteed to give unbiased estimates. 

The stimulus presented to the fly was comprised of B = 10 bars, each covering 1.29° by 30° of the visual field. The row of bars was oriented perpendicular to the local preferred direction of the H1 neuron (Krapp et al., 2001). The stimulus was centered at 0° elevation and 45° azimuth and covered an angle of 13° along the horizontal direction and 30° along the vertical. The horizontal extent of the stimulus is shown with respect to the fly's hexagonal eye lattice structure in Figure 1a. The contrast of the frames was updated at 120 Hz. At the low average luminance conditions of 20 cd/m² there was no phase locking of the spikes to the stimulus presentation. For the model simulations we used a finer spatial resolution (B = 20 bars), as we are not limited by the time of presentation. 

Figure 1

 

Stimulation. (a) An example stimulus frame. In the fly experiments, each bar had a width of 1.29° and height of 30°. The frame is displayed over the eye facet (ommatidial) structure of the compound eye. (b) A segment of the stimulus and cartoon response shown across time. The regions enclosed by gray dots indicate stimulus histories, _i, leading to a spike. (c) A collection of these histories forms the spike-triggered ensemble. We apply STC and ST-ICA analysis methods to this ensemble. (d) Illustration of the calculation of STA and covariance axes (C₁ and C₂) relating to a set of spike-generating data points, K, indicated in red (Note: In this illustration, the STA is not projected before finding the covariance axes out, as described in Equation 4).

We refer to the stimulus frames preceding any point in time as the “stimulus history” (see Figure 1b). With B bars and a history of T frames, the stimuli ( s _i, where i refers to the index of the stimulus) are points in a B × T dimensional space. We ensure that the stimulus ensemble is unbiased in all of these dimensions; that is, the average across all stimuli is zero, and the variance is the same along every dimension. 

Linear–nonlinear–linear Poisson (LNLP) model

To test the applicability of the ST-ICA method we simulated the response of a motion sensitive neuron using a linear–nonlinear–linear Poisson (LNLP) model. The LNLP model consists of a set of orthogonal linear filters operating on the visual input, symmetric nonlinear functions, a linear summation, and a Poisson spike generator as illustrated in Figure 2. The input passes through the parallel series of orthogonal linear filters ( l _j) and cascades down to the symmetric nonlinear function ( f _j(.)). The outputs of the nonlinear functions are weighted ( w _j) and linearly summed before entering the Poisson generator. The Poisson generator mimics the stochastic nature of a neuron's spike generation mechanism. The probability for an LNLP model to spike, given a certain stimulus, can be written as:  

$P(\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{k}\mathrm{e}|\mathbf{s})=\sum _j{w}_j{f}_j\left({\mathit{l}}_j^T\mathbf{s}\right).$

(1)

 

Figure 2

 

Illustration of the linear–nonlinear–linear Poisson (LNLP) model. The input first goes through a set of parallel linear filters ( l _j) and through nonlinear functions f _j(.) to calculate the outputs. The final spike generator signal is a weighted sum of these outputs. A Poisson process is used to simulate the stochastic firing nature of the neuron.

The LNLP class of models with orthogonal linear filters and symmetric nonlinear functions have been used previously to model the responses of neurons. The most popular among these is the energy model of complex cells originally described by Adelson and Bergen (1985). Some cascade models of neurons in primate visual areas V1 and MT also fall under this class of models (Heeger, Simoncelli, & Movshon, 1996; Rust, Mante, Simoncelli, & Movshon, 2006). We used a simulation of the LNLP model with orthogonal linear filters and symmetric nonlinear functions to test the ST-ICA method for its computational feasibility and effectiveness in recovering the components of the model. The functional subunits that we are trying to find using the ST-ICA method are the linear–nonlinear (LN) stages of the LNLP model. The response of the fly H1 neuron is modeled well by a Reichardt correlation detector, which, depending on the filters used, can be equivalent to the energy model of motion sensitive complex cells (Adelson & Bergen, 1985; Borst & Egelhaaf, 1993; Egelhaaf, Borst, & Reichardt, 1989). As the energy model falls under the LNLP class of models (Adelson & Bergen, 1985), we fit the data from electrophysiology recordings of the fly H1 neuron to the LNLP class of models. 

Fly preparation and electrophysiology

The ST-ICA method was applied to the H1 neuron of the fly visual system. Experiments were carried out in 1- to 2-week-old female blowflies, Calliphora vicina, at a temperature between 18 and 22°C. The flies were immobilized in wax and the head was aligned with the horizontal plane of the projection screen using the symmetrical deep pseudopupil in the frontal region of both eyes (Franceschini, 1975). A small opening was made in the left head capsule and fat tissue and air sacs were cleared to gain access to the lobula plate. The nervous tissue was kept moist using a Ringer solution (Karmeier, Tabor, Egelhaaf, & Krapp, 2001). We recorded extracellular activity from the H1 neurons of four flies for a period of 30 to 90 minutes, with an average stimulus time of 45 min. The recordings were made using 2-MΩ tungsten electrodes (FHC Inc.) amplified by an A-M Systems amplifier and sampled at 25 kHz. Spike sorting was carried out online using a template-matching algorithm on Spike2 (CED Ltd.). 

Spike-triggered methods

We denote the spike-triggered stimulus histories as k _i, where i is the index over spikes—see Figure 1 for a schematic depiction. The spike-triggered ensemble, K = [ k ₁ k ₂…], consists of these spike-triggered stimulus histories as its column vectors. Spike-triggered methods analyze the statistical distribution of this ensemble to understand the relationship between stimuli and response. We present the STA and STC methods as background before introducing the ST-ICA method as a potential tool for recovering putative functional units in the underlying neural circuit. 

Spike-triggered average (STA)

When a white noise stimulus is used, the average of all points in the spike-triggered ensemble is its STA:  

$\mathbf{a}=\frac{1}{m}\sum _{i=1}^m{\mathbf{k}}_i,$

(2)

where m is the number of spikes. The STA is illustrated in Figure 1d. For nonwhite stimuli, the points in the spike-triggered ensemble need to be whitened with respect to the original stimulus covariance before starting the analysis. This is to prevent the recovered properties of the neuron being biased towards dominant features in the stimulus distribution. To do so, we first carry out the eigenvalue decomposition of the covariance matrix of the stimulus,  

$\stackrel{⌣}{\mathbf{S}}{\stackrel{⌣}{\mathbf{S}}}^T={\mathbf{E}}_{\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{m}}{\mathbf{D}}_{\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{m}}{\mathbf{E}}_{\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{m}}^T,$

(3)

where

$\stackrel{⌣}{\mathbf{S}}$

= [

$\stackrel{⌣}{\mathbf{s}}$

₁

$\stackrel{⌣}{\mathbf{s}}$

₂…] is a matrix of all stimuli, E _stim is the matrix whose columns are the eigenvectors, and D _stim is the diagonal matrix of the corresponding eigenvalues, D _stim = diag( d _stim1, d _stim2,…). The whitened stimulus is:  

$\mathbf{S}={\mathbf{E}}_{\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{m}}{\mathbf{D}}_{\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{m}}^{-1/2}{\mathbf{E}}_{\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{m}}^T\stackrel{⌣}{\mathbf{S}},$

(4)

where S = [ s ₁ s ₂ s ₃…] and s _i are the individual stimulus histories. The matrix D _stim ^−1/2 is computed by an element-wise operation, D _stim ^−1/2 = diag( d _stim1 ^−1/2, d _stim2 ^−1/2,…). The stimuli in whitened space that elicit a spike are columns vectors of the spike-triggered ensemble. Theunissen et al. (2000) and Touryan et al. (2005) give a more detailed description of this process. 

The STA analyses the first-order moment in the spike-triggered ensemble. Hence it is effective only when the receptive field of a neuron can be characterized by a single linear filter followed by an asymmetric nonlinearity. 

Spike-triggered covariance (STC)

As described above, the probability of a spike being generated in an LNLP model is given by Equation 1. A column in the spike-triggered ensemble, K, is a draw from the distribution P( s∣spike). Because the stimulus, s, we use is white (or whitened), it has an equal variance along every stimulus direction. When f _j(.)'s are symmetric excitatory functions, the distribution of K would have a high variance along the corresponding l _j's direction, and symmetric suppressive f _j(.)'s cause a low variance along the corresponding l _j's. By looking for the directions in the spike-triggered distribution with a high or low variance, we can find excitatory or suppressive filters of the neuron. 

To find the directions with high or low variance, we calculate the covariance matrix ( C) of the spike-triggered distribution  

$\mathbf{C}=\frac{\tilde{\mathbf{K}}{\tilde{\mathbf{K}}}^T}{m},$

(5)

where m is the number of spikes and  

$\tilde{\mathbf{K}}=[{\tilde{\mathbf{k}}}_1{\tilde{\mathbf{k}}}_2\dots {\tilde{\mathbf{k}}}_m],$

(6)

with the spike-triggered stimulus histories as the columns vectors of the matrix.  

${\tilde{\mathbf{k}}}_i={\mathbf{k}}_i-\widehat{\mathbf{a}}\left({\mathbf{k}}_i^T\widehat{\mathbf{a}}\right),$

(7)

 

$\mathrm{a}\mathrm{n}\mathrm{d}\widehat{\mathbf{a}}=\frac{\mathbf{a}}{\left|\mathbf{a}\right|}.$

(8)

 

$\widehat{\mathbf{a}}$

is the unit vector in the direction of the STA. We use

$\tilde{\mathbf{K}}$

so that the spike-triggered ensemble is normalized with respect to the STA, and the calculated STC axes are orthogonal to the STA (Schwartz et al., 2006). Applying eigenvalue decomposition to the covariance matrix (C), we calculate a set of orthogonal eigenvectors and their corresponding eigenvalues. The excitatory and suppressive filters can be chosen by a significance criterion described in the subsection below. 

The filters obtained by eigenvalue decomposition of this covariance matrix give us the subspace of stimulus directions along which the neuron is most or least responsive, which we will hereafter refer to as its receptive field (Schwartz et al., 2006). Using the filters as the linear stages (_j), we fit the neuronal response to an LNLP model. We use a histogram-based approach to recover the nonlinear function of the LNLP model as described in the Recovering nonlinearities and weights section. 

STC is effective when the stimulus is unbiased across the dimensions of the stimulus. This requires the stimulus to have a zero mean and be Gaussian of equal variance (Paninski, 2003). However, empirically biases due to non-Gaussian nature are often found to be small. Other stimuli, like natural scenes, have been used with the STC analysis by correcting for any stimulus-introduced covariance (Touryan et al., 2005). However, the correction does not guarantee unbiased estimates. STC recovers stimulus dimensions along which the variance is maximal or minimal. The recovered linear–nonlinear stages cannot be attributed to individual functional units but are the combined effect of all units contributing to the properties of the circuit. 

Significance testing

The eigenvectors are considered significant if the eigenvalues deviate substantially from the trend followed by other eigenvalues, shown in Figure 3a. To identify these points, we rank them (largest first) and find their linear trend. We then subtract the corresponding point along the line for each rank to get the corrected eigenvalues for each eigenvector ( Figure 3b). Eigenvectors with a corrected eigenvalue greater than twice the standard deviation are considered to be significantly different from the trend. 

Figure 3

 

Finding significant eigenvectors and recovering nonlinear functions. (a) The original eigenvalues of the covariance matrix, ranked with the greatest eigenvalue given the lowest rank. The red line is the best linear fit to these points. (b) The corrected eigenvalues, which are calculated by subtracting the trend from the data points. (c) The corrected eigenvalues plotted with units of standard deviation ( SD); red lines bound the region of significant eigenvalues (significance criterion: 2 SD). (d) The process of recovering the nonlinear function (red curve): the histogram of the output of the linear filter to all stimuli (lighter gray) is compared with the output from the spike-generating stimuli (darker gray).

Pre-processing for ST-ICA

We need to pre-process the data before carrying out the ST-ICA. Firstly, to ensure the data has a zero mean we project out the STA from each of the points in the spike-triggered ensemble as described in Equation 7. Secondly, we need to project the data onto the STC subspace, Ψ, and whiten it. For this, we carry out the eigenvalue decomposition (EVD) of the covariance matrix

$\tilde{\mathbf{K}}\tilde{\mathbf{K}}$

^T/ m = EDE ^T, where E is a matrix with the eigenvectors, e _i, as its columns, E = [ e ₁, e ₂,…], D is a diagonal matrix of corresponding eigenvalues, D = diag( d ₁, d ₂,…), and m is the total number of spikes. The pre-processed spike-triggered ensemble, K _Ψ, is calculated as  

${\mathbf{K}}_{\Psi }={\mathbf{D}}_n^{-1/2}{\mathbf{E}}_n^T\tilde{\mathbf{K}}.$

(9)

where E _n = [ e ₁,…, e _n] is the matrix whose columns are the n significant eigenvectors and D _n is the diagonal matrix of the corresponding n eigenvalues, D _n = diag( d ₁,…, d _n). The matrix D _n ^−1/2 is computed by an element-wise operation as D _n ^−1/2 = diag( d ₁ ^−1/2,…, d _n ^−1/2). n is calculated using the significance criterion described in the previous section. 

Spike-triggered independent component analysis (ST-ICA)

In this subsection we define the terminology we use before describing the motivation behind using ICA to characterize a neuron's input elements from its recorded response. This is followed by a brief description of ICA and its implementation. Using ST-ICA we attempt to find the LNLP model ( Figure 2) with orthogonal linear filters and symmetric nonlinear functions that best describes the functional properties of a neuron. The model has n independent original sources, where n is the number of eigenvectors we find to be significant. STC analysis of the spike-triggered ensemble, K, reveals the subspace, Ψ, along which the ensemble variance is highest. The linear filters of the n original sources will lie within Ψ and can be described as l′ _j, where l′ _j = E _n ^T l _j with E _n containing the significant eigenvectors along its columns. The directions of the linear filters correspond to the original source directions, as the linear filters of an LNLP model are its only directional components. 

Projecting the spike-triggered ensemble onto Ψ and whitening, we obtain the observations K _Ψ, as described by Equation 9 in the Pre-processing for ST-ICA section. These observations can be considered as mixtures of an original distribution, Z, with the mixing matrix L = [ l′ ₁ l′ ₂… l′ _n]:  

${\mathbf{K}}_{\Psi }=\mathbf{L}\mathbf{Z}.$

(10)

As the linear filters are orthogonal to each other, L ^T L = I and Z = L ^T K _Ψ. The elements of the jth row of Z, z _j are the projections of the observations K _Ψ along the jth linear filter l′ _j,  

${\mathbf{z}}_j=l_j^{\prime T}{\mathbf{K}}_{\Psi }.$

(11)

Given the observations K _Ψ, our goal is to recover the mixing matrix L. 

If the original sources of the neuron were just linear filters and we use a Gaussian stimulus, the observations K _Ψ projected onto Ψ would also have a Gaussian distribution. However, as the original sources have symmetric nonlinear functions after the linear filter, the observations along the original source directions, z _j, will deviate from being Gaussian distributions. ST-ICA exploits the non-Gaussian nature of the distribution of observations along the original source directions to find the mixing matrix ( L). It recovers the optimal L that maximizes the non-Gaussian nature of the distributions z _j. 

ICA has been widely used to determine the mixture matrix from observations (Bell & Sejnowski, 1995; Hyvarinen & Oja, 2000; van Hateren & van der Schaaf, 1998). It relies on the central limit theorem, which states that the sum of independent, identically distributed random variables tends towards a normal distribution (see, e.g., Papoulis & Pillai, 2002). In multivariate space, if certain directions have identical independent distributions of points, directions that are rotations of them would tend towards having a normal distribution of points, since they can be represented as linear mixtures of the original directions. If distributions along the original source directions are non-Gaussian, distributions along rotations of them, the other directions in Ψ, will be more Gaussian. This implies that the observations have the most non-Gaussian distribution along the original source directions (this is strictly true only for identical distributions; Hyvarinen & Oja, 2000). ICA can be used to recover for the directions in Ψ along which the observations have the most non-Gaussian distribution. 

We use the FastICA package (available for download at www.cis.hut.fi/projects/ica/fastica) to implement ICA. The package is based on the FastICA algorithm, a fixed-point iteration scheme to maximize non-Gaussian nature of the distributions (Hyvarinen, 1999). We chose neg-entropy as the measure of non-Gaussian nature and implemented the Gaussian contrast function in the FastICA package (Hyvarinen, 1997). Our approach recovers the same number of filters as the number of significant eigenvectors. 

The ideal stimulus for ST-ICA is a Gaussian white-noise stimulus because it ensures that the method does not attribute any stimulus features to the receptive field structure. However, it is sufficient for the stimulus to be Gaussian in Ψ, the STC subspace, as ICA is implemented within this subspace. Since our stimulus was a cubic symmetric and not Gaussian, we made sure that the distribution of the stimuli was Gaussian in Ψ before implementing the ST-ICA. 

The original source directions recovered by the ST-ICA method are the most independent directions in spike-triggered ensemble. Although the stimuli and sources are independent, our observations, the spike-triggered ensemble, are not generally independent. This is because the observations are the summed affect of all sources and characteristics of one unit can “explain away” some of the characteristics of the others (Wellman & Henrion, 1993). The performance of ICA thus relies on the degree of independence in the distribution introduced by the symmetric nonlinear functions present in the original sources. 

Recovering the nonlinearities and weights

The probability for the LNLP model firing a spike as described in Equation 1  

$P(\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{k}\mathrm{e}|\mathbf{s})=\sum _j{w}_j{f}_j\left({\mathit{l}}_j^T\mathbf{s}\right).$

(12)

In this equation we know the stimulus s, and using the STC or ST-ICA analysis we recover l _j. When we use a stimulus that is unbiased in every direction and the linear filters are orthogonal, we are able to obtain f _j(.) from (Schwartz et al., 2006) 

$f_j\left({\mathit{l}}_j^T\mathbf{s}\right)\propto P(\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{k}\mathrm{e}|{\mathit{l}}_j^T\mathbf{s})=\frac{P\left({\mathit{l}}_j^T\mathbf{s}\right|\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{k}\mathrm{e}\left)P\right(\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{k}\mathrm{e})}{P\left({\mathit{l}}_j^T\mathbf{s}\right)}.$

(13)

Therefore, we observe that 

$f_j\left({\mathit{l}}_j^T\mathbf{s}\right)\propto \frac{P\left({\mathit{l}}_j^T\mathbf{s}\right|\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{k}\mathrm{e})}{P\left({\mathit{l}}_j^T\mathbf{s}\right)}.$

(14)

 

One example of this process is shown in Figure 3d. 

We use least squares regression to recover the weights, w _j, of each filter.  

$\mathbf{w}=\mathrm{O}\mathrm{u}\mathrm{t}\mathrm{p}\mathrm{u}\mathrm{t}\times {\mathbf{g}}^T\times (\mathbf{g}\times {\mathbf{g}}^T{)}^{-1},$

(15)

where g _j = f _j( l _j ^T s), “Output” is the firing rate obtained by convolving a Gaussian kernel with the spike train, g = [ g ₁ g ₂… g _R] and w = [ w ₁ w ₂… w _R]. 

Results

Using simulations of an LNLP model, we demonstrate the strengths of the ST-ICA method in recovering the functional structure of a neuron. We then apply the analysis to recordings from the fly H1 neuron and compare the results with those from the STC method. 

Model simulation

To understand the advantages or disadvantages of the ST-ICA and STC methods, we tested how well they recover the parameters of a simulated LNLP model with orthogonal linear filters and symmetric nonlinear functions. We tested the methods on two versions of the LNLP model illustrated in Figure 4. The models have two excitatory filters with receptive fields in different regions of space. We would therefore expect to see spatially separated receptive fields in the structure of the EMDs that input to the H1 neuron as they get inputs from distinct regions of the visual field (Buchner, 1984). Having the receptive fields in different regions of space also ensures that the filters are independent. The shapes of the spike-triggered stimulus distributions, along the direction of the linear filters, mainly vary as a result of the symmetric nonlinear function while the variance is also affected by the weights. We tested the two methods with spike-triggered stimulus distributions of different variance and/or symmetric nonlinear functions. 

Figure 4

 

Illustration of the LNLP model used for the simulation to test differences between STC and ST-ICA. There are two sources with each of them a having a linear (_j) and a nonlinear f_j(.) function.

In the first version of the model, we used different weights and the same nonlinearity with two filters, resulting in the same shape but different variance and density of points along the two filters. In this case, the filters obtained with both the ST-ICA and STC methods, shown in Figure 5, were good estimates of the original filters. STC differentiates filters based on differences in variance along the two axes and in this case, ST-ICA used the difference in the density of points. They both perform equally well at recovering the parameters of the model. 

Figure 5

 

Comparison of derived from the ST-PCA and ST-ICA methods. The model simulation in this case has unequal variance and density of points along the two original filter directions. (a) The projection of the points in the stimulus ensemble onto the two filters ( l ₁ and l ₂) being the axis. There is an asymmetry between the two axes with l ₂ having a lesser density of points. The plot also shows the projections of the recovered filters from the STC (blue) and ST-ICA (red) plotted on the same axis. Both methods give good predictions of the original filters. (b) The original filters used in the model and (c) The filters recovered using the two methods. Note that the predictions made by both the STC and ST-ICA methods are sign-invariant.

In the second version of the model, we kept the variance along the two axes similar but varied the nonlinear function ( Figure 6). The difference in the distribution of points along the two axes is clearly visible in Figure 6a. As the variance of the distribution along the two axes is equal, the STC method was unable to recover the original filters, but recovered the subspace within which they exist. Therefore the resulting STC filters were mixtures of the original filters ( Figures 6a and 6b). As ST-ICA is able to exploit the differences in higher-order moments, in this case kurtosis, its estimated filters matched well the original filters. 

Figure 6

 

Comparison of models derived from the STC and ST-ICA methods. Model simulation in this case has with equal variance along the two axes but different nonlinear functions. (a) The projection of the points in the stimulus ensemble onto the two filters ( l ₁ and l ₂) being the axis. A structure can be observed in the distribution of the points, especially along l ₂. The plot also shows the projections of the recovered filters from the STC (blue) and ST-ICA (red) plotted on the same axis. The ST-ICA axes are much closer to the original axes as compared to the STC axes. (b) The original filters used in the model, the filters recovered using the STC and ST-ICA methods. (c) The corresponding nonlinear function of panels a and b. Each of the STC filters is a mixture of the two original filters, while the ST-ICA is a good estimate of the original filters. The nonlinear functions recovered by the ST-ICA are also closer to the original than those recovered by STC. Note that the predictions made by both the STC and ST-ICA methods are sign-invariant.

The estimation of the nonlinear function, Equation 14, is based on the prediction of the linear stage; errors in the calculation of the linear stage propagate into the errors in the estimation of the nonlinear function. This is demonstrated in Figure 6c, where the ST-ICA estimates of the nonlinear functions are more accurate than the STC estimates. 

Fly experiments

By applying the ST-ICA method to recordings from the fly H1 neuron ( Figure 7), we found four excitatory filters together with their corresponding nonlinear functions and relative weights. The independent space-time filters that we found are shown on the left side of Figure 7. Gray regions in these filters indicate that the neuron did not produce a significant response, while regions that show shifts towards white (+1, ON) or black (−1, OFF) correspond to features of its receptive field. The STA of the neuron contains mostly gray regions, as we would expect from a phase averaging neuron. For each ST-ICA filter, we observe a latency of approximately 50 ms. The neuron responds most strongly from 50 ms to 90 ms. Within this period, we notice that with increasing time, the ON region of the receptive field shifts towards the left. In the stimulated eye region, this corresponds to horizontal back-to-front motion, which is in agreement with previous studies of H1's response (e.g., Krapp et al., 2001). The nonlinear functions calculated for all the filters are approximately squaring functions. 

Figure 7

 

Recovered LNLP model of the fly H1 neuron. (a) The ST-ICA method finds four sources (functional subunits) for H1. The figure shows the linear filter, nonlinear function and the weight corresponding to each source. (b) The STA of the neuron. As H1 is a phase averaging neuron, the STA does not have any significant features.

Each of the independent filters we recovered spans a slightly different region of space. The extent of each of these filters, as measured by the full width half maximum, is 3.7 ± 0.7° (11 filters across 4 flies) in the spatial dimension and 21 ± 4 ms (15 filters across 4 flies) in the temporal dimension. The filters that we derive are what we would expect to be the receptive fields of the dominant elementary movement detectors (EMDs) at low-light adapted conditions (Schuling, Mastebrock, Bult, & Lenting, 1989). The model for an EMD is illustrated in Figure 8a (Reichardt, 1987). The relative weights calculated (the w_j listed in Figure 7) indicate that the filters around 35° azimuth modulate the response of H1 more strongly than the those at 50° azimuth. This result is in agreement with previous electrophysiological studies (Krapp et al., 2001). Based on this information, the parameters measured in our experiment can be represented as shown in Figure 8b, with each filter representing one correlation detector, with a corresponding weight w_j. The suppressive filters, which are shown in gray, do not always have enough strength to affect firing beyond the significance criteria we set for filter recovery. 

Figure 8

 

Model for the fly H1. (a) An illustration of the correlation model for movement of stimulus across two neighboring eye facets. The gray segment is an inhibitory unit while the black segment is an excitatory unit. (b) An elaborate model where we illustrate a series of these detectors, each subunit is weighted before it is summed. w_p^r is the weight where p indicates whether it is excitatory (ex) or suppressive (in) and r is their rank according to spatial position.

Comparison of ST-ICA with STC

The different filters obtained using STC and ST-ICA are shown in Figure 9. In this example (the same neuron as shown in Figure 7) we applied the same significance criteria for both the methods. The ST-ICA filters thus lie within the subspace defined by the STC filters and are linear transformations of the STC filters. The STC filters in this case are observations based on a mixture of contributions from the independent filters. The difference between the two filters is most observable in the third and fourth filters; the ST-ICA filters seem better defined than their corresponding STC filters. The differences in the linear filters are better visualized in the spatial frequency—temporal frequency plots, where the STC filters have four peaks while the ST-ICA filters have only two peaks. This could be due to the mixing of some components from the first and second filters in the estimation of the third and fourth filters by STC. 

Figure 9

 

The comparison of the STC and ST-ICA methods on the fly H1 data. The left half of the figure shows the spatiotemporal representation of the filters recovered by the STC and the ST-ICA processes, with the corresponding spatial frequency—temporal frequency domain representations shown on the right. While the first two filters do not show much difference between the two methods, the third and fourth filters from the ST-ICA are better defined than the corresponding filters recovered from the STC method. The difference is more distinct in the frequency domain where we see two peaks in the ST-ICA filters as opposed to four in the STC method.

Discussion

We have here introduced the ST-ICA method, exploiting independent component analysis to characterize the receptive fields and functional structure of the input elements to a neuron that can be described by orthogonal linear filters and symmetric nonlinearities. We applied this method to the fly H1 neuron, whose functional substructure and receptive fields are fairly well understood (Buchner, 1984; Krapp & Hengstenberg, 1997). When applied to the fly H1 neuron, the ST-ICA method recovered filters selective for movement in the back-to-front direction within an area of approximately 3.7° along the azimuth (Figure 7). These filters agree with the known response properties of EMDs that input to the H1 neuron (Schuling et al., 1989). It should be mentioned, however, that the filters characterize the neuron's local receptive field properties only for the region where the stimulus was presented. As the H1 neuron is a wide-field neuron that responds to visual motion within an entire visual hemisphere, the filters obtained are likely to vary slightly across the neuron's receptive field (Krapp et al., 2001). Using the filters recovered from the ST-ICA method as the linear stages of an LNLP model, we calculate the nonlinear function operating on their output and the relative weights of each filter. Hence, we can use the method to build an LNLP model of the neuron being studied. 

To test the efficacy of recovering the structure from spike-triggered data, we simulated data using an LNLP model and compared the results obtained from the ST-ICA and STC methods to the original model. The filters recovered using the ST-ICA method match the original filters that are used in the simulation, while filters recovered by the STC method are mixed observations of the original filters. The occurrence of mixed observations is more frequent in cases where the variance of the response is similar in the directions of the different filters of the original model ( Figure 6). As they are mixed observations of the original filters, the STC filters are in the same subspace as the original filters. STC analysis can accurately predict the subspace of a neuron's receptive field in a high dimensional stimulus space. It has been successfully used to characterize the receptive field properties of V1 cortical cells using random stimuli (Rust et al., 2005) and natural scenes (Touryan et al., 2005). STC analysis uses variance along the different directions to select filters; while variance is effective at finding the subspace of relevant stimulus dimensions, it cannot be used as a criterion to differentiate functional units, as we have demonstrated in the example in Figure 6. The ST-ICA method instead uses independence as a criterion to make such a classification. The results we obtain show that the method effectively reconstructs simulated LNLP neural models with orthogonal linear filters and symmetric nonlinear functions (Figures 5 and 6). The ST-ICA method also recovers the filters underlying the functional input structure of the fly H1 neuron (Figure 9), demonstrating the validity of using the independence criterion. 

The response latency observed in the ST-ICA filters recovered for the H1 neuron is about 50 ms. Studies with high-contrast grating stimuli at high temperatures report response latencies around 25 ms (Krapp & Hengstenberg, 1997; Warzecha & Egelhaaf, 2000). Low temperature has been shown to increase the response latency to >50 ms (Warzecha & Egelhaaf, 2000; Warzecha, Horstmann, & Egelhaaf, 1999). Our range of temperatures (18–22°C) and the random noise stimulus used are the most likely explanation for the higher response latency observed. 

The STA of the H1 neuron has no significant features relative to its STC or ST-ICA filters ( Figure 7). This shows that H1 is a phase averaging neuron similar to the complex cells in the mammalian visual cortex (Hubel & Wiesel, 1962). However, there could be situations where a neuron may have both significant STA and significant ST-ICA filters. Such neurons were found in the mammalian visual cortex (Rust et al., 2005). Functional units that have the linear stage followed by an asymmetric nonlinearity would be included in the STA. As we normalize the spike-triggered stimulus ensemble with respect to the STA, before starting the ST-ICA method, we do not recover any of these features in the ST-ICA filters. To reconstruct a complete model of a neuron, which has a significant STA, we need to consider the STA as an additional independent source to the neuron we are characterizing. 

The ST-ICA method looks for the directions in stimulus space along which the distributions of points are most independent. In the case of excitatory filters, we look for the presence of independent distributions of points in different directions of stimulus space. For suppressive filters, the relevant directions are those along which there is an absence of points. Estimates of the shape of distributions along the suppressive filter directions would be greatly under-sampled as they are sparse distributions. Hence, the ST-ICA method will not be able to recover independent suppressive units of the circuit, and we will have to rely on the estimates of the subspace that can be calculated using the STC analysis. 

When looking for multiple sources, the FastICA algorithm repeats the iterative search starting at different points. To avoid converging to the same direction each time, it projects the previously found directions out of the data before starting the next iteration. This constrains the recovered filters to be orthogonal, limiting the ST-ICA analysis. There exist some techniques that attempt nonorthogonal ICA (Lee, Vrins, & Verleysen, 2006; Yeredor, 2002) that could make the analysis less constrained and may reduce the bias in estimating the functional units. 

The ST-ICA method can recover an LNLP model for a neuron when linear–nonlinear functions are good approximations of the corresponding functional units. It requires that the functional units operate in independent regions of stimulus space and have symmetric nonlinearities. ST-ICA strictly works only for identical distributions. The efficacy of the recovered model can be assessed from the accuracy of predicting the response of the neuron to test stimuli. But how do we verify the existence of recovered functional units? Functional units are subunits carrying out components of processing in the system. The biophysical implementation of this processing could be either within the dendritic arborization in pre-synaptic neurons or in much larger networks of neurons that form the neuron's input. To verify whether we have indeed recovered real functional units, we would have to find and simultaneously record from the individual units, which may be difficult to achieve using current experimental techniques. The ST-ICA method gives us an idea of what the functional structure of such input units looks like, in case we are using established experimental techniques like single unit recordings to study them. 

The spatial extent of the EMDs is reported to cover about the extent of 3.7° (Schuling et al., 1989) under low-light adapted conditions. We also know the directional preference of the EMDs that feed onto the tangential cells from the local preferred directions (Krapp et al., 2001). Using just an assumption of independent sources, we were able to recover a model that matches the previously reported functional substructure of the H1 neuron. In the last decade there have been many studies that apply system identification methods to analyze neural systems (reviewed in Wu et al., 2006) with the goal of estimating a function that relates neural response to sensory stimuli. These methods are limited with regard to understanding how the computation of this stimulus–response function is achieved by the functional structure of the system. The ST-ICA method that we propose here could provide additional details on the functional organization of neural circuits that map sensory stimuli onto cellular responses. Even though results obtained with the ST-ICA method may be difficult to interpret if applied to less well-understood sensory stems, it will help us to reveal their functional structure. 

Appendix A

Table: Variables used in the paper

Variable	Dimensions	Description
B	B	No. of bars in the stimulus
T	T	No. of stimulus time frames considered
s i	N = B × T	Stimulus history
l j	N	jth linear filter
f j(.)	scalar	jth nonlinear function
w j	scalar	jth weight
k i	N × 1	jth spike-triggered stimulus history
K	N × m	Spike-triggered ensemble
m	scalar	No. of spikes
a	N × 1	STA
S ⌣	N × no. of stimuli	Original stimuli
S	N × no. of stimuli	Stimuli corrected for stimulus covariance
E stim	N × N	Stimulus covariance eigenvectors
D stim	N × N	Stimulus covariance eigenvalues
d stim	scalar	Stimulus covariance eigenvalue
C	N × N	Spike-triggered covariance matrix
K ˜	N × m	Spike-triggered ensemble—STA
k ˜ i	N × 1	Spike-triggered observation—STA
a ^	N × 1	STA unit vector
E	N × N	STC eigenvector matrix
D	N × N	STC eigenvalue diagonal matrix
e i	N × 1	STC eigenvector
d i	scalar	STC eigenvalue
n	scalar	No. of significant eigenvalues
K Ψ	n × m	Pre-possessed spike-triggered ensemble
E n	N × n	n significant eigenvectors
D n	n × n	n significant eigenvalues diagonal matrix
Ψ	n-dim. space	STC subspace
l′ j	n × 1	The jth linear filter in Ψ
L	n × n	Linear filters in Ψ
Z	n × m	Original distribution
z j	1 × m	Distribution of spike-triggered ensemble along the jth filter
w	1 × n	Weight vector
g	1 × n	Outputs of the n linear–nonlinear stages
g j	scalar	Output of the jth linear–nonlinear stage

Acknowledgments

This research was funded by a grant from the Gatsby Charitable Trust (GAT2830) to SRS. We would like to thank Rasmus Peterson for critically reviewing the manuscript and Kit D. Longden, Julija Krupic, Naveed Ejaz, and Philip Bream for useful discussions and help with writing the manuscript. 

Commercial relationships: none. 

Corresponding author: Aman B. Saleem. 

Email: [e-mail]. 

Address: Department of Bioengineering, Imperial College London, South Kensington, London, SW7 2AZ, UK.