Definition of the network model

Our model consists of a network of spiking neurons, arranged in a WTA circuit, which is one of the most frequently studied connectivity patterns (or network motifs) of cortical microcircuits [16]. The input of the circuit is represented by the excitatory neurons 👁 Image
. This input projects to a population of excitatory neurons 👁 Image
that are arranged in a WTA circuit (see Fig. 1). We model the effect of lateral inhibition, which is the competition mechanism of a WTA circuit [33], by a common inhibitory signal 👁 Image
that is fed to all 👁 Image
neurons and in turn depends on the activity of the 👁 Image
neurons. Evidence for such common local inhibitory signals for nearby neurons arises from numerous experimental results, see e.g. [16], [34]–[36]. We do not a priori impose a specific functional relationship between the common inhibition signal and the excitatory activity. Instead we will later derive necessary conditions for this relationship, and propose a mechanism that we use for the experiments.

The individual units 👁 Image
are modeled by a simplified Spike Response Model [37] in which the membrane potential is computed as the difference between the excitatory input 👁 Image
and the common inhibition term 👁 Image
. 👁 Image
sums up the excitatory inputs from neurons 👁 Image
as👁 Image
(2)👁 Image
models the EPSPs evoked by spikes of the presynaptic neuron 👁 Image
, and 👁 Image
models the intrinsic excitability of the neuron 👁 Image
. In order to simplify our analysis we assume that the EPSP can be modeled as a step function with amplitude 👁 Image
, i.e., 👁 Image
it takes on the value 1 in a finite time window of length 👁 Image
after a spike and is zero before and afterwards. Further spikes within this time window do not contribute additively to the EPSP, but only extend the time window during which the EPSP is in the high state. We will later show how to extend our results to the case of realistically shaped and additive EPSPs.

We use a stochastic firing model for 👁 Image
, in which the firing probability depends exponentially on the membrane potential, i.e.,👁 Image
(3)which is in good agreement with most experimental data [38]. We can thus model the firing behavior of every neuron 👁 Image
in the WTA as an independent inhomogeneous Poisson process whose instantaneous firing rate is given by 👁 Image
.

In order to understand how this network model generates samples from a probability distribution, we first observe that the combined firing activity of the neurons 👁 Image
in the WTA circuit is simply the sum of the 👁 Image
independent Poisson processes, and can thus again be modeled as an inhomogeneous Poisson process with rate 👁 Image
. Furthermore, in any infinitesimally small time interval 👁 Image
, the neuron 👁 Image
spikes with probability 👁 Image
. Thus, if we know that at some point in time 👁 Image
, i.e. within 👁 Image
, one of the neurons 👁 Image
produces an output spike, the conditional probability 👁 Image
that this spike originated from neuron 👁 Image
can be expressed as👁 Image
(4)Every single spike from the WTA circuit can thus be seen as an independent sample from the instantaneous distribution in Eq. (4) at the time of the spike. Although the instantaneous firing rate of every neuron directly depends on the value of the inhibition 👁 Image
, the relative proportion of the rate 👁 Image
to the total WTA firing rate 👁 Image
is independent of the inhibition, because all neurons receive the same inhibition signal 👁 Image
. Note that 👁 Image
determines only the value of the sample at time 👁 Image
, but not the time point at which a sample is created. The temporal structure of the sampling process depends only on the overall firing rate 👁 Image
.

This implementation of a stochastic WTA circuit does not constrain in any way the kind of spike patterns that can be produced. Every neuron fires independently according to a Poisson process, so it is perfectly possible (and sometimes desirable) that there are two or more neurons that fire (quasi) simultaneously. This is no contradiction to the above theoretical argument of single spikes as samples. There we assumed that there was only one spike at a time inside a time window, but since we assumed these windows to be infinitesimally small, the probability of two spikes occurring exactly at the same point in continuous time is zero.

Example 1: Learning of probabilistic models with STDP

Fig. 3 demonstrates the emergence of Bayesian computation in the generic network motif of Fig. 1A in a simple example. Spike inputs 👁 Image
(top row of Fig. 3D) are generated through four different hidden processes (associated with four different colors). Each of them is defined by a Gauss distribution over a 2D pixel array with a different center, which defines the probability of every pixel to be on. Spike trains encode the current value of a pixel by a firing rate of 25 Hz or 0 Hz for 40 ms. Each pixel was encoded by two input neurons 👁 Image
via population coding, exactly one of them had a firing rate of 25 Hz for each input image. A 10 ms period without firing separates two images in order to avoid overlap of EPSPs for input spikes belonging to different input images.

After unsupervised learning with STDP for 500 s (applied to continuous streams of spikes as in panel D of Fig. 3) the weight vectors shown in Fig. 3F (projected back into the virtual 2D input space) emerged for the four output neurons 👁 Image
, demonstrating that these neurons had acquired internal models for the four different processes that were used to generate inputs. The four different processes for generating the underlying 2D input patterns had been used with different prior probabilities (👁 Image
). Fig. 3G shows that this imbalance resulted in four different priors 👁 Image
encoded in the biases 👁 Image
of the neurons 👁 Image
. When one compares the unequal sizes of the colored areas in Fig. 3H with the completely symmetric internal models (or likelihoods) of the four neurons shown in panel F, one sees that their firing probability approximates a posterior over hidden causes that results from multiplying their learned likelihoods with their learned priors. As a result, the spike output becomes sparser, and almost all neurons only fire when the current input spikes are generated by that one of the four hidden processes on which they have specialized (Fig. 3D, bottom row). In Fig. 3I the performance of the network is quantified over time by the normalized conditional entropy 👁 Image
, where 👁 Image
is the correct hidden cause of each input image 👁 Image
in the training set, and 👁 Image
denotes the discrete random variable defined by the firing probabilities of output neurons 👁 Image
for each image under the currently learned model. Low conditional entropy indicates that each neuron learns to fire predominantly for inputs from one class. Fig. 3E as well as the dashed blue line in Fig. 3I show that the learning process is improved when a common background oscillation at 20 Hz is superimposed on the firing rate of input neurons and the membrane potential of the output neurons, while keeping the average input and output firing rates constant. The reason is that in general it may occur that an output neuron 👁 Image
receives during its integration time window (40 ms in this example) no information about the value of a pixel (because neither the neuron 👁 Image
that has a high firing rate for 40 ms if this pixel is black, nor the associated neuron 👁 Image
that has a high firing rate if this pixel is white fire during this time window). A background oscillation reduces the percentage of such missing values by driving presynaptic firing times together (see top row of Fig. 3E). Note that through these oscillations the overall output firing rate 👁 Image
fluctuates strongly, but since the same oscillation is used consistently for all four types of patterns, the circuit still learns the correct distribution of inputs.

This task had been chosen to become very fast unsolvable if many pixel values are missing. Many naturally occurring input distributions, like the ones addressed in the subsequent computer experiments, tend to have more redundancy, and background oscillations did not improve the learning performance for those.

STDP approximates Expectation Maximization

In this section we will develop the link between the unsupervised learning of the generative probabilistic model in Fig. 1B and the learning effect of STDP as defined in our spiking network model in Fig. 1A. Starting from a learning framework derived from the concept of Expectation Maximization [31], we show that the biologically plausible STDP rule from Fig. 2 can naturally approximate a stochastic, online version of this optimization algorithm. We call this principle SEM (spike-based EM).

SEM can be viewed as a bootstrapping procedure. The relation between the firing probabilities of the neurons within the WTA circuit and the continuous updates of the synaptic weights with our STDP rule in Eq. (5) drive the initially random firing of the circuit in response to an input 👁 Image
towards learning the correct generative model of the input distribution. Whenever a neuron 👁 Image
fires in response to 👁 Image
, the STDP rule increases the weights 👁 Image
of synapses from those presynaptic neurons 👁 Image
that had fired shortly before 👁 Image
. In absence of a recent presynaptic spike from 👁 Image
the weight 👁 Image
is decreased. As a consequence, when next a pattern similar to 👁 Image
is presented, the probability for the same 👁 Image
to fire and further adapt its weights, is increased. Since 👁 Image
becomes more of an “expert” for one subclass of input patterns, it actually becomes less likely to fire for non-matching patterns. The competition in the WTA circuit ensures that other 👁 Image
-neurons learn to specialize for these different input categories.

In the framework of Expectation Maximization, the generation of a spike in a 👁 Image
-neuron creates a sample from the currently encoded posterior distribution of hidden variables, and can therefore be viewed as the stochastic Expectation, or 👁 Image
-step. The subsequent application of STDP to the synapses of this neuron can be understood as an approximation of the Maximization, or 👁 Image
-step. The online learning behavior of the network can be understood as a stochastic online EM algorithm.

The role of the inhibition

We have previously shown that the output spikes of the WTA circuit represent samples from the posterior distribution in Eq. (11), which only depends on the ratios between the membrane potentials 👁 Image
. The rate at which these samples are produced is the overall firing rate 👁 Image
of the WTA circuit and can be controlled by modifying the common inhibition 👁 Image
of the neurons 👁 Image
.

Although any time-varying output firing rate 👁 Image
produces correct samples from the posterior distribution in Eq. (11) of 👁 Image
, for learning we also require that the input patterns 👁 Image
observed at the spike times are unbiased samples from the true input distribution 👁 Image
. If this is violated, some patterns coincide with a higher 👁 Image
, and thus have a stronger influence on the learned synaptic weights. In Methods we formally show that 👁 Image
acts as a multiplicative weighting of the current input 👁 Image
, and so the generative model will learn a slightly distorted input distribution.

An unbiased set of samples can be obtained if 👁 Image
is independent of the current input activation 👁 Image
, e.g. if 👁 Image
is constant. This could in theory be achieved if we let 👁 Image
depend on the current values of the membrane potentials 👁 Image
, and set 👁 Image
. Such an immediate inhibition is commonly assumed in rate-based soft-WTA models, but it seems implausible to compute this in a spiking neuronal network, where only spikes can be observed, but not the presynaptic membrane potentials.

However, our results show that a perfectly constant firing rate is not a prerequisite for convergence to the right probabilistic model. Indeed we can show that it is sufficient that 👁 Image
and 👁 Image
are stochastically independent, i.e. 👁 Image
is not correlated to the appearance of any specific value of 👁 Image
. Still this might be difficult to achieve since the firing rate 👁 Image
is functionally linked to the input 👁 Image
by 👁 Image
, but it clarifies the role of the inhibition 👁 Image
as de-correlating 👁 Image
from the input 👁 Image
, at least in the long run.

One possible biologically plausible mechanism for such a decorrelation of 👁 Image
and 👁 Image
is an inhibitory feedback from a population of neurons that is itself excited by the neurons 👁 Image
. Such WTA competition through lateral inhibition has been studied extensively in the literature [16], [33]. In the implementation used for the experiments in this paper every spike from the 👁 Image
-neurons causes an immediate very strong inhibition signal that lasts longer than the refractory period of the spiking neuron. This strong inhibition decays exponentially and is overlaid by a noise signal with high variability that follows an Ornstein-Uhlenbeck process (see “Inhibition Model in Computer Simulations” in Methods). This will render the time of the next spike of the system almost independent of the value of 👁 Image
.

It should also be mentioned that a slight correlation between 👁 Image
and 👁 Image
may be desirable, and 👁 Image
might also be externally modulated (for example through attention, or neuromodulators such as Acetylcholin), as an instrument of selective input learning. This might lead e.g. to slightly higher firing rates for well-known inputs (high 👁 Image
), or salient inputs, as opposed to reduced rates for unknown arbitrary inputs. In general, however, combining online learning with a sampling rate 👁 Image
that is correlated to 👁 Image
may lead to strange artifacts and might even prohibit the convergence of the system due to positive feedback effects. A thorough analysis of such effects and of possible learning mechanisms that cope with positive feedback effects is the topic of future research.

Our theoretical analysis sheds new light on the requirements for inhibition in spiking WTA-like circuits to support learning and Bayesian computation. Inhibition does not only cause competition between the excitatory neurons, but also regulates the overall firing rate 👁 Image
of the WTA circuit. Variability in 👁 Image
does not influence the performance of the circuit, as long as there is no systematic dependence between the input and 👁 Image
.

Continuous-time interpretation with realistically shaped EPSPs

In our previous analysis we have assumed a simplified non-additive step-function model for the EPSP. This allowed us to describe all input evidence within the last time window of length 👁 Image
by one binary vector 👁 Image
, but required us to assume that no two neurons within the same group 👁 Image
fired within that period. We will now give an intuitive explanation to show that this restriction can be dropped and present an interpretation for additive biologically plausibly shaped EPSPs as inference in a generative model.

The postsynaptic activation 👁 Image
under an additive EPSPs is given by the convolution👁 Image
(17)where 👁 Image
describes an arbitrarily shaped kernel, e.g. an 👁 Image
-shaped EPSP function which is the difference of two exponential functions (see [37]) with different time constants. We use 1 ms for the rise and 15 ms for the decay in our simulations. 👁 Image
replaces 👁 Image
in Eq. (2) in the computation of the membrane potential 👁 Image
of our model neurons. We can still understand the firing of neurons in the WTA circuit according to the relative firing probabilities 👁 Image
in Eq. (4) as Bayesian inference. To see this, we imagine an extension of the generative probabilistic model 👁 Image
in Fig. 1B, which contains multiple instances of 👁 Image
, exactly one for every input spike from all input neurons 👁 Image
. For a fixed common hidden cause 👁 Image
, all instances of 👁 Image
are conditionally independent of each other, and have the same conditional distributions for each 👁 Image
(see Methods for the full derivation of the extended probabilistic model). According to the definition in Eq. (8) of the population code every input spike represents evidence that 👁 Image
in an instance 👁 Image
should take on a certain value. Since every spike contributes only to one instance, any finite input spike pattern can be interpreted as valid evidence for multiple instances of inputs 👁 Image
.

The inference of a single hidden cause 👁 Image
in such extended graphical model from multiple instances of evidence is relatively straightforward: due to the conditional independence of different instances, we can compute the input likelihood for any hidden cause simply as the product of likelihoods for every single evidence. Inference thus reduces to counting how often every possible evidence occurred in all instances 👁 Image
, which means counting the number of spikes of every 👁 Image
. Since single likelihoods are implicitly encoded in the synaptic weights 👁 Image
by the relationship 👁 Image
, we can thus compute the complete input likelihood by adding up step-function like EPSPs with amplitudes corresponding to 👁 Image
. This yields correct results, even if one input neuron spikes multiple times.

In the above model, the timing of spikes does not play a role. If we want to assign more weight to recent evidence, we can define a heuristic modification of the extended graphical model, in which contributions from spikes to the complete input log-likelihood are linearly interpolated in time, and multiple pieces of evidence simply accumulate. This is exactly what is computed in 👁 Image
in Eq. (17), where the shape of the kernel 👁 Image
defines how the contribution of an input spike at time 👁 Image
evolves over time. Defining 👁 Image
as the weight for the evidence of the assignment of 👁 Image
to value 👁 Image
, it is easy to see (and shown in detail in Methods) that the instantaneous output distribution 👁 Image
represents the result of inference over causes 👁 Image
, given the time-weighted evidences of all previous input spikes, where the weighting is done by the EPSP-function 👁 Image
. Note that this evidence weighting mechanism is not equivalent to the much more complex mechanism for inference in presence of uncertain evidence, which would require more elaborate architectures than our feed-forward WTA-circuit. In our case, past evidence does not become uncertain, but just less important for the inference of the instantaneous hidden cause 👁 Image
.

We can analogously generalize the spike-triggered learning rule in Eq. (5) for continuous-valued input activations 👁 Image
according to Eq. (17):👁 Image
(18)The update of every weight 👁 Image
is triggered when neuron 👁 Image
, i.e. the postsynaptic neuron, fires a spike. The shape of the LTP part of the STDP curve is determined by the shape of the EPSP, defined by the kernel function 👁 Image
. The positive part of the update in Eq. (18) is weighted by the value of 👁 Image
at the time of firing the postsynaptic spike. Negative updates are performed if 👁 Image
is close to zero, which indicates that no presynaptic spikes were observed recently. The complex version of the STDP curve (blue dashed curve in Fig. 1B), which resembles more closely to the experimentally found STDP curves, results from the use of biologically plausible 👁 Image
-shaped EPSPs. In this case, the LTP window of the weight update decays with time, following the shape of the 👁 Image
-function. This form of synaptic plasticity was used in all our experiments. If EPSPs accumulate due to high input stimulation frequencies, the resulting shape of the STDP curve becomes even more similar to previously observed experimental data, which is investigated in detail in the following section.

The question remains, how this extension of the model and the heuristics for time-dependent weighting of spike contributions affect the previously derived theoretical properties. Although the convergence proof does not hold anymore under such general conditions we can expect (and show in our Experiments) that the network will still show the principal behavior of EM under fairly general assumptions on the input: we have to assume that the instantaneous spike rate of every input group 👁 Image
is not dependent on the value of 👁 Image
that it currently encodes, which means that the total input spike rate must not depend on the hidden cause 👁 Image
. Note that this assumption on every input group is identical to the desired output behavior of the WTA circuit according to the conditions on the inhibition as derived earlier. This opens up the possibility of building networks of recursively or hierarchically connected WTA circuits. Note also that the grouping of inputs into different 👁 Image
is only a notational convenience. The neurons in the WTA circuit do not have to know which inputs are from the same group, neither for inference nor for learning, and can thus treat all input neurons equally.

Relationship to experimental data on synaptic plasticity

In biological STDP experiments that induce pairs of pre- and post-synaptic spikes at different time delays, it has been observed that the shape of the plasticity curve changes as a function of the repetition frequency for those spike pairs [40]. The observed effect is that at very low frequencies no change or only LTD occurs, a “classical” STDP window with timing-dependent LTD and LTP is observed at intermediate frequencies around 20 Hz, and at high frequencies of 40 Hz or above only LTP is observed, independently of which spikes comes first.

Although our theoretical model does not explicitly include a stimulation-frequency dependent term like other STDP models (e.g. [53]), we can study empirically the effect of a modification of the frequency of spike-pairing. We simulate this for a single synapse, at which we force pre- and post-synaptic spikes with varying time differences 👁 Image
, and at fixed stimulation frequencies 👁 Image
of either 1 Hz, 20 Hz, or 40 Hz. Modeling EPSPs as 👁 Image
-kernels with time constants of 1 ms for the rise and 15 ms for the decay, we obtain the low-pass filtered signals 👁 Image
as in Eq. (17), which grow as EPSPs start to overlap at higher stimulation frequencies. At the time of a post-synaptic spike we compute the synaptic update according to the rule in Eq. (18), but keep both the weight and the learning rate fixed (at 👁 Image
) to distinguish timing-dependent from weight-dependent effects.

In Fig. 4A we observe that, as expected, at low stimulation frequencies (1 Hz) the standard shape of the complex STDP rule in Eq. (18) from Fig. 2 is recovered, since there is no influence from previous spikes. The shift towards pure LTD that is observed in biology [40] would require an additional term that depends on postsynaptic firing rates like in [53], and is a topic of future research. However, note that in biology this shift to LTD was observed only in paired recordings, neglecting the cooperative effect of other synapses, and other studies have also reported LTP at low stimulation frequencies [43]. At higher stimulation frequencies (20 Hz in Fig. 4B) the EPSPs from different pre-synaptic spikes start to overlap, which results in larger 👁 Image
compared with isolated pre-synaptic spikes. We also see that the LTD part of the STDP window becomes timing-dependent (due to overlapping EPSPs), and thus the shape of the STDP curve becomes similar to standard models of STDP and observed biological data [43], [54]. For even higher stimulation frequencies the STDP window shifts more and more towards LTP (see Fig. 4B and C). This is in good accordance with observations in biology [40]. Also in agreement with biological data, the minimum of the update occurs around 👁 Image
, because there the new 👁 Image
-kernel EPSP is not yet effective, and the activation due to previous spikes has decayed maximally.

Another effect that is observed in hippocampal synapses when two neurons are stimulated with bursts, is that the magnitude of LTP is determined mostly by the amount of overlap between the pre- and post-synaptic bursts, rather than the exact timing of spikes [55]. In Fig. 4D we simulated this protocol with our continuous-time SEM rule for different onset time-differences of the bursts, and accumulated the synaptic weight updates in response to 50 Hz bursts of 5 pre-synaptic and 4 post-synaptic spikes. We performed this experiment for the same onset time differences used in Fig. 3 of [55], and found qualitatively similar results. For long time-differences, when EPSPs have mostly decayed, we observed an LTD effect, which was not observed in biology, but can be attributed to differences in synaptic time constants between biology and simulation.

These results suggest that our STDP rule derived from theoretical principles exhibits several of the key properties of synaptic plasticity observed in nature, depending on the encoding of inputs. This is quite remarkable, since these properties are not explicitly part of our learning rule, but rather emerge from a simpler rule with strong theoretical guarantees. Other phenomenological [56], [57] or mechanistic models of STDP [58] also show some of these characteristics, but come without such theoretical properties. The functional consequence of reproducing such key biological characteristics of STDP is that our new learning rule also exhibits most of the key functional properties of STDP, like e.g. strengthening synapses of inputs that are causally involved in firing the postsynaptic neuron, while pruning the connections that do not causally contribute to postsynaptic firing [10], [13]. At low and intermediate firing rates our rule also shifts the onset of postsynaptic firing towards the start of repeated spike patterns [49], [50], [59], while depressing synapses that only become active for a pattern following the one for which the post-synaptic neuron is responsive. If patterns change quickly, then the stronger depression for presynaptic spikes with small 👁 Image
in Fig. 4B enhances the capability of the WTA to discriminate such patterns. With simultaneous high frequency stimulation (Fig. 4C and D) we observe that only LTP occurs, which is due to the decay of EPSPs not being fast enough to allow depression. In this scenario, the learning rule is less sensitive to timing, and rather becomes a classical Hebbian measure of correlations between pre- and post-synaptic firing rates. However, since inputs are encoded in a population code we can assume that the same neuron is not continuously active throughout, and so even at high firing rates for active input neurons, the synapses that are inactive during postsynaptic firing will still be depressed, which means that convergence to an equilibrium value is still possible for all synapses.

It is a topic of future research which effects observed in biology can be reproduced with more complex variations of the spike-based EM rule that are also dependent on postsynaptic firing rates, or whether existing phenomenological models of STDP can be interpreted in the probabilistic EM framework. In fact, initial experiments have shown that several variations of the spike-based EM rule can lead to qualitatively similar empirical results for the learned models in tasks where the input spike trains are Poisson at average or high rates over an extended time window (such as in Fig. 3). These variations include weight-dependent STDP rules that are inversed in time, symmetrical in time, or have both spike timing-dependent LTD and LTP. Such rules can converge towards the same equilibrium values as the typical causal STDP rule. However, they will behave differently if inputs are encoded through spatio-temporal spike patterns (as in Example 4: Detection of Spatio-Temporal Spike Patterns). Further variations can include short-term plasticity effects for pre-synaptic spikes, as observed and modeled in [60], which induce a stimulation-frequency dependent reduction of the learning rate, and could thus serve as a stabilization mechanism.

Spike-timing dependent LTD

Current models of STDP typically assume a “double-exponential” decaying shape of the STDP curve, which was first used in [54] to fit experimental data. This is functionally different from the shape of the complex STDP curve in Fig. 2 and Eq. (5), where the LTD part is realized by a constant timing-independent offset.

Although not explicitly covered by the previously presented theory of SEM, the same analytical tools can be used to explain functional consequences of timing-dependent LTD in our framework. Analogous to our approach for the standard SEM learning rule, we develop (in Methods) an extension of the simple step-function STDP rule from Fig. 2 with timing-dependent LTD, which is easier to analyze. We then generalize these results towards arbitrarily shaped STDP curves. The crucial result is that as long as the spike-timing dependent LTD rule retains the characteristic inversely-exponential weight-dependent relationship between the strengths of LTP and LTD that was introduced for standard SEM in Eq. (5), an equilibrium property similar to Eq. (6) still holds (see Methods for details). Precisely speaking, the new equilibrium will be at the difference between the logarithms of the average presynaptic spiking probabilities before and after the postsynaptic spike. This shows that spike-timing dependent LTD also yields synaptic weights that can be interpreted in terms of log-probabilities, which can thus be used for inference.

The new rule emphasizes contrasts between the current input pattern and the immediately following activity. Still, the results of the new learning rule and the original rule from Eq. (5) in our experiments are qualitatively similar. This can be explained from a stochastic learning perspective: at any point in time the relative spiking probabilities of excitatory neurons in the WTA circuit in Eq. (4) depend causally on the weighted sums of preceding presynaptic activities 👁 Image
. However, they clearly do not depend on future presynaptic activity. Thus, the postsynaptic neuron will learn through SEM to fire for increasingly similar stochastic realizations of presynaptic input 👁 Image
, whereas the presynaptic activity pattern following a postsynaptic spike will become more variable. In the extreme case where patterns are short and separated by noise, there will be no big difference between input patterns following firing of any of the WTA neurons, and so their relevance for the competition will become negligible.

Experimental evidence shows that the time constants of the LTP learning window are usually smaller than the time constants of the LTD window ([40], [60]), which will further enhance the specificity of the LTP learning as opposed to the LTD part that computes the average over a longer window.

Note that the exponential weight dependence of the learning rule implies a certain robustness towards linearly scaling LTP or LTD strengths, which only leads to a constant offset of the weights. Assuming that the offset is the same for all synapses, this does not affect firing probabilities of neurons in a WTA circuit (see Methods “Weight offsets and positive weights”).

Example 2: Learning of probabilistic models for orientation selectivity

We demonstrated in this computer experiment the emergence of orientation selective cells 👁 Image
through STDP in the WTA circuit of Fig. 1A when the spike inputs encode isolated bars in arbitrary orientations. Input images were generated by the following process: Orientations were sampled from a uniform distribution, and lines of 7 pixels width were drawn in a 28×28 pixel array. We added noise to the stimuli by flipping every pixel with a 👁 Image
chance, see Fig. 5A. Finally, a circular mask was applied to the images to avoid artifacts from image corners. Spikes trains 👁 Image
were encoded according to the same population coding principle described in the previous example Fig. 3, in this case using a Poisson firing rate of 20 Hz for active units.

After training with STDP for 200 s, presenting 👁 Image
different images, the projection of the learned weight vectors back into the 2D input space (Fig. 5B) shows the emergence of 10 models with different orientations, which cover the possible range of orientations almost uniformly. When we plot the strongest responding neuron as a function of orientation (Fig. 5C, D), measured by the activity in response to 360 noise-free images of oriented bars in 👁 Image
steps, we can see no structure in the response before learning (Fig. 5C). However, after unsupervised learning, panel D clearly shows the emergence of continuous, uniformly spaced regions in which one of the 👁 Image
neurons fires predominantly. This can also be seen in the firing behavior in response to the input spike trains in Fig. 5E, which result from the example images in panel A. Fig. 5F shows that the output neurons initially fire randomly in response to the input, and many different 👁 Image
neurons are active for one image. In contrast, the responses after learning in panel G are much sparser, and only occasionally multiple neurons are active for one input image, which is the case when the angle of the input image is in between the preferred angles of two output neurons, and therefore multiple models have a non-zero probability of firing.

In our experiment the visual input consisted of noisy images of isolated bars, which illustrates learning of a probabilistic model in which a continuous hidden cause (the orientation angle) is represented by a population of neurons, and also provides a simple model for the development of orientation selectivity. It has previously been demonstrated that similar Gabor-like receptive field structures can be learned with a sparse-coding approach using patches of natural images as inputs [61]. The scenario considered here is thus substantially simplified, since we do not present natural but isolated stimuli. However, it is worth noting that experimental studies have shown that (in mice and ferret) orientation selectivity, but not e.g. direction selectivity, exists in V1 neurons even before eye opening [62], [63]. This initial orientation selectivity develops from innate mechanisms and from internally generated inputs during this phase [63], e.g. retinal waves, which have different, and very likely simpler statistics than natural stimuli. Our model shows that a WTA circuit could learn orientation selectivity from such simple bar-like inputs, but does not provide an alternative explanation to the results of studies like [61] using natural image stimuli. Although beyond the scope of this paper, we expect that later shaping of selectivity through exposure to natural visual experience would not alter the receptive fields by much, since the neurons have been primed to spike (and thereby trigger plasticity) only in response to a restricted class of local features.

Example 3: Emergent discrimination of handwritten digits through STDP

Spike-based EM is a quite powerful learning principle, as we demonstrate in Fig. 6 through an application to a computational task that is substantially more difficult than previously considered tasks for networks of spiking neurons: We show that a simple network of spiking neurons can learn without any supervision to discriminate handwritten digits from the MNIST benchmark dataset [64] consisting of 70,000 samples (30 are shown in Fig. 6A). This is one of the most frequently used benchmark tasks in machine learning. It has mostly been used to evaluate supervised or semi-supervised machine learning algorithms [27], [65], or to evaluate unsupervised feature learning approaches [66], [67]. Although the MNIST dataset contains labels (the intended digit) for each sample of a handwritten digit, we deleted these labels when presenting the dataset to the neural circuit of Fig. 1A, thereby forcing the 👁 Image
neurons on the output layer to self-organize in a completely unsupervised fashion. Each sample of a handwritten digit was encoded by 708 spike trains over 40 ms (and 10 ms periods without firing between digits to avoid overlap of EPSPs between images), similarly as for the task of Fig. 3. Each pixel was represented by two input neurons 👁 Image
, one of which produced a Poisson spike train at 40 Hz during these 40 ms. This yielded usually at most one or two spikes during this time window, demonstrating that the network learns and computes with information that is encoded through spikes, rather than firing rates. After 500 s of unsupervised learning by STDP almost all of the output neurons fired more sparsely, and primarily for handwritten samples of just one of the digits (see Fig. 6E).

The application to the MNIST dataset had been chosen to illustrate the power of SEM in complex tasks. MNIST is one of the most popular benchmarks in machine learning, and state-of-the-art methods achieve classification error rates well below 👁 Image
. The model learned by SEM can in principle also be used for classification, by assigning each neuron to the class for which it fires most strongly. However, since this is an unsupervised method, not optimized for classification but for learning a generative model, the performance is necessarily worse. We achieve an error rate of 👁 Image
on the 10-digit task on a previously unseen test set. This compares favorably to the 👁 Image
error that we obtained with a standard machine learning approach that directly learned the mixture-of-multinomials graphical model in Fig. 1B with a batch EM algorithm. This control experiment was not constrained by a neural network architecture or biologically plausible learning, but instead mathematically optimized the parameters of the model in up to 200 iterations over the whole training set. The batch method achieves a final conditional entropy of 👁 Image
, which is slightly better than the 👁 Image
final result of the SEM approach, and shows that better performance on the classification task does not necessarily mean better unsupervised model learning.

Example 4: Detection of Spatio-Temporal Spike Patterns

Our final application demonstrates that the modules for Bayesian computation that emerge in WTA circuits through STDP can not only explain the emergence of feature maps in primary sensory cortices like in Fig. 5, but could also be viewed as generic computational units in generic microcircuits throughout the cortex. Such generic microcircuit receives spike inputs from many sources, and it would provide a very useful computational operation on these if it could autonomously detect repeatedly occurring spatio-temporal patterns within this high-dimensional input stream, and report their occurrence through a self-organizing sparse coding scheme to other microcircuits. We have created such input streams with occasionally repeated embedded spike patterns for the computer experiment reported in Fig. 7. Fig. 7D demonstrates that sparse output codes for the 5 embedded spike patterns emerge after applying STDP in a WTA circuit for 200 s to such input stream. Furthermore, we show in the Supplement that these sparse output codes generalize (even without any further training) to time-warped versions of these spike patterns.

Even though our underlying probabilistic generative model (Fig. 1B) does not include time-dependent terms, the circuit in this example performs inference over time. The reason for this is that synapses that were active when a neuron fired become reinforced by STDP, and therefore make the neuron more likely to fire again when a similar spatial pattern is observed. Since we use EPSPs that smoothly decay over time, one neuron still sees a trace of previous input spikes as it fires again, and thus different spatial patterns within one reoccurring spatio-temporal pattern are recognized by the same neuron. The maximum length for such patterns is determined by the time constants of EPSPs. With our parameters (1 ms rise, 15 ms decay time constant) we were able to recognize spike patterns up to 50–100 ms. For longer spatio-temporal patterns, different neurons become responsive to different parts of the pattern. The neuron that responds mostly to noise in Figs. 7D did not learn a specific spatial pattern, and therefore wins by default when none of the specialized neurons responds. Similar effects have previously been described [59], [68], but for different neuron models, classical STDP curves, and not in the context of probabilistic inference.

For this kind of task, where also the exact timing of spikes in the patterns matters (which is not necessarily the case in the examples in Figs. 3, 5, and 6, where input neurons generate Poisson spike trains with different rates), we found that the shape of the STDP kernel plays a larger role. For example, a time-inverted version of the SEM rule, where pre-before-post firing causes LTD instead of LTP, cannot learn this kind of task, because once a neuron has learned to fire for a sub-pattern of the input, its firing onset is shifted back in time, rather than forward in time, which happens with standard SEM, but also with classical STDP [50], [59]. Instead, with a time-inverted SEM rule, different neurons would learn to fire stronger for the offsets of different patterns.

Such emergent compression of high-dimensional spike inputs into sparse low-dimensional spike outputs could be used to merge information from multiple sensory modalities, as well as from internal sources (memory, predictions, expectations, etc.), and to report the co-occurrence of salient events to multiple other brain areas. This operation would be useful from the computational perspective no matter in which cortical area it is carried out. Furthermore, the computational modules that we have analyzed can easily be connected to form networks of such modules, since their outputs are encoded in the same way as their inputs: through probabilistic spiking populations that encode for abstract multinomial variables. Hence the principles for the emergence of Bayesian computation in local microcircuits that we have exhibited could potentially also explain the self-organization of distributed computations in large networks of such microcircuits.

Prior related work

A first model for competitive Hebbian learning paradigm in non-spiking networks of neurons had been introduced in [77]. They analyzed a Hebbian learning rule in a hard WTA network and showed that there may exist equilibrium states, in which the average change of all weight values vanishes for a given set of input patterns. They showed that in these cases the weights adopt values that are proportional to the conditional probability of the presynaptic neuron being active given that the postsynaptic unit wins (rather than the log of this conditional probability, as in our framework). [78] showed that the use of a soft competition instead of a hard winner assignment and corresponding average weight updates lead to an exact gradient ascent on the log-likelihood function of a generative model of a mixture of Gaussians. However, these learning rules had not yet been analyzed in the context of EM.

Stochastic approximation algorithms for expectation maximization [31] were first considered in [79], incremental and on-line EM algorithms with soft-max competition in [80]–[82]. A proof of the stochastic approximation convergence for on-line EM in exponential family models with hidden variables was shown in [29]. They developed a sophisticated schedule for the learning rate in this much more general model, but did not yet consider individual learning rates for different weights.

[54] initiated the investigation of STDP in the context of unsupervised competitive Hebbian learning and demonstrated that correlations of input spike trains can be learned in this way. They also showed that this leads to a competition between the synapses for the control of the timing of the postsynaptic action potential. A similar competition can also be observed during learning in our model, since our learning rule automatically drives the weights towards satisfying the normalization conditions in Eq. (10).

[83] present a network and learning model that is designed to perform Independent Component Analysis (ICA) with spiking neurons through STDP and intrinsic plasticity. The mixture model of independent components can also be formulated as a generative model, and the goal of ICA is to find the optimal parameters of the mixing matrix. It has been shown that also this problem can be solved by a variant of Expectation Maximization [84], so there is some similarity to the identification of hidden causes in our model.

Recently, computer experiments in [85], [86] have used STDP in the context of WTA circuits to achieve a clustering of input patterns. Their STDP rules implements linear updates, independent of the current weight values, mixed with a homeostasis rule to keep the sum of all weights constant and every weight between 0 and 1. This leads to weights that are roughly proportional to the probability of the presynaptic neuron's firing given that the post-synaptic neuron fires afterwards. The competition between the output neurons is carried out as hard-max. In [85] the 4 output neurons learn to differentiate the 4 presented patterns and smoothly interpolate new rotated input patterns, whereas in [86] 48 neurons learn to differentiate characters in a small pixel raster. [86] uses a STDP rule where both LTP and LTD are modeled as exponentially dependent on the time difference. However, the very specific experimental setting with synchronous regular firing of the input neurons makes it difficult to generalize their result to more general input spike trains. No theoretical analysis is provided in [85] or [86], but their experimental results can be explained by our SEM approach. Instead of adding up logs of conditional probabilities and performing the competition on the exponential of the sums, they sum up the conditional probabilities directly and use this sum of probabilities for the competition. This can be seen as a linear approximation of SEM, especially under the additional normalization conditions that they impose by homeostasis rules.

It has previously been shown that spike patterns embedded in noise can be detected by STDP [49], [50], [59]. Competitive pattern learning through STDP has recently been studied in [68]. They simulate a deterministic version of a winner-take-all circuit consisting of a fixed number of neurons, all listening to the same spiking input lines and connected to each other with a strong inhibition. The STDP learning rule that they propose is additive and weight-independent. Just like our results, they also observe that different neurons specialize on different fixed repeated input pattern, even though the repeated patterns are embedded in spiking noise such that the mean activity of all inputs remains the same throughout the learning phase. Additionally they show that within each pattern the responsible neuron tries to detect the start of the pattern. In contrast to our approach they do not give any analysis of convergence guarantees, nor does their model try to build a generative probabilistic model of the input distribution.

[87]–[89] investigated the possibility to carry out Bayesian probabilistic computations in recurrent networks of spiking neurons, both using probabilistic population codes. They showed that the ongoing dynamics of belief propagation in temporal Bayesian models can be represented and inferred by such networks, but they do not exhibit any neuronal plausible learning mechanism. [7] presented another approach to Bayesian inference using probabilistic population codes, also without any learning result.

An interesting complementary approach is presented in [90], [91], where a single neuron is modeled as hidden Markov model with two possible states. This approach has the advantage, that the instantaneous synaptic input does not immediately decide the output state, but only incrementally influences the probability for switching the state. The weights and the temporal behavior can be learned online using local statistics. The downside of this approach is that this hidden Markov model can have only two states. In contrast, the SEM approach can be applied to networks with any number of output neurons.

In [92] it was shown that a suitable rule for supervised spike-based learning (the Tempotron learning rule) can be used to train a network to recognize spatio-temporal spike patterns. This discriminative learning scheme enables the recognizing neuron to focus on the most discriminative segment of the pattern. In contrast, our generative unsupervised learning scheme drives the recognizing neuron to generalize and spike many times during the whole pattern, and thus learns the spatial average activity pattern. The conductance based approach of [92] differs drastically from our method (and the results shown in the Supplement) insofar as here only STDP was used (focusing on average spatial patterns), no supervision was involved, and the time-warped input pattern had never been shown during training.

An alternative approach to implement the learning of generative probabilistic models in spiking neuronal networks is given in [93], [94]. Both approaches are based on the idea to model a sequence of spikes in a Hidden-Markov-Model-like probabilistic model and learn the model parameters through different variants of EM, in which a sequence of spikes represents one single sample of the model's distribution. Due to the explicit incorporation of inference over time, these models are more powerful than ours and thus require non-trivial, non-local learning mechanisms.

Experimentally testable predictions of the proposed model

Our analysis has shown, that STDP supports the creation of internal models and implements spike-based EM if changes of synaptic weights depend in a particular way on the current value of the weight: Weight potentiation depends in an inversely exponential manner on the current weight (see Eq. (5)). This rule for weight potentiation (see Fig. 8A) is consistent with all published data on this dependence: Fig. 5 in [43] and Fig. 5C in [40] for STDP, as well as Fig. 10 in [95] and Fig. 1 in [96] for other protocols for LTP induction. One needs to say, however, that these data exhibit a large trial-to-trial variability, so that it is hard to infer precise quantitative laws from them. On the other hand, the applications of STDP that we have examined in Fig. 3–7 work almost equally well if the actual weight increase varies by up to 100% from the weight increase proposed by our STDP rule (see open circles in Fig. 8A). The resulting distribution of weight increases matches qualitatively the above mentioned experimental data quite well.

The prediction of our model for the dependence of the amount of weight depression on the current weight is drastically different: Even though we make the strong simplification that the depression part of the STDP rule is independent of the time difference between pre- and postsynaptic spike, the formulation in Eq. (5) makes the assumption, that the amount of the depression should be independent of the current weight value. It is this contrast between an exponential dependency for LTP and a constant LTD which makes the weight converge to the logarithm of the conditional presynaptic firing probability in Eq. (6). In experiments this dependency has been investigated in-vitro [40]. There it has been found that the percentage of weight depression under STDP is independent of the current weight, which implies that the amount of depression is linear in the current weight value. This seems to contradict the presented learning rule. However, the key property that is needed for the desired equilibrium condition is the ratio between LTP and LTD. So the equilibrium proof in Eq. (28) remains unchanged if 👁 Image
is multiplied (for potentiation and depression) by some arbitrary function 👁 Image
of the current weight value. Choosing for example 👁 Image
yields a depression whose percentage is independent of the initial value, which would be consistent with the above mentioned in-vitro data [40]. The resulting dependence for potentiation is plotted in Fig. 8B. Since this curve is very similar to that of Fig. 8A, the above mentioned experimental data for potentiation are too noisy to provide a clear vote for one of these two curves. Thus more experimental data are needed for determining the dependence of weight potentiation on the initial weight. Whereas the relevance of this dependency had previously not been noted, our analysis suggests that such a contrast it is in fact essential for the capability of STDP to create internal models for high-dimensional spike inputs.

Our analysis has shown, that if the excitability of neurons is also adaptive, with a rule as in Eq. (7) that is somewhat analogous to that for synaptic plasticity, then neurons can also learn appropriate priors for Bayesian computation. Several experimental studies have already confirmed, that the intrinsic excitability of neurons does in fact increase when they are more frequently activated [45], see [97], [14] and [39] for reviews. But a quantitative study, which relates the resulting change in intrinsic excitability to its initial value, is missing.

Our model proposes that pyramidal neurons in cortical microcircuits are organized into stochastic WTA circuits, that together represent a probability distribution. This organization is achieved by a suitably regulated common inhibitory signal, where the inhibition follows the excitation very closely. Such instantaneous balance between excitation and inhibition was described by [34]. A resulting prediction of the WTA structure is that the firing activity of these neurons is highly de-correlated due to the inhibitory competition. In contrast to previous experimental results, that reported higher correlations, it has recently been confirmed in [35] for the visual cortex of awake monkey that nearby neurons, even though they share common input show extremely low correlations.

Another prediction is that neural firing activity especially for awake animals subject to natural stimuli is quite sparse, since only those neurons fire whose internal model matches their spike input. A number of experimental studies confirm this predictions (see [98] for a review). Our model also predicts, that the neural firing response to stimuli exhibits a fairly high trial-to-trial variability, as is typical for drawing repeated samples from a posterior distribution (unless the posterior probability is close to 0 or 1). A fairly high trial-to-trial variability is a common feature of most recordings of neuronal responses (see e.g. [99], Fig. 1B in [100]; a review is provided in [101]). In addition, our model predicts that this trial-to-trial variability decreases for repeatedly occurring natural stimuli (especially if this occurs during attention) and discrimination capability improves for these stimuli, since the internal models of neurons are becoming better fitted to their spike input during these repetitions (“sharpening of tuning”), yielding posterior probabilities closer to 1 or 0 for these stimuli. These predictions are consistent with a number of experimental data related to perceptual learning [102], [103], and with the evolution of neuronal responses to natural scenes that were shown repeatedly in conjunction with nucleus basalis stimulation [104].

In addition our model predicts that if the distribution of sensory inputs changes, the organization of codes for such sensory inputs also changes. More frequently occurring sensory stimuli will be encoded with a finer resolution (see [105] for a review of related experimental data). Furthermore in the case of sensory deprivation (see [106]) our model predicts that neurons that used to encode stimuli which no longer occur will start to participate in the encoding of other stimuli.

We have shown in Fig. 3 that an underlying background oscillation on neurons that provide input to a WTA circuit speeds up the learning process, and produces more precise responses after learning. This result predicts that cortical areas that collaborate on a common computational task, especially under attention, exhibit some coherence in their LFP. This has already been shown for neurons in close proximity [107] but also for neurons in different cortical areas [108], [109].

If one views the modules for Bayesian computation that we have analyzed in this article as building blocks for larger cortical networks, these networks exhibit a fundamental difference to networks of neurons: Whereas a neuron needs a sufficiently strong excitatory drive in order to reach its firing threshold, the output neurons 👁 Image
of a stochastic WTA circuit according to our model in Eq. (3) are firing already on their own - even without any excitatory drive from the input neuron 👁 Image
(due to assumed background synaptic inputs; modeled in our simulations by an Ornstein-Uhlenbeck process, as suggested by in-vivo data [110]). Rather, the role of the input from the 👁 Image
-neurons is to modulate which of the neurons in the WTA circuit fire. One consequence of this characteristic feature is that even relatively few presynaptic neurons 👁 Image
can have a strong impact on the firing of the 👁 Image
-neurons, provided the 👁 Image
-neurons have learned (via STDP) that these 👁 Image
-neurons provide salient information about the hidden cause for the total input 👁 Image
from all presynaptic neurons. This consequence is consistent with the surprisingly weak input from the LGN to area V1 [16], [111], [112]. It is also consistent with the recently found exponential distance rule for the connection strength between cortical areas [112]. This rule implies that the connection strength between distal cortical areas, say between primary visual cortex and PFC, is surprisingly weak. Our model suggests that these weak connections can nevertheless support coherent brain computation and memory traces that are spread out over many, also distal, cortical areas.

Apart from these predictions regarding aspects of brain computation on the microscale and macroscale, a primary prediction of our model is that complex computations in cortical networks of neurons - including very efficient and near optimal processing of uncertain information - are established and maintained through STDP, on the basis of genetically encoded stereotypical connection patterns (WTA circuits) in cortical microcircuits.

Equilibrium condition

We will now show that all equilibria of the stochastic update rule in Eq. (5) and Eq. (7), i.e., all points where 👁 Image
, exactly match the implicit solution conditions in Eq. (46), and vice versa:👁 Image
(28)Analogously, one can show that 👁 Image
. Note that this result implies that the learning rule in Eq. (5) and Eq. (7) has no equilibrium points outside the normalization conditions in Eq. (10), since all equilibrium points fulfill the implicit solutions condition in Eq. (46) and these in turn fulfill the normalization conditions.

Details to Learning the parameters of the probability model by EM

In this section we will analyze the theoretical basis for learning the parameters 👁 Image
of the generative probability model 👁 Image
given in Eq. (9) from a machine learning perspective. In contrast to the intuitive explanation of the Results section which was based on Expectation Maximization we will now derive an implicit analytical solution for a (locally) optimal weight vector 👁 Image
, and rewrite this solution in terms of log probabilities. We will later use this derivation in order to show that the stochastic online learning rule provably converges towards this solution.

For an exact definition of the learning problem, we assume that the input is given by a stream of vectors 👁 Image
, in which every 👁 Image
is drawn independently from the input distribution 👁 Image
. In principle, this stream of 👁 Image
's corresponds to the samples 👁 Image
that are observed at the spike times 👁 Image
of the circuit. However, in order to simplify the proofs in this and subsequent sections, we will neglect any possible temporal correlation between successive samples.

The learning task is to find parameter values 👁 Image
, such that the marginal 👁 Image
of the model distribution 👁 Image
approximates the actual input distribution 👁 Image
as accurately as possible. This is equivalent to minimizing the Kullback-Leibler divergence between the two distributions:👁 Image
(29)where 👁 Image
is the (constant) entropy of the input distribution 👁 Image
, and 👁 Image
denotes the expectation over 👁 Image
, according to the distribution 👁 Image
. Since 👁 Image
is constant, minimizing the right hand side of Eq. (29) is equivalent to maximizing the expected log likelihood 👁 Image
.

There are many different parametrizations 👁 Image
that define identical generative distributions 👁 Image
in Eq. (24). There is, however, exactly one 👁 Image
in this sub-manifold of the weight space that fulfills the normalization conditions in Eq. (10).

We thus redefine the goal of learning more precisely as the constrained maximization problem👁 Image
(30)👁 Image
(31)

This maximization problem never has a unique solution 👁 Image
, because any permutation of the values of 👁 Image
and their corresponding weights leads to different joint distributions 👁 Image
, all of them having identical marginals 👁 Image
. The local maxima of Eq. (30) can be found using the Lagrange multiplier method.

Note that we do at no time enforce normalization of 👁 Image
during the learning process, nor do we require normalized initialization of 👁 Image
. Instead, we will show that the learning rule in Eq. (5,7) automatically drives 👁 Image
towards a local maximum, in which the normalization conditions are fulfilled.

Under the constraints in Eq. (31) the normalization constant 👁 Image
in Eq. (21) equals 👁 Image
, thus 👁 Image
simplifies to 👁 Image
- with 👁 Image
- and we can define a Lagrangian function 👁 Image
for the maximization problem in Eq. (30,31) by👁 Image
(32)Setting the derivatives to zero we arrive at the following set of equations in 👁 Image
and 👁 Image
:👁 Image
(33)👁 Image
(34)Summing over those equations that have the same multiplier 👁 Image
or 👁 Image
, resp., leads to👁 Image
(35)👁 Image
(36)where 👁 Image
is the shorthand notation for the equivalent expression 👁 Image
. The identity 👁 Image
, the identity 👁 Image
the fact that 👁 Image
, which follows from the definition of population encoding, and the constraints in Eq. (31) are used in order to derive the explicit solution for the Lagrange multipliers👁 Image
(37)in dependence of 👁 Image
. We insert this solution for 👁 Image
into the gradient Eq. (33,34) and get👁 Image
(38)👁 Image
from which we derive an implicit solution for 👁 Image
:👁 Image
(39)👁 Image
It is easily verified that all fixed points of this implicit solution satisfy the normalization constraints:👁 Image
(40)👁 Image
(41)

Finally, in order to simplify the notation we use the augmented input distribution 👁 Image
. The expectations in Eq. (39) nicely evaluate to👁 Image
(42)👁 Image
(43)👁 Image
(44)👁 Image
(45)which allows us to rewrite the implicit solution in a very intuitive form as:👁 Image
(46)Any weight vector 👁 Image
that fulfills Eq. (46) is either a (local) maximum, a saddle point or a (local) minimum of the log likelihood function 👁 Image
under the normalization constraints.

An obvious numerical approach to solve this fixed point equation is the repeated application of Eq. (39). According to the derivations in the Results section this corresponds exactly to the Expectation Maximization algorithm. But every single iteration asks for the evaluation of expectations with respect to the input distribution 👁 Image
, which theoretically requires infinite time in an online learning setup.

Details to Spike-based Expectation Maximization

We derive the update rule in Eq. (5) from the statistical perspective that each weight can be interpreted as 👁 Image
, where 👁 Image
and 👁 Image
correspond to counters of the events 👁 Image
and 👁 Image
. Every new event 👁 Image
leads to a weight update👁 Image
(47)👁 Image
(48)👁 Image
(49)👁 Image
(50)where the log-function is linearly approximated around 1 as 👁 Image
. The factor 👁 Image
is understood as learning rate 👁 Image
in the additive update rule 👁 Image
. If 👁 Image
, i.e. if there is no postsynaptic spike, the update 👁 Image
. In the case of a postsynaptic spike, i.e., 👁 Image
, the update 👁 Image
decomposes in the two cases 👁 Image
and 👁 Image
as it is stated explicit in Eq. (5).

As a side note, we observe that by viewing our STDP rule as an approximation to counting statistics, the learning rate 👁 Image
can be understood as the inverse of the equivalent sample size from which the statistics was gathered. If the above rule is used with a small constant learning rate we will get a close approximation to an exponentially decaying average. If the learning rate decays like 👁 Image
we will get an approximation to an online updated average, where all samples are equally weighted. We will come back to a regulation mechanism for the learning rate in the section ‘Variance Tracking’.

Details to Proof of convergence

In this section we give the proof of Theorem 1. Formally, we define the sequences 👁 Image
, 👁 Image
, 👁 Image
, 👁 Image
and 👁 Image
for 👁 Image
: For all 👁 Image
we assume that 👁 Image
is drawn independently from 👁 Image
. The value of 👁 Image
is drawn from the posterior distribution of the model 👁 Image
(see Eq. (11)), given the input 👁 Image
and the current model parameters 👁 Image
. The weight updates 👁 Image
, and 👁 Image
, are calculated according to Eq. (5) and (7) with 👁 Image
. The sequence of weight vectors 👁 Image
is determined by the randomly initialized vector 👁 Image
, and by the iteration equation👁 Image
(51)The projection function 👁 Image
represents a coordinate-wise clipping of 👁 Image
to a hyper-rectangle 👁 Image
such that 👁 Image
(52)The bound 👁 Image
is assumed to be chosen so that all (finite) maxima of 👁 Image
are inside of 👁 Image
. For the sequence of learning rates 👁 Image
we assume that👁 Image
(53)

Under these assumptions we can now restate the theorem formally:

Theorem 1: The sequence 👁 Image
converges with probability 1 to the set 👁 Image
of all points within the hyper-rectangle 👁 Image
that fulfill the equilibrium conditions in Eq. (6). The stable convergence points among 👁 Image
are the (local) maxima of 👁 Image
, subject to the normalization constraints in Eq. (10).

The iterative application of the learning rule in Eq. (5) and (7) is indeed a stochastic approximation algorithm for learning a (locally) optimal parameter vector 👁 Image
. We resort to the theory of stochastic approximation algorithms as presented in [52] and use the method of the “mean limit” ordinary differential equation (ODE). The goal is to show that the sequence of the weight vector 👁 Image
under the stochastic learning rule in Eq. (5) and (7) converges to one of the local maxima of Eq. (30) with probability one, i.e., the probability to observe a non-converging realization of this sequence is zero. The location of the local maximum to which a single sequence of 👁 Image
converges depends on the starting point 👁 Image
as well as on the concrete realization of the stochastic noise sequence. We will not discuss the effect of this stochasticity in more detail, except for stating that a stochastic approximation algorithm is usually less prone to get stuck in small local maxima than its deterministic version. The stochastic noise introduces perturbations that decrease slowly over time, which has an effect that is comparable to simulated annealing.

We will use the basic convergence theorem of [52] to establish the convergence of the sequence 👁 Image
to the limit set of the mean limit ODE. Then it remains to show that this limit set is identical to the desired set of all equilibrium points and thus, particularly, does not contain limit cycles.

Proof: In the notation of [52], the mean update of the stochastic algorithm in Eq. (51) is 👁 Image
. The bounds 👁 Image
imply that 👁 Image
for all 👁 Image
and 👁 Image
.

For any set 👁 Image
we define 👁 Image
as the positive limit set of the mean limit ODE 👁 Image
for all initial conditions 👁 Image
:👁 Image
(54)

According to Theorem 3.1 in Chapter 5 of [52], the sequence 👁 Image
under the algorithm in Eq. (51) converges for all start conditions 👁 Image
to the limit set 👁 Image
with probability one in the sense that👁 Image
(55)

We will now show that the limit set 👁 Image
of 👁 Image
is identical to the set of stationary points 👁 Image
and does not contain limit cycles. It is obvious that 👁 Image
is a subset of 👁 Image
since for all initial conditions 👁 Image
the trajectory of 👁 Image
fulfills 👁 Image
for all 👁 Image
. Thus it remains to be shown that there are no other points in 👁 Image
(like e.g. limit cycles).

We split the argument into two parts. In the first part we will show that for 👁 Image
all trajectories of 👁 Image
converge asymptotically to the manifold 👁 Image
defined by the normalization constraints 31. This leads to the conclusion that 👁 Image
. In the second part we will show that all trajectories within 👁 Image
converge to the stationary points 👁 Image
, i.e., 👁 Image
. Both parts together yield the desired result that 👁 Image
are the only limit points of the ODE 👁 Image
.

The first part we start by defining the set of functions 👁 Image
and 👁 Image
for all 👁 Image
to represent the deviation of the current 👁 Image
from each of the normalization constraints 31, i.e.,👁 Image
(56)The manifold 👁 Image
is the set of all points 👁 Image
where 👁 Image
and 👁 Image
for all 👁 Image
. Furthermore, we calculate the gradient vectors 👁 Image
and 👁 Image
for each of these functions with respect to the argument 👁 Image
. Note that many entries of these gradient vectors are 👁 Image
, since every single function 👁 Image
and 👁 Image
only depends on a few entries of its argument 👁 Image
. The nonzero entries of these gradients are👁 Image
(57)We can now show that the trajectory of 👁 Image
in any point 👁 Image
always points in direction of decreasing absolute values for all deviations 👁 Image
and 👁 Image
:👁 Image
(58)👁 Image
(59)👁 Image
(60)👁 Image
(61)This shows that 👁 Image
for all 👁 Image
and 👁 Image
. This implies that the limit set of all trajectories with initial conditions outside 👁 Image
is contained in 👁 Image
, or more formally 👁 Image
. Note that the continuity and the boundedness of 👁 Image
on 👁 Image
implies 👁 Image
and 👁 Image
if 👁 Image
for all 👁 Image
. Therefore we can now conclude as the result of the first part👁 Image
(62)i.e. the limit set of all trajectories starting outside the manifold of normalized weights is contained in the limit set of all trajectories starting within the normalization constraints. The equations (61) also prove that any trajectory with initial condition 👁 Image
stays within 👁 Image
, since all components of 👁 Image
with directions orthogonal to the tangent space of 👁 Image
in 👁 Image
are 👁 Image
for all 👁 Image
, thus 👁 Image
is in the tangent space 👁 Image
in 👁 Image
.

This immediately leads to the second part of the proof, which is based on the gradient 👁 Image
of the Lagrangian 👁 Image
as given in Eq. (33, 34). For any 👁 Image
let 👁 Image
be the linear projection matrix that orthogonally projects any vector 👁 Image
into the tangent space of 👁 Image
in 👁 Image
. The projection 👁 Image
of the gradient of 👁 Image
at any 👁 Image
points towards the strongest increase of the value of the objective function 👁 Image
under the constraints of the normalization conditions. Thus, the value of 👁 Image
increases in the direction of any vector within the tangent space of 👁 Image
in 👁 Image
that has a positive scalar product with 👁 Image
. As 👁 Image
is a tangent vector of 👁 Image
in 👁 Image
for all 👁 Image
, the orthogonal component 👁 Image
of the gradient is orthogonal to 👁 Image
. Thus, the value of the scalar product with the projected gradient 👁 Image
is identical to the value of the scalar product with the gradient itself 👁 Image
:👁 Image
(63)👁 Image
with equality if and only if 👁 Image
, which is equivalent to 👁 Image
. This shows that all trajectories with initial condition 👁 Image
stay within 👁 Image
forever and converge to the set of stationary points 👁 Image
, i.e. 👁 Image
. Combining the results of both parts as👁 Image
(64)establishes the stochastic convergences of any sequence 👁 Image
to the set 👁 Image
with probability one.

Details to Role of the Inhibition

Details to Continuous-Time Interpretation with Realistically Shaped EPSPs

Let the external input vector 👁 Image
consist of multiple discrete-valued functions in time 👁 Image
, and let us assume that for every input 👁 Image
there exists an independent Poisson sampling process with rate 👁 Image
which generates spike times for the group of neurons 👁 Image
with 👁 Image
. At every spike time 👁 Image
there is exactly one neuron in the group that fires a spike, and this is the neuron that is associated with the value 👁 Image
. First, we analyze additive step-function EPSPs, i.e. the postsynaptic activation 👁 Image
is given by the convolution in Eq. (17) where 👁 Image
is a step-function kernel with 👁 Image
for 👁 Image
for a fixed EPSP-duration 👁 Image
and 👁 Image
otherwise. In order to understand the resulting distribution 👁 Image
in Eq. (4) as Bayesian inference we extend our underlying generative probabilistic model 👁 Image
such that it contains multiple instances of the variable vector 👁 Image
, called 👁 Image
, where 👁 Image
is the total number of spikes from all input neurons 👁 Image
within the time window 👁 Image
. We can see every spike as a single event in continuous time. The full probabilistic model is defined as👁 Image
(80)which defines that the multiple instances are modeled as being conditionally independent of each other, given 👁 Image
. Let the vectors 👁 Image
describe the corresponding spike “patterns” in which every binary vector 👁 Image
has exactly one 👁 Image
entry 👁 Image
. All other values are zero, thus it represents exactly one evidence for 👁 Image
, i.e. 👁 Image
, with 👁 Image
, s.t. 👁 Image
, according to the decoding in Eq. (8).

Due to the conditional independences in the probabilistic model every such evidence, i.e. every spike, contributes one factor 👁 Image
to the likelihood term in the inference of the hidden node 👁 Image
. The inference is expressed as 👁 Image
(81)The identity 👁 Image
reveals that the above posterior distribution is realized by the relative spike probability 👁 Image
of the network model according to Eq. (4), where 👁 Image
replaces 👁 Image
in the computation of the membrane potential 👁 Image
. Due to the step function 👁 Image
the result of the convolution in 👁 Image
equals the number of spikes within the time window 👁 Image
from neuron 👁 Image
. The factor 👁 Image
, which has the meaning 👁 Image
in the network model, is multiplied 👁 Image
times to the likelihood.

The above discrete probabilistic model gives an interpretation only for integer values of 👁 Image
, i.e. for functions 👁 Image
such that 👁 Image
is 👁 Image
or any positive integer at any time 👁 Image
. For an interpretation of arbitrarily shaped EPSPs 👁 Image
- especially for continuously decaying functions - in the context of our probabilistic model, we now extend this weighting mechanism from integer valued weights to real valued weights by a linear interpolation of the likelihood in the log-space.

The obvious restrictions on the EPSP function 👁 Image
are that it is non-negative, zero for 👁 Image
, and 👁 Image
, in order to avoid acausal or nondecaying behavior, and unboundedly growing postsynaptic potentials at constant input rates. We assume the normalization 👁 Image
. Let again 👁 Image
be the times of the past spiking events and 👁 Image
be the indices of the corresponding input neurons. The output distribution 👁 Image
can be written as👁 Image
(82)which nicely illustrates that every single past spike at time 👁 Image
is seen as an evidence in the inference, but that evidence is weighted with a value 👁 Image
, which is between 0 and 1.

The analogous interpolation for continuous-valued input activations 👁 Image
yields the learning rule in Eq. (18), which is illustrated in Fig. 2 as the “Complex STDP rule” (blue dashed curve). The resulting shape of the LTP part of the STDP curve is determined by the EPSP shape defined by 👁 Image
. The positive part of the update in Eq. (18) is weighted by the value of 👁 Image
at the time of firing the postsynaptic spike. Negative updates are performed if 👁 Image
is close to zero, which indicates that no presynaptic spikes were observed recently.

The proof of stochastic convergence does not explicitly assume that 👁 Image
is a binary vector, but is valid for any (positive) random variable vector 👁 Image
with finite variance. Further, the proof assumes the condition that in every group 👁 Image
the sum of the input activities 👁 Image
is 👁 Image
at all times or at least at those points in time at which one 👁 Image
neuron of the WTA-circuit fires. The condition can be relaxed such that the sum per group does not have to be equal to 👁 Image
but to any arbitrary (positive) constant if the corresponding normalization constraint is adapted accordingly. Due to the decaying character of the EPSP shape, this sum will never stay constant, even for very regular input patterns. If we only assumed a constant average activation within a group, allowing for stochastic fluctuations around the target value, it turns out that this condition alone is not enough. We need to further assume that these stochastic fluctuations in the sum of every input group 👁 Image
are stochastically independent of the circuit's response 👁 Image
. This assumption is intricate and may depend on the data and the learning progress itself, so it will usually not be exactly fulfilled. We can, however, argue that we are close to independence if at least the sum of activity in every group 👁 Image
is independent of the value of the underlying abstract variable 👁 Image
.

In our simulations we obtain the input activations 👁 Image
by simulating biologically realistic EPSPs at every synapse, using 👁 Image
-kernels with plausible time constants to model the contributions of single input spikes.

Details to Spike-timing dependent LTD

We formalize the presynaptic activity of neuron 👁 Image
after a postsynaptic spike at time 👁 Image
by 👁 Image
, s.t. 👁 Image
if there is a spike from neuron 👁 Image
within the time window 👁 Image
and 👁 Image
otherwise. This trace is used purely for mathematical analysis, and cannot be known to the postsynaptic neuron at time 👁 Image
, since the future input activity is unknown. Mechanistically, however, 👁 Image
can be implemented as a trace updated by postsynaptic firing, and utilized for plasticity at the time of presynaptic firing [113]. Let us now consider the STDP rule illustrated by the red curve in Fig. 9, where a depression of the synapse happens only if there is a presynaptic spike within the short time window of length 👁 Image
after the postsynaptic spike, i.e. if 👁 Image
. The application of this STDP-rule in our neuronal circuit is equivalent to the circuit-spike triggered update rule👁 Image
(83)which replaces Eq. (5). In analogy to Eq. (6) the equilibrium of this new update rule can be derived as👁 Image
(84)👁 Image
(85)under the assumption that 👁 Image
and 👁 Image
are sampled from a stationary distribution 👁 Image
. This shows that the synaptic weights can be interpreted as the log-likelihood ratio of the presynaptic neuron firing before instead of after the postsynaptic neuron. In other words, the neuron's synaptic weights learn the contrast between the current input pattern 👁 Image
that caused firing, and the following pattern of activity 👁 Image
. Note that any factor 👁 Image
(for LTP) or 👁 Image
(for LTD) only leads to a constant offset of the weight which - under the assumption that the offset is the same for all synapses - can be neglected due to the WTA circuit (see Methods “Weight offsets and positive weights”).

Similarly to our analysis for the standard SEM rule, we can derive a continuous-time interpretation of the timing-dependent LTD rule. As we did in Eq. (17), we can define👁 Image
(86)where 👁 Image
is the same convolution kernel as in Eq. (17), and 👁 Image
is an arbitrary but time-inversed kernel, such that 👁 Image
for positive 👁 Image
and 👁 Image
for negative 👁 Image
. The value of 👁 Image
thus reflects a time-discounted sum of presynaptic activity immediately after the postsynaptic spike.

The complex STDP rule from Fig. 2, which models LTD as a constant time-independent depression, can be seen as an extreme case of the spike-timing dependent LTD rule. If 👁 Image
is a step function with 👁 Image
in the interval 👁 Image
and 👁 Image
everywhere else, then 👁 Image
is just the average rate of presynaptic activity in the time interval 👁 Image
following a postsynaptic spike. In the limit of 👁 Image
this is equivalent to the overall spiking rate of the neuron 👁 Image
, which is proportional to the marginal 👁 Image
in the probabilistic model. Precisely, 👁 Image
, where 👁 Image
is the base firing rate of an active input in our input encoding model. The equilibrium point of every weight 👁 Image
becomes 👁 Image
, neglecting the offsets induced by the constants 👁 Image
,👁 Image
and 👁 Image
. It is easy to see that the probabilistic interpretation of the neuronal model from Eq. (4) is invariant under the transformation 👁 Image
, since👁 Image
(87)👁 Image
(88)👁 Image
(89)which proves that in our network model the complex STDP rule from Fig. 2 is equivalent to an offset-free STDP rule in the limit of an arbitrarily long window for LTD. In practice, of course, we can assume that the times between pre- and post-synaptic spikes are finite, and we have shown in Fig. 4 that as a result, very realistic shapes of STDP curves emerge at intermediate stimulation frequencies.