Optimal Localist and Distributed Coding of Spatiotemporal Spike Patterns Through STDP and Coincidence Detection

Repeating spatiotemporal spike patterns exist and carry information. Here we investigated how a single spiking neuron can optimally respond to one given pattern (localist coding), or to either one of several patterns (distributed coding, i.e., the neuron's response is ambiguous but the identity of the pattern could be inferred from the response of multiple neurons), but not to random inputs. To do so, we extended a theory developed in a previous paper (Masquelier, 2017), which was limited to localist coding. More specifically, we computed analytically the signal-to-noise ratio (SNR) of a multi-pattern-detector neuron, using a threshold-free leaky integrate-and-fire (LIF) neuron model with non-plastic unitary synapses and homogeneous Poisson inputs. Surprisingly, when increasing the number of patterns, the SNR decreases slowly, and remains acceptable for several tens of independent patterns. In addition, we investigated whether spike-timing-dependent plasticity (STDP) could enable a neuron to reach the theoretical optimal SNR. To this aim, we simulated a LIF equipped with STDP, and repeatedly exposed it to multiple input spike patterns, embedded in equally dense Poisson spike trains. The LIF progressively became selective to every repeating pattern with no supervision, and stopped discharging during the Poisson spike trains. Furthermore, tuning certain STDP parameters, the resulting pattern detectors were optimal. Tens of independent patterns could be learned by a single neuron using a low adaptive threshold, in contrast with previous studies, in which higher thresholds led to localist coding only. Taken together these results suggest that coincidence detection and STDP are powerful mechanisms, fully compatible with distributed coding. Yet we acknowledge that our theory is limited to single neurons, and thus also applies to feed-forward networks, but not to recurrent ones.


Introduction
In a neural network, either biological or artificial, two forms of coding can be used: localist or distributed. With localist coding, each neuron codes (i.e. maximally responds) for one and only one category of stimulus (or stimulus feature). As a result, the category of the stimulus (or the presence of a certain feature) can be inferred from the response of this sole neuron, ignoring the other neurons' responses. Conversely, with distributed coding each neuron responds to multiple stimulus categories (or features) in a similar way. Therefore the response of each neuron is ambiguous, and the category of the stimulus, or the presence of a certain feature, can only be inferred from the responses of multiple neurons. Thus the distinction between the two schemes is the number of different stimuli to which a given neuron responds -not the number of neurons which respond to a given stimulus. Indeed, a localist network can have redundancy, and use multiple "copies" of each category specific neuron [Bowers, 2009, Thorpe, 1989.
Does the brain use localist or distributed coding? This question has been, and still is, intensively 1 This manuscript is now published in Frontiers in Computational Neuroscience. Please cite it as: Masquelier, T. and Kheradpisheh, S. R. (2018). Optimal localist and distributed coding of spatiotemporal spike patterns through STDP and coincidence detection. debated. In practice, discriminating between the two schemes from electrophysiological recordings is tricky [Quian Quiroga and Kreiman, 2010], since the set of tested stimuli is always limited, the responses are noisy, the thresholds are arbitrary and the boundaries between categories are fuzzy. Here we do not attempt to do a complete review of the experimental literature; but rather to summarize it. It is commonly believed that distributed coding is prevalent [Hung et al., 2005, O'Reilly, 1998, Quiroga et al., 2008, Rolls et al., 1997, but there is also evidence for localist coding, at least for familiar stimuli, reviewed in [Bowers, 2009, Roy, 2017, Thorpe, 2009. The question of localist vs. distributed coding is also relevant for artificial neural networks, and in particular for the recently popular deep neural networks. Most of the time, these networks are trained in a supervised manner, using the backpropagation algorithm [LeCun et al., 2015]. The last layer contains exactly one neuron per category, and backpropagation forces each neuron to respond more strongly when the stimulus belongs to the neuron's category. In other words, localist coding is imposed in the last layer. Conversely, the hidden layers are free to choose their coding scheme, which is supposedly optimal for the categorization task at hand. It is thus very interesting to analyze the chosen coding scheme. It is not easy to do such analysis on the brain (as explained above), but we can do it rigorously for computational models by computing the responses to huge amounts of images, and even synthesizing images that maximize the responses. Results indicate that some hidden neurons respond to one object category only [Nguyen et al., 2016, Olah et al., 2017, Zhou et al., 2015, while others respond to multiple different objects [Nguyen et al., 2016, Olah et al., 2017. Thus it appears that both localist and distributed codes can be optimal, depending on the task, the layer number, and the network parameters (number of layers, neurons, etc.).
Let us come back to the brain, in which computation is presumably implemented by spiking neurons performing coincidence detection [Abeles, 1982, Brette, 2015, König et al., 1996. This observation raises an important question, which we tried to address in this theoretical paper: can coincidence detector neurons implement both localist and distributed codes? In this context, different stimuli correspond to different spatiotemporal in-put spike patterns. Here each pattern was generated randomly, leading to chance-level overlap between patterns. In addition, each pattern was jittered at each presentation, resulting in categories of similar, yet different, stimuli. Can a neuron respond to one, or several of these patterns, and not to random inputs? What is the required connectivity to do so in an optimal way? And finally, can this required connectivity emerge with spiketiming-dependent plasticity (STDP), in an unsupervised manner?
To address these questions, we extended a theory that we developed in a previous paper, but which was limited to one pattern only, i.e. localist coding [Masquelier, 2017], to the multi-pattern case. Briefly, we derived analytically the signal-to-noise ratio (SNR) of a multi-pattern detector, and investigated the conditions for its optimality. In addition, using numerical simulations, we showed that a single neuron equipped with STDP can become selective to multiple repeating spike patterns, even without supervision and that the resulting detectors can be close to the theoretical optimum. Surprisingly, a single neuron could robustly learn up to ∼ 40 independent patterns (using parameters arguably in the biological range). This was not clear from previous simulations studies, in which neurons equipped with STDP only learned one pattern (localist coding) [Gilson et al., 2011, Humble et al., 2012, Hunzinger et al., 2012, Kasabov et al., 2013, Klampfl and Maass, 2013, Krunglevicius, 2015, Masquelier, 2017, Masquelier et al., 2008, Nessler et al., 2013, Sun et al., 2016, or two patterns [Yger et al., 2015]. This shows that STDP and coincidence detection are compatible with distributed coding.

Formal description of the problem
The problem we addressed is similar to the one of [Masquelier, 2017], but extended to the multipattern case. For the reader's convenience, we fully describe it below.
We addressed the problem of detecting one or several spatiotemporal spike patterns with a single LIF neuron. Intuitively, one should connect the neurons that are active during the patterns (or during subsections of them) to the LIF neuron. That way, the LIF will tend to be more activated by The LIF is connected to the neurons that fire in some subsections of the patterns with duration ∆t ≤ L (these emit red spikes) (Bottom) The LIF potential peaks for patterns, and the double arrow indicates the peak height.
the patterns than by some other inputs. More formally, we note P the number of spike patterns, and assume that they all have the same duration L. We note N the number of neurons involved. For each pattern, we chose a subsection with duration ∆t ≤ L, and we connect the LIF to the M neurons that emit at least one spike during at least one of these subsections (Fig. 1).
We hypothesize that all afferent neurons fire according to a homogeneous Poisson process with rate f , both inside and outside the patterns. That is the patterns correspond to some realizations of the Poisson process, which can be repeated (this is sometimes referred to a "frozen noise"). At each repetition a random time lag (jitter) is added to each spike, drawn from a uniform distribution over [−T, T ] (a normal distribution is more often used, but it would not allow analytical treatment [Masquelier, 2017]).
We also assume that synapses are instantaneous, which facilitates the analytic calculations.
For now we ignore the LIF threshold, and we want to optimize its signal-to-noise ratio (SNR), defined as: where V max is the maximal potential reached during the pattern presentations, V noise is the mean value for the potential with Poisson input (noise period), and σ noise is its standard deviation. Obviously, a higher SN R means a larger difference between the LIF membrane potential during the noise periods and its maximum value, which occurs during the selected ∆t window of each pattern. Therefore, the higher the SN R the lower the probability of missing patterns, and of false alarms. We consider that P , L, N , f and T are imposed variables, and that we have the freedom to choose ∆t ≤ L and the membrane time constant τ in order to maximize the SN R.
We note that this problem is related to the synfire chain theory [Abeles, 1991]. A synfire chain consists of a series of pools of neurons linked together in a feed-forward chain, so that volleys of synchronous spikes can propagate from pool to pool in the chain. Each neuron can participate in several of such chains. The number of different chains that can coexist in a network of a given size has been termed capacity. This capacity can be optimized [Herrmann et al., 1995]. To do so, a given neuron should respond to certain spike volleys, but not to others, which is similar to our optimization of a multi-pattern SN R. Yet it is also different: we use homogeneous Poisson activity, not spike volleys, and we ignore the threshold, while synfire chains require thresholds.

Deriving the SNR analytically
Here we are to find the optimum SN R of the LIF for P patterns. To this end we should first calculate the SN R analytically. Again, the derivations are similar to the ones in [Masquelier, 2017], but extended to the multi-pattern case (which turned to mainly impact Equation 7).
In this section, we assume non-plastic unitary synaptic weights. That is an afferent can be either connected (w = 1) or disconnected (w = 0) (in the Appendix we estimate the cost of this constraint 3 A theoretical optimum 4 on the SN R). Thus the LIF obeys the following differential equation: where the t i are the presynaptic spike times of all the connected afferents.
Since synapses are instantaneous and firing is Poissonian, during the noise periods and outside the ∆t windows we have: V noise = τ f M and σ noise = τ f M/2 [Burkitt, 2006], where M is the number of connected input neurons (with unitary weights).
To compute V max , it is convenient to introduce the reduced variable: where V ∞ = τ r is the mean potential of the steady regime that would be reached if ∆t was infinite, and r is the input spike rate during the ∆t window, resulting from the total received spikes from all input neurons during this window.
v max can be calculated by exact integration of the LIF differential equation [Masquelier, 2017]. Here we omit the derivation and present the final equation: Using the definition of v max in Equation 3, we can rewrite the SN R equation as: Obviously, different Poisson pattern realizations will lead to different values for M and r that consequently affect each of the terms V ∞ , V noise and σ noise . Here we want to compute the expected value of the SN R across different Poisson pattern real-izations: In Section 3.2 we justify this last approximation through numerical simulations, and we also show that this average SN R is not much different from the SN R of particular Poisson realizations. The last step to compute SN R in Equation 6 is to calculate M and r . Since firing is Poissonian with rate λ = f ∆t, the probability that a given afferent fires at least once in a given pattern subsection of length ∆t is p = 1 − e −f ∆t . Here, we consider independent patterns, i.e. with chancelevel overlap. Hence the probability that a given afferent fires at least once in at least one of the P pattern subsection is 1−(1−p) P . Thus the number of selected afferents M is on average: Finally, the expected effective input spike rate during the ∆t window is the expected total number of spikes, f N ∆t, divided by ∆t, thus: We note that the SN R scales with √ N . In the rest of this paper we used N = 10 4 afferents, which is in the biological range.

Numerical validations
We first checked if the variability of the SN R across Poisson realizations is small, and also if the approximation we made to compute the average SN R in Equation 6 is reasonable. To this aim, we generated 10 5 Poisson patterns, and computed M , r and the reduced SN R, snr = ( r − f M )/ M , for each of them (i.e. the right factor of the SN R in Equation 6, which is the only one that depends on the Poisson realization). As can be seen on Figure 2 left, M and r are strongly correlated, and the data points lie near a line which corresponds to nearly constant snr values (see the colored background). In other words, the snr does not change much for different Poisson pattern realizations and the average snr well represents the snr distribution even for the worst and best cases.
In addition, as can be seen on Figure 2 right, the average snr across different Poisson patterns is very close to the snr corresponding to the average-case scenario, i.e. M = M and r = r (as defined by Equations 7 and 8 respectively). Note that this Figure was done with relatively small values for the parameters P , ∆t and f (respectively 1, 2ms, and 1Hz). Our simulations indicate that when increasing these parameter values, the approximation becomes even better (data not shown).
Next, we verified the complete SN R formula (Eq. 6), which also includes v max , through numerical simulations. We used a clock-based approach, and integrated the LIF equation using the forward Euler method with a 0.1ms time bin. We used P = 1 and P = 5 patterns, and performed 100 simulations with different random Poisson patterns of duration L = 20ms with rate f = 5Hz. We chose ∆t = L = 20ms, i.e. the LIF was connected to all the afferents that emitted at least once during one of the patterns. In order to estimate V max , each pattern was presented 1000 times, every 400ms. The maximal jitter was T = 5ms. Between pattern presentations, the afferents fired according to a Poisson process, still with rate f = 5Hz, which allowed to estimate V noise and σ noise . We could thus compute the SN R from Equation 1 (and its standard deviation across the 100 simulations), which, as can be seen on Figure 3, matches very well the theoretical values, for P = 1 and 5. Note that the SN R standard deviation is small, which confirms that the average SN R, i.e SN R , represents well the individual ones.

Optimizing the SNR
We now want to optimize the SN R given by Equation 6, by tuning τ and ∆t. We also add the constraint τ f M ≥ 10 (large number of synaptic inputs), so that the distribution of V is approximately Gaussian [Burkitt, 2006]. Otherwise, it would be positively skewed 1 , thus a high SN R would not guarantee a low false alarm rate. We assume that 4 Simulations show that STDP can be close-to-optimal L is sufficiently large so that an upper bound for ∆t is not needed. We used the Matlab R2017a Optimization Toolbox (MathWorks Inc., Natick, MA, USA) to compute the optimum numerically. Figure 4 illustrates the results with P = 2. One can make the following observations (similar to our previous paper which was limited to P = 1 [Masquelier, 2017]): • Unless f and T are both high, the optimal τ and ∆t have the same order of magnitude (see Figure 4 left).
• Unless T is high (>10ms), or f is low (<1Hz), then these timescales should be relatively small (at most a few tens of ms; see Figure 4 middle). This means that even long patterns (hundreds of ms or more) are optimally detected by a coincidence detector working at a shorter timescale, and which thus ignores most of the patterns. One could have thought that using τ ∼ L, to integrate all the spikes from the pattern would be the best strategy. But a long τ also decreases the detector's temporal resolution, thus patterns and random inputs elicit more similar responses, decreasing the SN R.
Hence there is a trade-off, and it turns out that it is often more optimal to have τ < L, that is to use subpatterns as signatures for the whole patterns.
• Unsurprisingly, the optimal SN R decreases with T (see Figure 4 right). What is less trivial, is that it also decreases with f . In other words, sparse activity is preferable. We will come back to this point in the discussion.
What is the biological range for T , which corresponds to the spike time precision? Millisecond precision in cortex has been reported [Havenith et al., 2011, Kayser et al., 2010, Panzeri and Diamond, 2010. We are aware that other studies found poorer precision, but this could be due to uncontrolled variable or the use of inappropriate reference times [Masquelier, 2013].
In the rest of the paper we focus, as an example, on the point on the middle of the T × f plane -T = 3.2ms and f = 3.2Hz. When increasing P , the optimal τ and ∆t decrease (Fig. 5). Unsurprisingly, the resulting SN R also decreases, but only slowly. It thus remains acceptable for several tens of independent patterns (e.g. SN R ∼ 7 for P = 40).

Simulations show that STDP can be close-to-optimal
Next we investigated, through numerical simulations, if STDP could turn a LIF neuron into an optimal multi-pattern detector. More specifically, since STDP does not adjust the membrane time constant τ , we set it to the optimal value and investigated whether STDP could learn all the patterns with an optimal ∆t 2 . Here, unlike in the previous section, we had to introduce a threshold, in order to have postsynaptic spikes, which are required for STDP. As a result, the optimal V max , computed in the previous section, was never reached. Yet a high V max guarantees a low miss rate, and a low V noise guarantees a low false alarm rate. Optimizing the previously defined SN R thus makes sense. Again, we used a clock-based approach, and the forward Euler method with a 0.1ms time bin. The Matlab R2017a code for these simulations has been made available in ModelDB [Hines et al., 2004] at http://modeldb.yale.edu/244684.

Input spikes
The setup we used was similar to the one of our previous studies [Gilson et al., 2011, Masquelier, 2017, Masquelier et al., 2008. Between pattern presentations, the input spikes were gener-ated randomly with a homogeneous Poisson process with rate f . The P spike patterns with duration L = 100ms were generated only once using the same Poisson process (frozen noise). The pattern presentations occurred every 400ms (in previous studies, we demonstrated that irregular intervals did not matter [Gilson et al., 2011, Masquelier et al., 2008, so here regular intervals were used for simplicity). The P patterns were presented alternatively, over and over again. Figure 6 shows an example with P = 2 patterns. At each pattern presentation, all the spike times were shifted independently by some random jitters uniformly distributed over [−T, T ].

A LIF neuron with adaptive threshold
We simulated a LIF neuron connected to all of the N afferents with plastic synaptic weights w i ∈ [0, 1], thus obeying the following differential equation: where t ij is the time of the j th spike of afferent i. We used an adaptive threshold (unlike in our previous studies [Gilson et al., 2011, Masquelier, 2017, Masquelier et al., 2008, in which a fixed threshold was used). This adaptive threshold was increased by a fixed amount (1.8θ 0 ) at each postsynaptic spike, and then exponentially decayed towards its baseline value θ 0 with a time constant τ θ = 80ms. This is a simple, yet good model of cortical cells, in the sense that it predicts very well we plotted the neuron's potential as a function of time. Colored rectangles indicate pattern presentations. Next, we plotted the two spike patterns, coloring the spikes as a function of the corresponding synaptic weights: blue for low weight (0), purple for intermediate weight, and red for high weight (1). Initial weights were uniform (here at 0.7, so the initial color is close to red).
(Middle) During learning. Selectivity progressively emerges. (Bottom) After convergence. STDP has concentrated the weights on the afferents which fire at least once in at least one of the pattern subsections, located at the beginning of each pattern, and whose duration roughly matches the optimal ∆t (shown in green). This results in one postsynaptic spike each time either one of the two pattern is presented. Elsewhere both V noise and σ noise are low, so the SN R is high. In addition V noise roughly matches the theoretical value V opt noise (shown in green), corresponding to the optimal SN R. We also show in green V opt max , the theoretical optimal value for V max . However, the potential never reaches it, because the adaptive threshold is reached before.
the spikes elicited by a given input current Naud, 2009, Kobayashi et al., 2009]. Here, such an adaptive threshold is crucial to encourage the neuron to learn multiple patterns, as opposed to fire multiple successive spikes to the same pattern. Since the theory developed in the previous sections ignored the LIF threshold, using an adaptive one is not worse than a fixed one, in the sense that it does not make the theory less valid. We did not know which value for θ 0 could lead to the optimum. We thus performed and exhaustive search, using a geometric progression with a ratio of 2.5%.

Synaptic plasticity
Initial synaptic weights were all equal. Their value was computed so that V noise = θ + σ noise (leading to an initial firing rate of about 4Hz, see Figure 6 top). They then evolved in [0, 1] with all-to-all spike STDP. Yet, we only modeled the Long Term Potentiation part of STDP, ignoring its Long Term Depression (LTD) term. As in Song et al. [2000], we used a trace of presynaptic spikes at each synapse i, A i pre , which was incremented by δA pre at each presynaptic spike, and then exponentially decayed towards 0 with a time constant τ pre = 20ms. At each postsynaptic spike this trace was used for LTP at each synapse: pre . Here LTD was modeled by a simple homeostatic mechanism. At each postsynaptic spike, all synapses were depressed: w i → w i + w i (1 − w i )w out where w out < 0 is a fixed parameter [Kempter et al., 1999].
Note that for both LTP and LTD we used the multiplicative term w i (1 − w i ), in contrast with additive STDP, with which the ∆w is independent of the current weight value [Kempter et al., 1999, Song et al., 2000. This multiplicative term ensures that the weights remain in the range [0,1], and the weight dependence creates a soft bound effect: when a weight approaches a bound, weight changes tend toward zero. Here it was found to increase performance (convergence time and stability), in line with our previous studies [Kheradpisheh et al., 2016, 2018, Masquelier and Thorpe, 2007, Mozafari et al., 2017.
The ratio between LTP and LTD, that is between δA pre and w out is crucial: the higher, the more synapses are maximally potentiated (w = 1) after convergence. Here we chose to keep δA pre = 0.1 and to systematically vary w out , using again a geometric progression with a ratio of 2.5%.

Results
For each θ 0 × w out point, 100 simulations were performed with different random pattern realizations, and we computed the proportion of "optimal" ones (see below), and reported it in Table 1. After 12,000s of simulated time, the synaptic weights had all converged by saturation. That is synapses were either completely depressed (w = 0), or maximally potentiated (w = 1). A simulation was considered optimal if 1. all the patterns were learned, and 2. in an optimal way, that is if all patterns exhibited a subsection in which all spikes corresponded to maximally potentiated synapses (w = 1), and whose duration roughly matched the theoretical optimal ∆t. In practice, we used the total number of potentiated synapses as a proxy of the mean subsection duration 5 Discussion

10
Tab. 1: Performance as a function of the number of patterns P . The first four lines are computed from the theoretical optimum. The next two lines are the optimal values found through exhaustive search (see text). The last four lines are performance indicators, estimated during the last 100 presentations of each pattern. P learned is the mean number of "learned patterns", that is by convention patterns which elicit at least one postsynaptic spike. The following line is the mean hit rate for those patterns. The subsequent line gives the false alarm rate, but we never observed any here. Finally P(opt) is the proportion of optimal cases. P 5 10 20 40 ∆t opt (ms) 11 8.1 5.7 3.7 τ opt (ms) 8.9 6.8 5.6 5.1 M opt 1600 2300 3100 3800 SN R opt 31 20 12 6.7 θ 0 190 140 110 92 w out −6.2 10 −3 −6.3 10 −3 −6.5 10 −3 −6.7 10 −3 P learned 5 10 20 39.5 Hit rate (%) 98.9 98.6 97.9 96.5 False alarms (Hz) 0 0 0 0 P(opt) (%) 100 100 100 58 (since there is a non-ambiguous mapping between the two variables, given by Equation 7), and checked if this number matched the theoretical optimal M (Eq. 7) with a 5% margin.
Note that this second condition alone would be easy to satisfy: the total amount of potentiated synapses is determined by the LTP/LTD ratio which we adjusted by fine-tuning w out . However, satisfying the two conditions is harder, especially when P increases (Table 1).
It is worth mentioning that the learned subsections always corresponded to the beginning of the patterns, because STDP tracks back through them [Gilson et al., 2011, Masquelier et al., 2008], but this is irrelevant here since all the subsections are equivalent for the theory. Figure 6 shows an optimal simulation with P = 2 patterns.
As can be seen in Table 1, the proportion of optimal simulations decreases with P , as expected. But more surprisingly, several tens of patterns can be optimally learned with reasonably high probability. With P = 40 the probability of optimal simulations is only 58%, but the average number of learned patterns is high: 39.5! This means that nearly all patterns are learned in all simulations, yet sometimes in a suboptimal manner. Finally, Figure 7 shows that convergence time increases with P .
Yet one issue with having one neuron selective to multiple patterns is stability. If one of the learned pattern does not occur for a long period during which the other patterns occur many times, causing postsynaptic spikes, the unseen pattern will tend to be forgotten. This is not an issue with localist coding: if the learned pattern does not occur, the threshold is hardly ever reached so the weights are not modified, and the pattern is retained indefinitely, even if STDP is "on" all the time.
Another issue with distributed coding is how the readout is done, that is how the identity of the stimulus can be inferred from multiple neuron responses, given that each response is ambiguous? This is out of the scope of the current paper, but we suspect that STDP could again help. As shown in this study, each neuron equipped with STDP can learn to fire to multiple independent stimuli. Let's suppose that stimuli are shown one at a time. When stimulus A is shown, all the neurons that learned this stimulus (among others) will fire synchronously. Let us call S this set of neurons. A downstream neuron equipped with STDP could easily become selective to this synchronous volley of spikes from neurons in S [Brette, 2012]. With an appropriate threshold, this neuron would fire if and only if all the neurons in S have fired. Does that necessarily mean that A is there? Yes, if the intersection of the sets of stimuli learned by neurons in S only contains A. In the general case, the intersection is likely to be much smaller than the typical sets of stimuli learned by the S neurons, so much of the ambiguity should be resolved.
What could determine the set of patterns to which a neuron responds? Here, we used independent, unrelated, patterns (i.e. with chance-level overlap), and yet several of these patterns could be learned by a single neuron. Of course, patterns with more overlap would be easier to group. So in the presence of multiple postsynaptic neurons, each one would tend to learn a cluster of similar patterns. Another factor is the time at which the patterns are presented: those presented at the same period are more likely to be learned by the same neuron -a neuron which was still unselective at that period. Indeed, neurons equipped with STDP have some sort of critical period, before convergence, during which they can learn new pattern easily. Conversely, after convergence, neurons tend to fire if and only if the patterns they have learned are presented (Fig.6), and thus can hardly learn any new pattern. This is interesting, because patterns presented at the same period are likely to be somewhat related. For example, a neuron could fire to the different people you have met on your first day at work. In the presence of neurogenesis, newborn neurons could handle the learning of other patterns during the subsequent periods of your life. Finally, here we did not use any reward signal. But such a signal, if available, could modulate STDP (leading to some form of supervised learning), and encourage a given neuron to fire to a particular, meaningful, set of patterns [Mozafari et al., 2018a,b], as opposed to a random set like here. For e.g. a single neuron could learn to fire to any animal, even if dif-ferent animals cause very different sensory inputs.
Here the STDP rule we used always led to binary weights after learning. That is an afferent could be either selected or discarded. We thus could use our SN R calculations derived with binary weights, and checked that the selected set was optimal given the binary weight constraint. Further calculations in the Appendix suggest that removing such a constraint could lead to a modest increase in SN R, of about 10%. More research is needed to see if a multiplicative STDP rule, which does not converge towards binary weights [Gütig et al., 2003, van Rossum et al., 2000, could lead to the optimal graded weights.
Our theoretical study suggests, together with others [Brette, 2012, Gütig andSompolinsky, 2006], that coincidence detection is computationally powerful. In fact, it could be the main function of neurons [Abeles, 1982, König et al., 1996. In line with this proposal, neurons in vivo appear to be mainly fluctuation-driven, not mean-driven [Brette, 2012, Rossant et al., 2011. This is the case in particular in the balanced regime [Brette, 2015], which appears to be the prevalent regime in the brain [Denève and Machens, 2016]. Several other points suggest that coincidence detection is the main function of neurons. Firstly, strong feedforward inhibitory circuits throughout the central nervous system often shorten the neurons' effective integration windows [Bruno, 2011]. Secondly, the effective integration time constant in dendrites might be one order of magnitude shorter than the soma's one [König et al., 1996]. Finally, recent experiments indicate that a neuron's threshold quickly adapts to recent potential values [Fontaine et al., 2014, Mensi et al., 2016, Platkiewicz and Brette, 2011, so that only a sudden potential increase can trigger a postsynaptic spike. This enhances coincidence detection. It remains unclear if other spike time aspects such as ranks [Thorpe and Gautrais, 1998] also matter.
Our results show that lower firing rates lead to better signal-to-ratio. It is worth mentioning that mean firing rates are probably largely overestimated in the electrophysiological literature, because extracellular recordings -by far the most popular technique -are totally blind to cells that do not fire at all [Thorpe, 2011]. Even a cell that fire only a handful of spikes will be ignored, because spike sorting algorithms need tens of spikes from a given cell before they can create a new cluster corresponding to that cell. Furthermore, experimentalists tend to search for stimuli that elicit strong responses, and, when they can move the electrode(s), tend to look for most responsive cells, introducing strong selection biases. Mean firing rates, averaged across time and cells, are largely unknown, but they could be smaller than 1 Hz [Shoham et al., 2006]. It seems like coding is sparse: neurons only fire when they need to signal an important event, and that every spike matters [Wolfe et al., 2010].
Finally, we see an analogy between our theory, and the one of neural associative memory (NAM), in which an output (data) vector is produced by multiplying an input (address) vector by a weight matrix. Unlike NAM, our framework is dynamic, yet after learning, to a first approximation, our STDP neurons count the number of input spikes arriving through reinforced synapses in a short integration window, and each one outputs a 1 (i.e. a spike) if this count exceeds a threshold, and a 0 otherwise, leading to a binary output vector, much like in a binary NAM. It is thus unsurprising that sparsity is desirable both in our theory, and in NAMs [Palm, 2013].
should put strong weights on the synapses corresponding to the most recent pattern spikes, since these weights will increase V max more than V noise . Conversely, very old pattern spikes that fall outside the integration window (if any) should be associated to nil weights: any positive value would only increase V noise , not V max . But between those two extremes, it might be a good idea to use intermediate weight values.
To check this intuition, we used numerical optimizations using a simplified setup. We used a single pattern (P = 1), that was repeated in the absence of jitter (T = 0). We divided the pattern into n different periods ∆t 1 , ...∆t n (in reverse chronological order), each one corresponding to a different synaptic weight w 1 , ...w n (see Figure 8 left for an example with n = 2). More specifically: the M 1 afferents that fire in the ∆t 1 window are connected with weight w 1 . The M 2 afferents that fire in the ∆t 2 window, but not in the ∆t 1 one, are connected with weight w 2 . More generally, the M i afferents that fire in the ∆t i window, but not in the ∆t 1 ...∆t i−1 ones, are connected with weight w i .
With this simple set up, the SN R can be computed analytically. For example, if n = 2 (Fig. 8 left), we have: The asymptotic steady regimes for the two time windows are: Let's call V i the potential at the end of window ∆t i , and V n+1 = V noise . Then V max = V 1 can be computed iteratively: Furthermore [Burkitt, 2006], σ noise = τ f (w 2 1 M 1 + w 2 2 M 2 )/2.
So we have everything we need to compute the SN R.
So the SN R can be computed for any n, and, importantly, it is differentiable with respect to the w i . We can thus use efficient numerical methods to optimize these weights. Since scaling the weights does not change the SN R, we imposed w 1 = 1. Figure 8 right gives an example with n = 70. Here the ∆t i were all equal to 5τ /n, and we optimized the corresponding w i . We chose τ = 10ms, and f = 1, 5, and 10Hz. The gain w.r.t. binary weights for the SN R were modest: 10.5%, 9.6% and 8.9% respectively. As f tends towards 0, the optimal weights appears to converge towards e t/τ (even if we could not prove it): the f = 1Hz curve (solid blue) is almost identical to e t/τ (dashed red). Didactic example with n = 2 weight values: w 1 for all the afferents that fire in the ∆t 1 window, and w 2 < w 1 for all the afferents that fire in the ∆t 2 window but not in the ∆t 1 one. V ∞ 1 and V ∞ 2 are the asymptotic potentials for the two periods. V max can be computed from those two values (see text). (Right) Numerical optimization of the weights with n = 70. With small f , the optimal solution appears to be close to e t/τ .