PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of ploscompComputational BiologyView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS Comput Biol. 2010 April; 6(4): e1000757.
Published online 2010 April 22. doi:  10.1371/journal.pcbi.1000757
PMCID: PMC2858697

Independent Component Analysis in Spiking Neurons

Abigail Morrison, Editor

Abstract

Although models based on independent component analysis (ICA) have been successful in explaining various properties of sensory coding in the cortex, it remains unclear how networks of spiking neurons using realistic plasticity rules can realize such computation. Here, we propose a biologically plausible mechanism for ICA-like learning with spiking neurons. Our model combines spike-timing dependent plasticity and synaptic scaling with an intrinsic plasticity rule that regulates neuronal excitability to maximize information transmission. We show that a stochastically spiking neuron learns one independent component for inputs encoded either as rates or using spike-spike correlations. Furthermore, different independent components can be recovered, when the activity of different neurons is decorrelated by adaptive lateral inhibition.

Author Summary

How the brain learns to encode and represent sensory information has been a longstanding question in neuroscience. Computational theories predict that sensory neurons should reduce redundancies between their responses to a given stimulus set in order to maximize the amount of information they can encode. Specifically, a powerful set of learning algorithms called Independent Component Analysis (ICA) and related models, such as sparse coding, have emerged as a standard for learning efficient codes for sensory information. These algorithms have been able to successfully explain several aspects of sensory representations in the brain, such as the shape of receptive fields of neurons in primary visual cortex. Unfortunately, it remains unclear how networks of spiking neurons can implement this function and, even more difficult, how they can learn to do so using known forms of neuronal plasticity. This paper solves this problem by presenting a model of a network of spiking neurons that performs ICA-like learning in a biologically plausible fashion, by combining three different forms of neuronal plasticity. We demonstrate the model's effectiveness on several standard sensory learning problems. Our results highlight the importance of studying the interaction of different forms of neuronal plasticity for understanding learning processes in the brain.

Introduction

Independent component analysis is a well-known signal processing technique for extracting statistically independent components from high-dimensional data. For the brain, ICA-like processing could play an essential role in building efficient representations of sensory data [1][4]. However, although many algorithms have been proposed for solving the ICA problem [5], only few consider spiking neurons. Moreover, the existing spike-based models [6], [7] do not answer the question how this type of learning can be realized in networks of spiking neurons using local, biologically plausible plasticity mechanisms (but see [8]).

Classic ICA algorithms often exploit the non-Gaussianity principle, which allows the ICA model to be estimated by maximizing some non-Gaussianity measure, such as kurtosis or negentropy [5]. A related representational principle is sparse coding, which has been used to explain various properties of V1 receptive fields [9]. Sparse coding states that only a small number of neurons are activated at the same time, or alternatively, that each individual unit is activated only rarely [10]. In the context of neural circuits, it offers a different interpretation of the goal of the ICA transform, from the perspective of metabolic efficiency. As spikes are energetically expensive, neurons have to operate under tight metabolic constraints [11], which affect the way information is encoded. Moreover, experimental evidence supports the idea that the activity of neurons in V1 is sparse. Close to exponential distributions of firing rates have been reported in various visual areas in response to natural scenes [12].

Interestingly, certain homeostatic mechanisms are thought to regulate the distribution of firing rates of a neuron [13]. These intrinsic plasticity (IP) mechanisms adjust ionic channel properties, inducing persistent changes in neuronal excitability [14]. They have been reported for a variety of systems, in brain slices and neuronal cultures [14], [15] and they are generally thought to play a role in maintaining system homeostasis. Moreover, IP has been found to occur in behaving animals, in response to learning (see [14] for review).

From a computational perspective, it is believed that IP may maximize information transmission of a neuron, under certain metabolic constraints [13]. Additionally, we have previously shown for a rate neuron model that, when interacting with Hebbian synaptic plasticity, IP allows the discovery of heavy-tailed directions in the input [16]. Here, we extend these results for a network of spiking neurons. Specifically, we combine spike-timing dependent plasticity (STDP) [17][19], synaptic scaling [20] and an IP rule similar to [16], which tries to make the distribution of instantaneous neuronal firing rates close to exponential.

We show that IP and synaptic scaling complement STDP learning, allowing single spiking neurons to learn useful representations of their inputs for several ICA problems. First, we show that output sparsification by IP together with synaptic learning is sufficient for demixing two zero mean supergaussian sources, a classic formulation of ICA. When using biologically plausible inputs and STDP, complex tasks, such as Foldiák's bars problem [21], and learning oriented receptive fields for natural visual stimuli, can be tackled. Moreover, a population of neurons learns to extract several independent components if the activity of different neurons are decorrelated by adaptive lateral inhibition. When investigating the mechanisms how learning occurs in our model, we show that IP is necessary for learning, as it enforces a sparse output, guiding learning towards heavy-tailed directions in the input. Lastly, for specific STDP implementations, we show that IP shifts the threshold between potentiation and depression, similar to a sliding threshold for Bienenstock-Cooper-Munro (BCM) learning [22].

The underlying assumption behind our approach, implicit in all standard models of V1 receptive field development, is that both input and output information are encoded in rates. In this light, one may think of our current work as a translation of the model in [16] to a spike-based version. However, the principles behind our model are more general than suggested by our work with rate neurons. We show that the same rule can be applied when inputs are encoded as spike-spike correlation patterns, where a rate-based model would fail.

Results

A schematic view of the learning rules is shown in Fig. 1. The stochastically spiking neuron [23] generates spikes as an inhomogeneous Poisson process, with mean expressed as a function of the total incoming current to the neuron An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e001.jpg, parametrized by variables An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e002.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e003.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e004.jpg. This transfer function is optimized by adapting the three parameters to make the distribution of instantaneous firing rates of the neuron approximatively exponential (for a complete mathematical formulation, see the Methods section). Additionally, Hebbian synaptic plasticity, implemented by nearest-neighbor STDP [24], changes incoming weights and a synaptic scaling mechanism keeps the sum of all incoming weights constant over time.

Figure 1
Overview of plasticity rules used for ICA-like learning.

A simple demixing problem

To illustrate the basic mechanism behind our approach, we first ask if enforcing a sparse prior by IP and Hebbian learning can yield a valid ICA implementation for the classic problem of demixing two supergaussian independent sources. In the standard form of this problem, zero mean, unit variance inputs ensure that the covariance matrix is the identity, such that simple Hebbian learning with a linear unit (equivalent to principal component analysis) would not be able to exploit the input statistics and would just perform a random walk in the input space. This is, however, a purely mathematical formulation, and does not make much sense in the context of biological neurons. Inputs to real neurons are bounded and —in a rate-based encoding— all positive. Nonetheless, we chose this standard formulation to illustrate the basic underlying principles behind our model. Below, we will consider different spike-based encodings of the input and learning with STDP.

As a special case of a demixing problem, we use two independent Laplacian distributed inputs, with unit variance: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e017.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e018.jpg. For the linear superposition, we use a rotation matrix An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e019.jpg:

equation image
(1)

where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e021.jpg is the angle of rotation, resulting in a set of inputs An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e022.jpg. Samples are drawn at each time step from the input distribution and are mapped into a total input to the neuron as An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e023.jpg, with the weight vector An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e024.jpg normalized). The neuron's transfer functions An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e025.jpg, the same as for our spiking model (Fig. 1), is adapted based on our IP rule, to make the distribution of firing rates exponential. For simplicity, here weights change by classic Hebbian learning: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e026.jpg, with An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e027.jpg being the synaptic learning rate (see Methods for details). Similar results can be obtained when synaptic changes follow the BCM rule.

In Fig. 2A we show the evolution of synaptic weights for different starting conditions. As our IP rule adapts the neuron parameters to make the output distribution sparse (Fig. 2B,C), the weight vector aligns itself along the direction of one of the sources. With this simple model, we are able to demix a linear combination of two independent sources for different mixing matrices and different weights constraints (Fig. 2D), as any other single-unit implementation of ICA.

Figure 2
A demixing problem: two rotated Laplace directions.

One neuron learns an independent component

After showing that combining IP and synaptic learning can solve a classical formulation of ICA, we focus on spike-based, biologically plausible inputs. In the following, STDP is used for implementing synaptic learning, while the IP and the synaptic scaling implementations remain the same.

Demixing with spikes

The demixing problem above can be solved also in a spike-based setting, after a few changes. First, the positive and negative inputs have to be separated into on- and off- channels (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e042.jpg) and converted into Poisson spike trains of a certain duration An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e043.jpg (see Methods), with the corresponding rates An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e044.jpg (note that the inputs are no longer white in this four-dimensional space). Secondly, to avoid the unbiological situation of having few very strong inputs, such that single presynaptic spikes always elicit a spike in the postsynaptic neuron, each channel consists of several (here, 25) synapses, with independent inputs having the same firing rate. Synapses are all positive and adapt by STDP, under the normalization An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e045.jpg. A more detailed description of the parameters can be found in the Methods section.

The evolution of the weights is shown in Fig. 3A. The corresponding receptive field of the neuron can be obtained by projecting the weight vector back onto the original two-dimensional space (details of this procedure are described in the Methods and Text S1). As in the original formulation of the problem, the neuron receptive field slowly aligns itself along one of the independent components (Fig. 3B).

Figure 3
The demixing problem with inputs encoded as spike trains.

Foldiák's bars: input encoded as firing rates

As a second test case for our model, we consider Foldiák's well-known bars problem [21]. This is a classic non-linear ICA problem and is interesting from a biological perspective, as it mimics relevant nonlinearities, e.g. occlusions. In the classic formulation, for a two-dimensional input An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e047.jpg, of size An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e048.jpg, a single bar must be learned after observing samples consisting of the nonlinear superposition of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e049.jpg possible individual bars (see Fig. 4A), each appearing independently, with probability An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e050.jpg (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e051.jpg). The superposition is non-linear, as the intersection of two bars has the same intensity as the other pixels in the bars (a binary OR).

Figure 4
Learning a single independent component for the bars problem.

In our implementation, the input vector is normalized and the value of each pixel An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e053.jpg is converted into a corresponding Poisson spike train of duration on the order of a fixation duration [25]. The details of the experimental setup are described in the Methods section. As the IP mechanism begins to take effect, making neuronal activity sparse, the receptive field of the neuron slowly adapts to one of the bars (Fig. 4B). This effect is robust for a wide range of parameters (see Text S2) and, as suggested by our previous results for rate neurons [16], does not critically depend on the particular implementation of the synaptic learning. We obtain similar results with an additive [26] and a simple triplet [27] model of STDP.

The input normalization makes a bar in an input sample containing a single IC stronger than a bar in a sample containing multiple ICs. This may suggest that a single component emerges by preferentially learning ‘easy’ examples, i.e. examples with a single bar. However, this is not the case and one bar is learned even when single bars never appear in the input, as in [28]. Specifically, we use a variant of the bars problem, in which the input always consists of 4 distinct bars, selected at random. In this case, the neuron correctly learns a single component. Moreover, similar results can be obtained for 2–5 bars (see Text S2), using the same set of parameters.

Foldiák's bars: input encoded by spike-spike correlations

In the previous experiments, input information was encoded as firing rates. In the following, we show that this stimulus encoding is not critical and the presence of spike-spike correlations is sufficient to learn an independent component, even in the absence of any presynaptic rate modulation. To demonstrate this, we consider a slight modification of the above bars problem, with input samples consisting of always two bars, each two pixel wide, with N distinct bars. This is a slightly more difficult task, as wide bars emphasize non-linearities due to overlap, but can be solved by our model with a rate encoding (see Text S2).

In this case, all inputs have the same mean firing rate and the information about whether a pixel belongs to a bar or not is encoded in correlation patterns between different inputs (see Methods). Specifically, background inputs are uncorrelated, while inputs belonging to bars are all pairwise correlated, with a correlation coefficient An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e054.jpg (see Fig. 5B).

Figure 5
Bars in a correlation-based encoding.

As for the rate-based encoding, the neuron is able to reliably learn a bar (Fig. 5C, D). Similar results were obtained for the version with always two bars, each one pixel wide, but with slower convergence. The fact that our approach also works for correlated inputs rests on the properties of STDP and IP. In the original rate formulation, strong inputs lead to higher firing in the postsynaptic neuron, causing the potentiation of the corresponding synapses. Similarly, if presynaptic inputs fire all with the same rate, correlated inputs are more successful in driving the neuron and hence their weights are preferentially strengthened [26]. Moreover, as before, IP enforces sparse post-synaptic responses, guiding learning towards a heavy-tailed direction in the input (a more formal analysis of the interaction between IP and STDP is presented below).

Due to its stochastic nature, our neuron model is not particularly sensitive to input correlations, hence An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e056.jpg cannot be too small. Stable receptive fields with a single 2 pixel wide bar are still obtained for lower correlation coefficients (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e057.jpg), but with slower convergence. We expect better results with a deterministic neuron model, such as the leaky integrate-and-fire. An approximation of IP, based on moment matching could be used in this case. Additionally, STDP-induced competition (due to the relatively small number of inputs, in some cases one weight grows big enough to elicit a spike in the postsynaptic neuron) [26] enforces some constraints on the model parameters to ensure the stability of a solution with multiple non-zero weights. This could be done by increasing the size of the input, restricting the overall mean firing of the neuron An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e058.jpg or by slightly changing the STDP parameters (see Methods). These parameter changes do not affect learning with the rate-based encoding, however.

Natural scenes input

The third classical ICA problem we consider is the development of V1-like receptive fields for natural scenes [3]. Several computational studies have emphasized that simple cell receptive fields in V1 may be learned from the statistics of natural images by ICA or other similar component extraction algorithms [9], [29], [30]. We hypothesized that the same type of computation could be achieved in our spiking neuron model, by combining different forms of plasticity. Only a rate encoding was used for this problem, partly for computational reasons and partly because it is not immediately obvious how a correlation-based encoding would look like in this case.

We use a set of images from the van Hateren database [31], with standard preprocessing (see Methods). The rectified values of the resulting image patches are linearly mapped into a firing frequency for an on- and off-input population, as done for the bars. The STDP, IP and other simulation parameters are the same as before.

As shown in Fig. 6A, the receptive field of the neuron computed as the difference between the weights of the on- and the off- input populations (depicted in Fig. 6B) evolves to an oriented filter, similar to those obtained by other ICA learning procedures [9], [29], [30], [32], [33]. A similar receptive field can be obtained by reverse correlation from white noise stimuli. The least-mean-square-error fit to a Gabor wavelet [34] is shown in Fig. 6C. As vertical edges are usually over-represented in the input, the neuron will typically learn a vertical edge filter, with a phase shift depending on the initial conditions. The receptive field has low spatial frequency, but more localized solutions result for a neural population (see below).

Figure 6
Learning a Gabor-like receptive field.

ICA in a neuron population

So far, learning has been restricted to a single neuron. For learning multiple independent components, we implement a neuron population in which the activities of different neurons are decorrelated by adaptive lateral inhibition. This approach is standardly used for feature extraction methods based on single-unit contrast functions [30]. Here, we consider a simple scheme for parallel (symmetrical) decorrelation. The all-to-all inhibitory weights (Fig. 7A) change by STDP and are subject to synaptic scaling, as done for the input synapses. We only use a rate-based encoding for this case, due to computational overhead, which also limits the size of networks we can simulate.

Figure 7
Learning multiple ICs.

We consider a population of 10 neurons. In order to have a full basis set for the bars problem, we use 2 pixel wide bars. For this case, our learning procedure is able to recover the original basis (Fig. 7B). As lateral inhibition begins to take effect, the average correlation coefficient between the responses of different neurons in the population decreases (Fig. 7C), making the final inhibitory weights unspecific (Fig. 7D). As decorrelation is not a sufficient condition for independence, we show that, simultaneously, the normalized mutual information decreases (see Methods for details). Using the same network for the image patches, we obtain oriented, localized receptive fields (Fig. 7E).

Due to the adaptive nature of IP, the balance between excitation and inhibition does not need to be tightly controlled, allowing for robustness to changes in parameters. However, the inhibition strength influences the time required for convergence (the stronger the inhibition, the longer it takes for the system to reach a stable state). A more important constraint is that the adaptation of inhibitory connections needs to be faster than that of feedforward connections to allow for efficient decorrelation (see Methods for parameters).

Is IP necessary for learning?

We wondered what the role of IP is in this learning procedure. Does IP simply find an optimal nonlinearity for the neuron's transfer function, given the input, something that could be computed offline (as for InfoMax [29]), or is the interaction between IP and STDP critical for learning? To answer this question, we go back to Foldiák's bars. We repeat our first bars experiment (Fig. 4) for a fixed gain function, given by the parameters obtained after learning (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e061.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e062.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e063.jpg). In this case, the receptive field does not evolve to an IC (Fig. 8). This suggests that ICA-like computation relies on the interplay between weight changes and the corresponding readjustment of neuronal excitability, which forces the output to be sparse. Note that this result holds for simulation times significantly larger than in the experiment before, where a bar emerged after An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e064.jpg, suggesting that, even if the neuron would eventually learn a bar, it would take significantly longer to do so.

Figure 8
IP is critical for learning.

We could assume that the neuron failed to learn a bar for the fixed transfer function just because the postsynaptic firing was too low, slowing down learning. Hence, it may be that a simpler rule, regulating just the mean firing rate of the neuron, would suffice to learn an IC. To test this hypothesis, we construct an alternative IP rule, which adjusts just An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e065.jpg to preserve the average firing rate of the neuron (see Methods). With the same setup as before and the new IP rule, no bar is learned and the output distribution is Gaussian, with a small standard deviation around the target value An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e066.jpg (Fig. 9A). However, after additional parameter tuning, a bar can sometimes be learned, as shown in Fig. 9B. In this case, the final output distribution is highly kurtotic, due to the receptive field. The outcome depends on the variance of the total input, which has to be large enough to start the learning process (variance was regulated by the parameter An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e067.jpg, see Methods). Most importantly, this dependence on model parameters shows that regulating the mean of the output distribution is not sufficient for reliably learning a bar and higher order moments need to be considered as well.

Figure 9
Mean firing constraint is not sufficient for reliable learning.

Interaction between IP and STDP

A good starting point for elucidating the mechanism by which the interaction between STDP and IP facilitates the discovery of an independent component is our initial problem of a single unit receiving a two dimensional input. We have previously shown in simulations that for a bounded, whitened, two dimensional input the weight vector tends to rotate towards the heavy-tailed direction in the input [16]. Here, we extend these results both analytically and in simulations. Our analysis focuses on the theoretical formulation of zero mean, unit variance inputs used for the demixing problem before and is restricted to expected changes in weights given the input and output firing rates, ignoring the time of individual spikes.

We report here only the main results of these experiments, while a detailed description is provided as supplemental information (see Text S3). Firstly, for conveniently selected pairs of input distributions, it is possible to show analytically that the weight vector rotates towards the heavy-tailed direction in the input, under the assumption that IP adaptation is faster than synaptic learning (previously demonstrated numerically in [16]). Secondly, due to the IP rule, weight changes mostly occur on the tail of the output distribution and are significantly larger for the heavy-tailed input. Namely, IP focuses learning to the heavy tailed direction in the input. When several inputs are supergaussian, the learning procedure results in the maximization of the output kurtosis, independent of the shape of the input distributions. Most importantly, we show that, for simple problems when a solution can be obtained by nonlinear PCA, our IP rule significantly speeds up learning of an independent component.

One way to understand these results could be in terms of nonlinear PCA theory. Given that for a random initial weight vector, the total input distribution is close to Gaussian, in order to enforce a sparse output, the IP has to change the transfer function in a way that ‘hides’ most of the input distribution (for example by shifting An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e072.jpg somewhere above the mean of the Gaussian). As a result, the nonlinear part of the transfer function will cover the ‘visible’ part of the input distribution, facilitating the discovery of sparse inputs by a mechanism similar to nonlinear PCA. In this light, IP provides the means to adapt the transfer function in a way that makes the nonlinear PCA particularly efficient.

Lastly, from an information-theoretic perspective, our approach can be linked to previous work on maximizing information transmission between neuronal input and output by optimizing synaptic learning [23]. This synaptic optimization procedure was shown to yield a generalization of the classic BCM rule [22]. We can show that, for a specific family of STDP implementations, which have a quadratic dependence on postsynaptic firing, IP effectively acts as a sliding threshold for BCM learning (see Text S4).

Discussion

Although ICA and related sparse coding models have been very successful in describing sensory coding in the cortex, it has been unclear how such computations can be realized in networks of spiking neurons in a biologically plausible fashion. We have presented a network of stochastically spiking neurons that performs ICA-like learning by combining different forms of plasticity. Although this is not the only attempt at computing ICA with spiking neurons, in previous models synaptic changes were not local, depending on the activity of neighboring neurons within a population [6], [7]. In this light, our model is, to our knowledge, the first to offer a mechanistic explanation of how ICA-like computation could arise by biologically plausible learning mechanisms.

In our model, IP, STDP and synaptic scaling interact to give rise to robust receptive field development. This effect does not depend on a particular implementation of STDP, but it does require an IP mechanism which enforces a sparse output distribution. Although there are very good theoretical arguments why this should be the case [11], [13], [16], the experimental evidence supporting this assumption is limited [12]. A likely explanation for this situation is the fact that it is difficult to map the experimentally observable output spikes into a probability of firing. Spike count estimates cannot be used directly, as they critically depend on the bin size. Additionally, the inter-spike interval (ISI) of an inhomogeneous Poisson process with exponentially distributed mean An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e073.jpg is indistinguishable from the ISI of a homogeneous Poisson distribution with mean An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e074.jpg. Hence, more complex statistical analyses are required for disentangling the two (see [35]).

From a computational perspective, our approach is reminiscent of several by-now classic ICA algorithms. As mentioned before, IP enforces the output distribution to be heavy-tailed, like in sparse coding [9]. Our model also shares conceptual similarities to InfoMax [29], which attempts to maximize output entropy (however, at the population level) by regulating the weights and a neuron threshold parameter. Maximizing information transmission between pre- and post-synaptic spike trains under the constraint of a fixed mean postsynaptic firing rate links our method to previous work on synaptic plasticity. A spike-based synaptic rule optimizing the above criterion [23] yields a generalization of the BCM rule [22], a powerful form of learning, which is able to discover heavy-tailed directions in the input [36], [37] and to learn Gabor receptive fields [38] in linear neurons. We have shown that, sliding threshold BCM can be viewed as a particular case of IP learning, for a specific family of STDP models.

It is interesting to think of the mechanism presented here in relation to projection pursuit [39], which tries to find good representations of high-dimensional spaces by projecting data on a lower dimensional space. The algorithm searches for interesting projection directions, a typical measure of interest being the non-Gaussianity of the distribution of data in the lower dimensional space. The difference here is that, although we do not explicitly define a contrast function maximizing kurtosis or other similar measure, our IP rule implicitly yields highly kurtotic output distributions. By sparsifying the neuron output, IP guides the synaptic learning towards the interesting (i.e. heavy-tailed) directions in the input.

From a different perspective, we can relate our method to nonlinear PCA. It is known that, for zero mean whitened data, nonlinear Hebbian learning in a rate neuron can successfully capture higher order correlations in the input [40], [41]. Moreover, it has been suggested that the precise shape of the Hebbian nonlinearity can be used for optimization purposes, for example for incorporating prior knowledge about the sources' distribution [40]. IP goes one step further in this direction, by adapting the transfer function online, during learning. From a biological perspective, there are some advantages in adapting the neuron's excitability during computation. Firstly, IP speeds up the nonlinear decorrelation of inputs. Secondly, the system gains great robustness to changes in parameters (as demonstrated in Text S2). Additionally, IP regulation plays a homeostatic role, making constraints on the input mean or second order statistics unnecessary. In the end, all the methods we have mentioned are closely related and, though conceptually similar, our approach is another distinct solution.

Our previous work was restricted to the a rate model neuron [16]. Beyond translating our results to a spiking neuron model, we have shown here that similar principles can be applied when information is encoded as spike-spike correlations, where a model relying just on firing rates would fail. It is a interesting challenge for future work to further investigate the exact mechanisms of receptive field development for different types of input encoding.

Methods

Neuron model

We consider a stochastically spiking neuron with refractoriness [23]. The model defines the neuron's instantaneous probability of firing as a function of the current membrane potential and the refractory state of the neuron,which depends on the time since its last spike. More specifically, the membrane potential is computed as An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e075.jpg where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e076.jpg is the resting potential, while the second term represents the total incoming drive to the neuron, computed as the linear summation of post-synaptic potentials evoked by incoming spikes. Here, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e077.jpg gives the strength of synapse An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e078.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e079.jpg is the time of a presynaptic spike, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e080.jpg the corresponding evoked post-synaptic potential, modeled as a decaying exponential, with time constant An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e081.jpg (for GABA-ergic synapses, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e082.jpg) and amplitude An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e083.jpg.

The refractory state of the neuron, with values in the interval An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e084.jpg, is defined as a function of the time of the last spike An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e085.jpg, namely:

equation image

where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e087.jpg gives the absolute refractory period, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e088.jpg is the relative refractory period and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e089.jpg is the Heaviside function.

The probability An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e090.jpg of the stochastic neuron firing at time An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e091.jpg is given as a function of its membrane potential and refractory state [23] An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e092.jpg where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e093.jpg is the time step of integration, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e094.jpg is a gain function, defined as: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e095.jpg Here An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e096.jpg , An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e097.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e098.jpg are model parameters, whose values are adjusted by intrinsic plasticity, as described below.

Intrinsic plasticity

Our intrinsic plasticity model attempts to maximize the mutual information between input and output, for a fixed energy budget [16], [42]. More specifically, it induces changes in neuronal excitability that lead to an exponential distribution of the instantaneous firing rate of the neuron [13]. The specific shape of the output distribution is justified from an information theoretic perspective, as the exponential distribution has maximum entropy for a fixed mean. This is true for distributions defined on the interval An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e099.jpg, but, under certain assumptions, can be a good approximation for the case where the interval is bounded, as it happens in our model due to the neuron's refractory period (see below). Optimizing information transmission under the constraint of a fixed mean is equivalent to minimizing the Kullback-Leibler divergence between the neuron's firing rate distribution and that of an exponential with mean An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e100.jpg:

equation image

with An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e102.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e103.jpg denoting the entropy and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e104.jpg the expected value. Note that the above expression assumes that the instantaneous firing rate of the neuron is proportional to An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e105.jpg, that is that An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e106.jpg. When taking into account the refractory period of the neuron, which imposes an upper-bound An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e107.jpg on the output firing rate, the maximum entropy distribution for a specific mean An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e108.jpg is a truncated exponential [43]. The deviation between the optimal exponential for the infinite and the bounded case depends on the values of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e109.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e110.jpg, but it is small in cases in which An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e111.jpg. Hence, our approximation is valid as long as the instantaneous firing rate An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e112.jpg is significantly lower than An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e113.jpg, that is when the mean firing rate of the neuron is small. In our case, we restrict An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e114.jpg. If not otherwise stated, all simulations have An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e115.jpg. Note also that the values considered here are in the range of firing rates reported for V1 neurons [44].

Computing the gradient of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e116.jpg for An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e117.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e118.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e119.jpg, and using stochastic gradient descent, the optimization process translates into the following update rules [42]:

equation image

Here, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e121.jpg is a small learning rate. Here, the instantaneous firing rate is assumed to be directly accessible for learning. Alternatively, it could be estimated based on the recent spike history.

Additionally, as a control, we have considered a simplified rule, which adjusts a single transfer function parameter in order to maintain the mean firing rate of the neuron to a constant value An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e122.jpg. More specifically, a low-pass-filtered version of the neuron firing rate is used to estimate the current mean firing rate of the neuron An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e123.jpg where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e124.jpg is the Dirac function and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e125.jpg is the time of firing of the post-synaptic neuron and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e126.jpg. Based on this estimate, the value of the parameter An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e127.jpg is adjusted as An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e128.jpg Here, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e129.jpg is the goal mean firing rate, as before and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e130.jpg is a learning rate, set such that, for a fixed Gaussian input distribution, convergence is reached as fast as for our IP rule described before (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e131.jpg).

Synaptic learning

The STDP rule implemented here considers only nearest-neighbor interactions between spikes [24]. The change in weights is determined by:

equation image

where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e133.jpg is the amplitude of the STDP change for potentiation and depression, respectively (default values An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e134.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e135.jpg), An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e136.jpg are the time scales for potentiation and depression (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e137.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e138.jpg; for learning spike-spike correlations An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e139.jpg)[19], and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e140.jpg is the time difference between the firing of the pre- and post-synaptic neuron. For the lateral inhibitory connections, the STDP learning is faster, namely An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e141.jpg. In all cases, weights are always positive and clipped to zero if they become negative.

This STDP implementation is particularly interesting as it can be shown that, under the assumption of uncorrelated or weakly correlated pre- and post-synaptic Poisson spike trains, it induces weight changes similar to a BCM rule, namely [24]:

equation image

where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e143.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e144.jpg are the firing rates of the pre- and post-synaptic neuron, respectively.

For the above expression, the fixed BCM threshold can be computed as:

equation image

which is positive when potentiation dominates depression on the short time scale, while, overall, synaptic weakening is larger than potentiation:

equation image

equation image

In some experiments we also consider the classical case of additive all-to-all STDP [26], which acts as simple Hebbian learning, the induced change in weight being proportional to the product of the pre- and post-synaptic firing rates (see [24] for comparison of different STDP implementations). The parameters used in this case are: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e148.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e149.jpg and the same time constants as for the nearest neighbor case. Additionally, the simple triplets STDP model used as an alternative BCM-like STDP implementation is described in Text S4.

Synaptic scaling

As in approaches which directly maximize kurtosis or similar measures [16], [30], [40], the weight vector is normalized: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e150.jpg, with An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e151.jpg, with weights being always positive. This value is arbitrary, as it represents a scaling factor of the total current, which can be compensated for by IP. It was selected in order to keep the final parameters close to those in [23]. Additionally, for the natural image patches the normalization was done independently for the on- and off- populations, using the same value for An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e152.jpg in each case.

In a neural population, the same normalization is applied for the lateral inhibitory connections. As before, weights do not change sign and are constrained by the An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e153.jpg norm: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e154.jpg, with An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e155.jpg.

Currently, the normalization is achieved by dividing each weight by An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e156.jpg, after the presentation of each sample. Biologically, this operation would be implemented by a synaptic scaling mechanism, which multiplicatively scales the synaptic weights to preserve the average input drive received by the neuron [20].

Setup for experiments

In all experiments, excitatory weights were initialized at random from the uniform distribution and normalized as described before. The transfer function was initialized to the parameters in [23] (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e157.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e158.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e159.jpg). Unless otherwise specified, all model parameters had the default values defined in the corresponding sections above. For all spike-based experiments, each sample was presented for a time interval An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e160.jpg, followed by the weight normalization.

For the experiments involving the rate-based model and a two-dimensional input, each sample was presented for one time step and the learning rates for IP and Hebbian learning were An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e161.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e162.jpg, respectively. In this case, the weight normalization procedure can influence the final solution. Namely, positive weights with constant An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e163.jpg norm always yield a weight vector in the first quadrant, but this limitation can be removed by a different normalization, which keeps the An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e164.jpg norm of the vector constant (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e165.jpg).

For the demixing problem, the input was generated as described for the rate-based scenario above. After the rectification, the firing rates of the input on the on- and off- channels were scaled by a factor of 20, to speed up learning. After convergence, the total weight of each channel was estimated as the sum of individual weights corresponding to that input. The resulting four-dimensional weight vector was projected back to the original two-dimensional input space using: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e166.jpg, with a sign given by that of the channel with maximum weight (positive for An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e167.jpg, negative otherwise). This procedure results by a minimum error projection of the weight vector onto the subspace defined by the constraint An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e168.jpg, see Text S1 for details.

For all variants of the bars problem, the input vector was normalized to An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e169.jpg, with An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e170.jpg defining the An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e171.jpg norm, as in [45]. Inputs were encoded using firing rates with mean An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e172.jpg, where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e173.jpg is the frequency of a background pixel and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e174.jpg gives the maximum input frequency, corresponding to a sample containing a single bar in the original bars problem.

When using the correlation-based encoding, all inputs had the same mean firing rate (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e175.jpg). Inputs corresponding to pixels in the background were uncorrelated, while inputs belonging to bars were all pairwise correlated, with a correlation coefficient An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e176.jpg. Poisson processes with such correlation structure can be generated in a computationally efficient fashion by using dichotomous Gaussian distributions [46].

When learning Gabor receptive fields, images from the van Hateren database [31] were convolved with a difference-of-gaussians filter with center and surround widths of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e177.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e178.jpg pixels, respectively. Random patches of size An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e179.jpg were selected from various positions in the images. Patches having very low contrast were discarded. The individual input patches were normalized to zero mean and unit variance, similar to the processing in [45]. The rectified values of the resulting image were mapped into a firing frequency for an on- and off-input population (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e180.jpg) and, as before, samples were presented for a duration An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e181.jpg.

For a neuronal population, input-related parameters were as for the single component, but with An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e182.jpg, to speed up learning. The initial parameters of the neuron transfer function were uniformly distributed around the default values mentioned above, with variance 0.1, 5, and 0.2 for An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e183.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e184.jpg, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e185.jpg, respectively. Additionally, the inhibitory weights were initialized at random, with no self-connections, and normalized as described before. The mutual information (MI), estimated within a window of 1000 s, was computed as An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e186.jpg, with An external file that holds a picture, illustration, etc.
Object name is pcbi.1000757.e187.jpg denoting the entropy (see [45]), applied for the average firing rate of the neurons for each input sample.

Supporting Information

Text S1

Receptive field estimation for the spike-based demixing problem

(0.04 MB PDF)

Text S2

Parameter dependence for learning one IC

(0.07 MB PDF)

Text S3

Learning with a two-dimensional input

(0.53 MB PDF)

Text S4

A link to BCM

(0.15 MB PDF)

Acknowledgments

We would like to thank Jörg Lücke and Claudia Clopath and for fruitful discussions. Additionally, we thank Jörg Lücke for providing the code for the Gabor l.m.s. fit and to Felipe Gerhard for his help with the natural images preprocessing.

Footnotes

The authors have declared that no competing interests exist.

The authors are supported by EC MEXT project PLICON and the German Federal Ministry of Education and Research (BMBF) within the Bernstein Focus: Neurotechnology through research grant 01GQ0840. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

1. Barlow H. Redundancy reduction revisited. Network. 2001;12:241–253. [PubMed]
2. Simoncelli E, Olshausen B. Natural image statistics and neural representations. Annu Rev Neuroscience. 2001;24:1193–1216. [PubMed]
3. Simoncelli E. Vision and the statistics of the visual environment. Curr Op Neurobiol. 2003;13:144–149. [PubMed]
4. Olshausen B, Field D. Sparse coding of sensory inputs. Curr Op Neurobiol. 2004;14:481–487. [PubMed]
5. Hyvärinen A, Karhunen J, Oja E. Independent Component Analysis. 2001. Wiley Series on Adaptive and Learning Systems for Signal Processing, Communication and Control.
6. Klampfl S, Legenstein R, Maass W. Spiking neurons can learn to solve information bottleneck problems and extract independent components. Neural Computation. 2009;21:911–959. [PubMed]
7. Parra L, Beck J, Bell A. On the maximization of information flow between spiking neurons. Neural Comp. 2009;21:1–19. [PubMed]
8. Clopath C, Büsing L, Vasilaki E, Gerstner W. Connectivity reflects coding: a model of voltage-based stdp with homeostasis. Nature Neuroscience 2010 [PubMed]
9. Olshausen B, Field D. Sparse coding with an overcomplete basis set: A strategy employed by V1? Vision Research. 1997;37:3311–3325. [PubMed]
10. Olshausen B, Field D. Emergence of simple-cell receptive field properties by learning a sparse code for natural images. Nature. 1996;381:607–609. [PubMed]
11. Lennie P. The cost of neural computation. Current Biology. 2003;13:493–497. [PubMed]
12. Baddeley R, Abbott L, Booth M, Sengpiel F, Freeman T, et al. Responses of neurons in primary and inferior temporal visual cortices to natural scenes. Proceedings Biological Sciences. 1997;264:1775–1783. [PMC free article] [PubMed]
13. Stemmler M, Koch C. How voltage-dependent conductances can adapt to maximize the information encoded by neuronal firing rate. Nature Neuroscience. 1999;2:521–527. [PubMed]
14. Zhang W, Linden D. The other side of the engram: experience-dependent changes in neuronal intrinsic excitability. Nature Reviews Neuroscience. 2003;4:885–900. [PubMed]
15. Cudmore R, Turrigiano G. Long-term potentiation of intrinsic excitability in LV visual cortical neurons. J Neurophysiol. 2004;92:341–348. [PubMed]
16. Triesch J. Synergies between intrinsic and synaptic plasticity mechanisms. Neural Computation. 2007;19:885–909. [PubMed]
17. Gerstner W, Kempter R, van Hemmen J, Wagner H. A neuronal learning rule for sub-millisecond temporal coding. Nature. 1996;383:76–78. [PubMed]
18. Markram H, Lübke J, Frotscher M, Sackmann B. Regulation of synaptic efficacy by coincidence of postsynaptic APS and EPSPS. Science. 1997;275:213–215. [PubMed]
19. Bi G, Poo M. Synaptic modifications in cultured hippocampal neurons: Dependence on spike timing, synaptic strength, and postsynaptic cell type. J Neurosci. 1998;18:10464–10472. [PubMed]
20. Turrigiano G, Nelson S. Homeostatic plasticity in the developing nervous system. Nature Reviews Neuroscience. 2004;5:97–107. [PubMed]
21. Foldiák P. Forming sparse representations by local anti-Hebbian learning. Biological Cybernetics. 1990;64:165–170. [PubMed]
22. Bienenstock E, Cooper L, Munro P. Theory for the development of neuron selectivity: Orientation specificity and binocular interactions in visual cortex. Journal of Neuroscience. 1982;2:32–48. [PubMed]
23. Toyoizumi T, Pfister J, Aihara K, Gerstner W. Generalized Bienenstock-Cooper-Munro rule for spiking neurons that maximizes information transmission. PNAS. 2005;102:5239–5244. [PubMed]
24. Izhikevich E, Desai N. Relating STDP to BCM. Neural Computation. 2003;15:1511–1523. [PubMed]
25. Martinez-Conde S, Macknik S, Hubel D. The role of fixational eye movements in visual perception. Nature Reviews Neuroscience. 2004;5:229–240. [PubMed]
26. Song S, Miller K, Abbott L. Competitive Hebbian learning through spike-timing-dependent synaptic plasticity. Nature Neurosci. 2000;3:919–926. [PubMed]
27. Pfister J, Gerstner W. Triplets of spikes in a model of spike timing-dependent plasticity. Journal of Neuroscience. 2006;26:9673–9682. [PubMed]
28. Lücke J, Sahani M. Maximal causes for non-linear component extraction. Journal of Machine Learning Research. 2008;9:1227–1267.
29. Bell A, Sejnowski T. The independent components of natural scenes are edge filters. Vision Research. 1997;37:3327–3338. [PMC free article] [PubMed]
30. Hyvärinen A. Survey on Independent Component Analysis. Neural Computing Surveys. 1999;2:94–128.
31. van Hateren J. Independent component filters of natural images compared with simple cells in primary visual cortex. Proceedings of the Royal Society B: Biological Sciences. 1998;265:359–366. [PMC free article] [PubMed]
32. Falconbridge M, Stamps R, Badcock D. A simple Hebbian/anti-Hebbian network learns the sparse, independent components of natural images. Neural Computation. 2006;18:415–429. [PubMed]
33. Weber C, Triesch J. A sparse generative model of V1 simple cells with intrinsic plasticity. Neural Computation. 2008;20:1261–1284. [PubMed]
34. Lücke J. Receptive field self-organization in a model of the fine-structure in V1 cortical columns. Neural Computation. 2009;21:2805–2845. [PubMed]
35. Cox D. Some statistical methods connected with series of events. Journal of the Royal Statistical Society, Series B. 1955;17:129–164.
36. Intrator N, Cooper L. Objective function formulation of the BCM theory of visual cortical plasticity: Statistical connections, stability conditions. Neural Networks. 1992;5:3–17.
37. Intrator N. Springer-Verlag, editor. Neuronal goals: efficient coding and coincidence detection. ICONIP Hong Kong: Progress in Neural Information Processing. 1998. pp. 29–34. volume 1.
38. Blais B, Intrator N, Shouval H, Cooper L. Receptive field formation in natural scene environments: comparison of single cell learning rules. Neural Computation. 1998;10:1797–1813. [PubMed]
39. Huber P. Projection pursuit. The Annals of Statistics. 1985;13:435–475.
40. Hyvärinen A. One-unit contrast functions for independent component analysis: a statistical analysis. Neural Networks for Signal Processing. 1997:388–397.
41. Oja E. The nonlinear PCA learning rule in independent component analysis. Neurocomputing. 1997;17:25–46.
42. Joshi P, Triesch J. Rules for information-maximization in spiking neurons using intrinsic plasticity. Proc. IJCNN. 2009:1456–1461.
43. Kapur N. Maximum Entropy Models in Science and Engineering. Wiley-Interscience; 1990.
44. Olshausen B, Simoncelli E. Natural image statistics and neural representation. Annu Rev Neuroscience. 2001;24:1193–1216. [PubMed]
45. Butko N, Triesch J. Learning sensory representations with intrinsic plasticity. Neurocomputing. 2007;70
46. Macke J, Berens P, Ecker A, Tolias A, Bethge M. Generating spike trains with specified correlation coefficients. Neural Comp. 2009;21:397–423. [PubMed]

Articles from PLoS Computational Biology are provided here courtesy of Public Library of Science