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.

^{2}, and intensities −1.0 and +1.0 reflect 0.5 and 40.0 cd/m

^{2}, 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).

^{2}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.

**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.

*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:

*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.).

*k*_{ i}, where

*i*is the index over spikes—see Figure 1 for a schematic depiction. The spike-triggered ensemble,

**K**= [

**k**

_{ 1}

**k**

_{ 2}…], 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.

*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,

_{1}

_{2}…] 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:

**S**= [

**s**

_{ 1}

**s**

_{ 2}

**s**

_{ 3}…] 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.

**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.

**C**) of the spike-triggered distribution

*m*is the number of spikes and

**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.

_{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.

^{ T}/

*m*=

**EDE**

^{ T}, where

**E**is a matrix with the eigenvectors,

**e**

_{ i}, as its columns,

**E**= [

**e**

_{1},

**e**

_{2},…],

**D**is a diagonal matrix of corresponding eigenvalues,

**D**= diag(

*d*

_{1},

*d*

_{2},…), and

*m*is the total number of spikes. The pre-processed spike-triggered ensemble,

**K**

_{Ψ}, is calculated as

**E**

_{ n}= [

**e**

_{1},…,

**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*

_{1},…,

*d*

_{ n}). The matrix

**D**

_{ n}

^{−1/2}is computed by an element-wise operation as

**D**

_{ n}

^{−1/2}= diag(

*d*

_{1}

^{−1/2},…,

*d*

_{ n}

^{−1/2}).

*n*is calculated using the significance criterion described in the previous section.

*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.

**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*_{1}

**′**

*l*_{2}…

**′**

*l*_{ n}]:

**L**

^{ T}

**L**=

**I**and

**Z**=

**L**

^{ T}

**K**

_{Ψ}. The elements of the

*j*th row of

**Z**,

**z**

_{ j}are the projections of the observations

**K**

_{Ψ}along the

*j*th linear filter

**′**

*l*_{ j},

**K**

_{Ψ}, our goal is to recover the mixing matrix

**L**.

**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}.

**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)

*w*

_{ j}, of each filter.

*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*

_{1}

*g*

_{2}…

*g*

_{ R}] and

**w**= [

*w*

_{1}

*w*

_{2}…

*w*

_{ R}].

*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.

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 |