Abstract

Polychronous neuronal group (PNG), a type of cell assembly, is one of the putative mechanisms for neural information representation. According to the reader-centric definition, some readout neurons can become selective to the information represented by polychronous neuronal groups under ongoing activity. Here, in computational models, we show that the frequently activated polychronous neuronal groups can be learned by readout neurons with joint weight-delay spike-timing-dependent plasticity. The identity of neurons in the group and their expected spike timing at millisecond scale can be recovered from the incoming weights and delays of the readout neurons. The detection performance can be further improved by two layers of readout neurons. In this way, the detection of polychronous neuronal groups becomes an intrinsic part of the network, and the readout neurons become differentiated members in the group to indicate whether subsets of the group have been activated according to their spike timing. The readout spikes representing this information can be used to analyze how PNGs interact with each other or propagate to downstream networks for higher-level processing.

1  Introduction

One of the fundamental goals of neuroscience is to understand how information is encoded and decoded in spike trains (Buonomano & Maass, 2009). A putative mechanism for information encoding and decoding is explained by the cell assembly theory, which is loosely characterized as a group of neurons having a strong mutual relationship and functioning as a representational entity. On one hand, information encoding is thought to be closely related to the dynamic formation of cell assemblies (Buzsáki, 2010). On the other hand, information decoding involves the activation of some members in the assembly according to a specific spatiotemporal pattern (Sutherland & McNaughton, 2000). The spike timing in the pattern is important for retrieving the desired piece of information, because only in this way do the downstream neurons receive synchronously arriving spikes and accumulate enough membrane potential to propagate the activity. A different spatiotemporal pattern leads to the activation of other assemblies, therefore retrieving a different piece of information.

Polychronous neuronal group (PNG) (Izhikevich, 2006) is one of the computational models for cell assembly. A PNG is defined as a group of neurons behaving as time-locked but not synchronized activities with millisecond precision due to the finite transmission delays in the network. In computer simulations, PNGs can be categorized into structural PNGs and activated (or dynamical) PNGs (Martinez & Paugam-Moisy, 2009). Structural PNGs are neuronal groups with prototypical time-locked firing patterns supported by the network anatomy (connectivity, weights, and delays). Activated PNGs refer to the actually activated structural PNGs during simulation, which can be partial and jittered compared to their prototypical patterns. Since spike-timing-dependent plasticity (STDP) constantly shapes the network, the structural PNGs dynamically emerge, grow, overlap, disappear, and compete for activation, leading to flexible and huge information representation capacity (Izhikevich, 2006; Szatmáry & Izhikevich, 2010).

However, it is difficult to detect PNGs from spike trains. An intuitive way is to enumerate all possible structural PNGs from network anatomy and use them as templates to find their actual activation in the spike trains (Izhikevich, 2006; Martinez & Paugam-Moisy, 2009). However, this approach is under the shadow of combination explosion. In contrast, response fingerprint (Guise, Knott, & Benuskova, 2014) is a nonenumeration-based algorithm that constructs a histogram of spikes from stimulation onset and then finds the peak time as the relative spike timing inside the PNG. However, constructing such a histogram requires the spikes to be aligned to an explicit time reference such as the stimulation onset, which can be unavailable in recurrent networks containing unknown recurrent stimulation. Further, these algorithms are feasible only in models, and currently no bioplausible rule such as STDP is used to detect PNGs.

To fill the gaps, we follow the reader-centric definition of cell assembly (Buzsáki, 2010) and consider PNGs as spatiotemporal patterns that can selectively activate downstream readout neurons. In fact, learning spatiotemporal patterns using readout neurons with STDP has already been studied under discrete spike volleys (Guyonneau, VanRullen, & Thorpe, 2005), continuous spike trains (Masquelier, Guyonneau, & Thorpe, 2009; Zenke, Agnes, & Gerstner, 2015) and oscillatory drive (Masquelier, Hugues, Deco, & Thorpe, 2009). The results are consistent: some readout neurons learn to fire if a particular pattern is presented in the spike trains. However, there is no way to know the spike timing inside the patterns according to the readout spikes, not to mention learning PNGs.

Here we show that the frequently activated structural PNGs can be learned using readout neurons with joint weight-delay STDP without enumeration or explicit time reference. During learning, the readout neurons adapt their incoming weights and delays jointly by strengthening and pulling together spikes arriving before the postsynaptic spike, and weakening and delaying others. Because of the delay plasticity, the delays become complementary to the spike timing in the PNG to maximize the arrival synchrony. In this way, the identity of neurons in the group and their expected spike timing, at millisecond scale, can be recovered from the incoming weights and delays of the readout neurons. The experiments were carried out under two realistic settings: (1) predefined patterns embedded in a continuous stream to simulate known structural PNGs (Masquelier, Guyonneau, et al., 2009) and (2) unknown self-organized structural PNGs from a sparsely connected recurrent network (Izhikevich, 2006). In the first experiment, it turned out that the pattern hit rate was higher and false alarm frequency lower than the conventional weight STDP. In the second experiment, we learned the frequently activated PNGs and showed improved performance using two layers of readout neurons.

Throughout this letter, all delays are assumed to be axonal. We use “weight STDP” to refer to the conventional STDP, which changes synaptic weights only, and use “joint STDP” for the introduced STDP, which changes synaptic weights and delays jointly. We use “structural PNGs” and “activated PNGs” as already defined; “learned PNG” or “receptive field of a readout neuron/spike” to refer to the expected neurons and their expected firing pattern recovered from the incoming weights and delays of a readout neuron; and “readout activation” or “activation of a readout neuron” to refer to the activation of the learned PNG associated with the readout neuron, which can be partial and jittered compared to the learned PNG. All simulations were done using NEural Simulation Tool (NEST) (Gewaltig & Diesmann, 2007). (The joint STDP synapse model is available at https://github.com/Hockey86/SNN.)

2  Joint Weight-Delay Spike-Timing-Dependent Plasticity

The joint weight-delay spike-timing-dependent plasticity learns weights indicating which neurons are contained in the PNG and delays complementary to the spike timing in the PNG to maximize the arrival synchrony. Like weight STDP, the learning rule is local; it requires only the information from the pair of pre- and postsynaptic neurons. It is also unsupervised, which strengthens the synapses carrying causal spikes while weakens the others without supervision.

Formally, the joint STDP is additive and uses nearest spike approximation.
formula
2.1
formula
2.2
Here wij and dij are the weight and axonal delay from a presynaptic neuron j to a postsynaptic neuron i; is the firing time of the neuron i; is the nearest firing time of the neuron j prior to ; and and are the non-hebbian homeostatic terms (Gilson, Masquelier, & Hugues, 2011).
The learning windows in equations 2.1 and 2.2 are (see Figure 1)
formula
2.3
formula
2.4
where and are the weight learning window time constants; and are the weight learning window amplitudes (); and are the effective weight learning window lengths; and are the delay learning window time constants on the same scale (); and are the delay learning window amplitudes on the same scale (); and are the weight-related gain offsets on the same scale. The weights are clipped at the boundaries , and the delays at . The values of the joint STDP parameters are shown in Table 1 in the appendix.
Figure 1:

(A) The weight learning window in equation 2.3. The shape is an exponential function with the negative branch dominating the positive branch (), which is the same with the weight STDP in the literature. Here the weight is a dimensionless scaling factor of postsynaptic potential (PSP). (B) The delay learning window D(s,w) in equation 2.4. The shape is a part of . Note that the amplitudes are controlled by weight-related gains and the effective length of the delay learning window is different from the one of the weight window.

Figure 1:

(A) The weight learning window in equation 2.3. The shape is an exponential function with the negative branch dominating the positive branch (), which is the same with the weight STDP in the literature. Here the weight is a dimensionless scaling factor of postsynaptic potential (PSP). (B) The delay learning window D(s,w) in equation 2.4. The shape is a part of . Note that the amplitudes are controlled by weight-related gains and the effective length of the delay learning window is different from the one of the weight window.

Under the dynamics, for readout neuron i, which learns a PNG involving neuron set P, its incoming weights and delays arrive at the following quasi-equilibrium state, with T becoming smaller slowly with speed proportional to (quasi-equilibrium) unless (fixed point) (the proof is in the appendix).
formula
2.5
formula
2.6

Here j runs over all the presynaptic neurons of neuron i; is the spike timing of the jth neuron in the PNG; is a subset of P containing synapses with spikes arriving before , defined as the postsynaptic firing time under and ; and T is an indeterminate variable to show that the delay is complementary to the spike timing in the PNG to maximize the arrival synchrony, .

Since Q is not unique, during learning, there are multiple quasi-equilibrium states under the dynamics, all of which are of the form of equation 2.5 and 2.6 with different Qs, that is, different subsets of P (the proof is in the appendix). Therefore the dynamics can be interpreted as each readout neuron, represented by a vector of its incoming weights and delays in the phase space, being attracted to multiple quasi-equilibrium states. On the other hand, the lateral inhibition between readout neurons decelerates or blocks the trajectories of some readout neurons. The noisy background activity pushes the trajectories away from quasi-equilibrium and makes them switch between multiple quasi-equilibrium states with different Qs (see Figure S1 in the online supplement). The switching behavior potentially leads to robust learning, since the readout neuron can switch to a different Q to avoid malfunctioned neurons. Note that the switching behavior is not observed in the weight STDP (see Figure S2 in the online supplement).

Since the joint STDP is additive, the part of delay STDP leads to bimodally distributed delays, which should be avoided for two reasons. First, PNGs can emerge and disappear due to the plasticity. Suppose a readout neuron has learned one PNG and its delays have become either or . After that PNG disappears, it is more difficult for the readout neuron to learn new PNGs since the delays are so divergent that some spikes arrive much later than others, being more difficult to accumulate enough postsynaptic potential (PSP) to fire and learn. Second, delays complementary to the PNG spike timing cannot be kept when the delays are either or . To avoid bimodally distributed delays, the weight-related gains and in the delay learning window act as gating factors (see proof part (A1.2) in appendix A.6). For strengthened synapses, the negative side of the delay learning window becomes small, so that the delays become smaller slowly, with speed proportional to and are complementary to the spike timing. For weakened synapses, the positive side of the delay learning window becomes small, so that the delays become larger slowly, with speed proportional to .

3  Learning Predefined Polychronous Neuronal Groups from Continuous Spike Trains

In realistic scenarios, the structural PNGs are generated from the network anatomy; therefore, their spike timing is unknown beforehand. However in the first experiment, three predefined spatiotemporal patterns embedded in continuous spike trains were used to simulate “known” structural PNGs. By knowing them exactly, we can validate the joint STDP by demonstrating the learned complementary delays, as well as computing pattern hit rate and false alarm frequency. Throughout this section, we use “patterns” to indicate structural PNGs since they are artificially predefined. Section 4 deals with learning unknown self-organized structural PNGs from the network anatomy.

3.1  Experiment Setup

The input continuous spike trains were similar to Masquelier, Guyonneau, et al. (2009) (see Figure 2), which had 100 spike sources and were 405 s long. There were three spatiotemporal spike patterns involving spike sources 1 to 50, 25 to 75, and 40 to 90, respectively. Each pattern was defined by a randomly selected 30 ms period from the spike trains and then copy-pasted at irregular intervals throughout the whole sequence, with a gaussian jitter of 1 ms standard deviation. The major difference with Masquelier, Guyonneau, et al. (2009) was that no pattern included all spikes from the selected period. Instead, only one randomly selected spike in each source was included in the pattern (called “wavefront”). The purpose of using wavefront patterns was that the delays complementary to the spike timing can be obtained without ambiguity. In fact, this made learning more difficult because each pattern became thinner than a whole period, so that the readout neurons were more difficult to accumulate enough PSPs from the patterns.

Figure 2:

One second of the input continuous spike trains. Three patterns, indicated by red, green, and blue, were embedded in the spike train; nonpattern spikes are gray. The bottom panel shows the population-averaged firing rates over 10 ms bins. The continuous spike trains had a firing rate at 38.5 Hz and standard deviation 5.4 Hz. From the firing rate alone, there was nothing characteristic to indicate the patterns, therefore, the information was encoded in the spike timing. The ratio of spikes involved in patterns was 14.3%, lower than 16.7% (1/3 pattern duration ratio × 1/2 sources in patterns) in Masquelier, Guyonneau, et al. (2009).

Figure 2:

One second of the input continuous spike trains. Three patterns, indicated by red, green, and blue, were embedded in the spike train; nonpattern spikes are gray. The bottom panel shows the population-averaged firing rates over 10 ms bins. The continuous spike trains had a firing rate at 38.5 Hz and standard deviation 5.4 Hz. From the firing rate alone, there was nothing characteristic to indicate the patterns, therefore, the information was encoded in the spike timing. The ratio of spikes involved in patterns was 14.3%, lower than 16.7% (1/3 pattern duration ratio × 1/2 sources in patterns) in Masquelier, Guyonneau, et al. (2009).

In the readout layer, there were 10 spike response model (SRM) neurons (Gerstner, Kistler, Naud, & Paninski, 2014). There was nonplastic lateral inhibition between readout neurons with no delay. To compare the performance between the weight STDP and joint STDP, we created two networks of the same structure, but the connections from spike sources to readout neurons were subject to different STDP types. For the joint STDP, the weights were initially , with delays at . For the weight STDP, the weights were initially , and the delays were constants at . The parameter values of the neuron model are shown in Table 3 in the appendix.

3.2  Network Dynamics under the Joint STDP

At the beginning of the simulation, one of the readout neurons first reached its firing threshold, where its incoming weights and delays were then updated by the joint STDP. For the weights, the update strengthened synapses with spikes arriving before the postsynaptic spike, while weakened others. For the delays, spikes arriving after the postsynaptic spike were pulled afterward to arrive even later, making it less possible to contribute to the postsynaptic spike. Spikes arriving before the postsynaptic spike were pulled forward and closer to each other to increase the synchrony, making it more possible to contribute to the postsynaptic spike. The pulling dynamics stemmed from the negative side of the delay learning window, which made the relative difference between spike arrival times of any synapses j and k smaller (, and (see proof part (A1.2) in appendix A.6).

At the end of the simulation, each pattern was optimally learned by some readout neuron(s), such as pattern 1 by readout 10 in Figure 3. Due to the joint STDP, the spikes from pattern 1 were pulled together and maximally strengthened. For other patterns, their spikes were scattered; therefore readout 10 could fire exclusively for pattern 1 (see Figure 3, right).

Figure 3:

(Top) The membrane potential trace of readout neuron 10 at the beginning and end of the simulation. Here, the rest potential is 0 mV for convenience and firing threshold mV. The top black vertical bars are the readout spikes. (Bottom) The presynaptic spike arrival times versus their weights. For each half of the figure, synaptic weights measured at the timing of the last readout spike in the respective half of the figure were used. At the beginning of the simulation, the arriving spikes from each pattern (color) were indistinguishable with nonpattern spikes (gray), leading to nonselective periodical firing of the readout neuron. At the end of the simulation, for this neuron that learned pattern 1, the spike arrival times for pattern 1 (red) were pulled together due to the delay STDP, leading to postsynaptic firing. In contrast, other patterns did not trigger postsynaptic firing due to the scattered arrival times.

Figure 3:

(Top) The membrane potential trace of readout neuron 10 at the beginning and end of the simulation. Here, the rest potential is 0 mV for convenience and firing threshold mV. The top black vertical bars are the readout spikes. (Bottom) The presynaptic spike arrival times versus their weights. For each half of the figure, synaptic weights measured at the timing of the last readout spike in the respective half of the figure were used. At the beginning of the simulation, the arriving spikes from each pattern (color) were indistinguishable with nonpattern spikes (gray), leading to nonselective periodical firing of the readout neuron. At the end of the simulation, for this neuron that learned pattern 1, the spike arrival times for pattern 1 (red) were pulled together due to the delay STDP, leading to postsynaptic firing. In contrast, other patterns did not trigger postsynaptic firing due to the scattered arrival times.

3.3  Learning Performance

The response latencies, hit rate, and false alarm frequency of the three patterns for both STDP types are shown in Figure 4. It turned out that the joint STDP had a higher pattern hit rate and lower false alarm frequency than the weight STDP for all three patterns. The reason again lay in the delay learning part. For the weight STDP, from Figure 5A we can see that it actually implemented delay selection, where synapses with similar spike timing were strengthened (with ), indicated by the horizontal lines. However, only a few synapses satisfied this condition; therefore, fewer synapses were strengthened for the weight STDP. This effect was more obvious for pattern 2 in Figure 5A in green. Due to fewer strengthened synapses, the firing of the readout neuron became less reliable, leading to a lower pattern hit rate and higher false alarm frequency. For the joint STDP, it explicitly pulled the spikes closer to maximize the arrival synchrony; therefore, more synapses were strengthened, and the readout response became more reliable.

Figure 4:

(A–C) The readout response latencies of the three patterns for the weight STDP. The latencies were computed with respect to the pattern occurrence times, with false alarm indicated by 0. (D–F) Similar to panels A–C, but for the joint STDP. We can observe a slowly decreasing trend due to the quasi-equilibrium state. The latencies for the weight STDP were lower because their delays were fixed at . (G–I) The hit rates of the three patterns under the weight and joint STDP, along with time. The values were computed using a 10 s sliding window with 5 s overlap. (J–L) The false alarm frequencies of the three patterns. The detailed method to compute pattern hit rate and false alarm frequency is available in the appendix.

Figure 4:

(A–C) The readout response latencies of the three patterns for the weight STDP. The latencies were computed with respect to the pattern occurrence times, with false alarm indicated by 0. (D–F) Similar to panels A–C, but for the joint STDP. We can observe a slowly decreasing trend due to the quasi-equilibrium state. The latencies for the weight STDP were lower because their delays were fixed at . (G–I) The hit rates of the three patterns under the weight and joint STDP, along with time. The values were computed using a 10 s sliding window with 5 s overlap. (J–L) The false alarm frequencies of the three patterns. The detailed method to compute pattern hit rate and false alarm frequency is available in the appendix.

Figure 5:

(A) Input and readout spike trains under the weight STDP. The readout spikes from the optimally learned readout neurons for each pattern are highlighted by colors and arrows. The horizontal lines in the input spike trains indicate the strengthened spike sources. (B) Input and readout spike trains under the joint STDP. The spike trains are the same as panel A but with a different timescale. The insets show the expected spike timing, represented by receptive fields with precision with respect to the colored readout spikes. (C, D) The weights at the end of the simulation under the weight and joint STDP, respectively. The dashed lines indicate the spike sources containing the pattern. (E) The pattern presynaptic arrival times for the strengthened synapses () under the joint STDP at the end of the simulation.

Figure 5:

(A) Input and readout spike trains under the weight STDP. The readout spikes from the optimally learned readout neurons for each pattern are highlighted by colors and arrows. The horizontal lines in the input spike trains indicate the strengthened spike sources. (B) Input and readout spike trains under the joint STDP. The spike trains are the same as panel A but with a different timescale. The insets show the expected spike timing, represented by receptive fields with precision with respect to the colored readout spikes. (C, D) The weights at the end of the simulation under the weight and joint STDP, respectively. The dashed lines indicate the spike sources containing the pattern. (E) The pattern presynaptic arrival times for the strengthened synapses () under the joint STDP at the end of the simulation.

We note that the performance of the weight STDP could be improved by having more spike sources involved in each pattern (50 in this work and 1000 in Masquelier, Guyonneau, et al., 2009). However, according to Figure 8 in Izhikevich (2006), for realistic structural PNGs generated from the network, their sizes are mostly less than 60, with the most probable value around 20. Therefore, for the purpose of learning PNGs, the joint STDP is more robust than the weight STDP.

3.4  Learned Patterns

The weights at the end of the simulation are shown in Figures 5C and 5D. Due to the additive rule, the weights became bimodally distributed for both weight STDP and joint STDP, from which we can identify which spike sources were contained in the pattern. The spike arrival times relative to each pattern occurrence under the joint STDP are shown in Figure 5E. Since the delays were complementary to the spike timing in the pattern, their arrival times should be the same (), up to the spike jitter.

Using the joint STDP, we recovered the identity of neurons in the group and their expected spike timing in the patterns, at millisecond scale, from the incoming weights and delays of the readout neurons, where the spike timing is unavailable in the weight STDP. The learned patterns are represented using receptive fields in Figure 5B. Formally, the receptive field with respect to a readout spike at is a set
formula
3.1
where is the precision of the receptive field, is approximated PSP rise time for multiple incoming spikes, which is typically shorter than . is the criterion for being strong synapse in this section.

4  Learning Self-Organized Polychronous Neuronal Groups from Recurrent Network

In contrast to using predefined spatiotemporal patterns as known structural PNGs in the previous experiment, in this experiment we dealt with learning unknown structural PNGs that were self-organized from a sparsely connected recurrent network, driven by the continuous spike train from the previous experiment. The difficulties in this experiment were that (1) the structural PNGs could emerge, grow, and disappear and (2) for the same structural PNG, different parts of it could be involved with some jitter when activated multiple times. We used two layers of readout neurons to improve learning performance and learn larger structural PNGs. The learning performance was measured by receptive field hit rate (RFHR)—the fraction of receptive field where spikes were actually received. To validate the learned PNGs, we compared them with structural PNGs enumerated from network anatomy (Izhikevich, 2006). We also validated the statistical significance of the readout activity using surrogate spike trains.

4.1  Experiment Setup

The main network was similar to Izhikevich (2006), consisting of 800 excitatory and 200 inhibitory SRM neurons. The first 100 excitatory neurons were input neurons that received only external spikes trains. Each excitatory neuron was connected to 100 neurons randomly chosen from the 900 noninput neurons, with synapses subject to the weight STDP. Each inhibitory neuron was connected to 100 neurons randomly chosen from the 700 noninput excitatory neurons, without STDP at the synapses. Initially all weights in the main network were . The delays in the main network were randomly generated from a uniform distribution ranging from to and kept fixed. Each neuron in the main network had a threshold so that at least three synchronous spikes from maximally strengthened synapses were needed to trigger postsynaptic firing.

There were two readout layers: the first layer had 150 SRM readout neurons connected to the 900 noninput neurons in the main network, and the second layer had 50 SRM readout neurons connected to the first layer. Within each layer, there was nonplastic lateral inhibition between readout neurons with no delay. All synapses from the main network to the first layer and from the first to second layer were subject to the joint STDP, with initial weights at and delays at . The amplitude of the joint STDP window was four times larger than the weight STDP in the main network, so that the readout dynamics is fast enough to adapt to the dynamic development of the structural PNGs. The detailed parameter values are in Tables 1 to 3 in the appendix.

The input spike trains were 1620 s long by repeating the spike trains from the previous experiment four times, so that the main network and both readout layers had longer time to evolve. The patterns in the input spike trains are marked with colors in Figures 6A and 7A. Note that the three patterns induced network activities involving different numbers of neurons in both the main network and the readout layers. This is recognized as a limitation in section 5.3.

Figure 6:

(A) The spike trains from the main network and the first readout layer from 1619.23 s to 1619.68 s of the simulation. The input spike trains (neuron 1–100) were distinguished from the other neurons. The two cyan-highlighted readout spikes indicate two activations of the same readout neuron, where the learned PNG was highlighted in orange/purple, that is, excitatory/inhibitory neurons (exclusive of inputs). The three patterns in the input spike trains were marked with red/green/blue. (B) The receptive fields of readout neuron 42 and their actual received spikes. A receptive field is said to receive a spike if there is a spike within around its expected timing (here ). (C) Causal firing relationship inside the actually received spikes in panel B. The causal relationship was obtained if a spike of a presynaptic neuron arrived less than before a postsynaptic spike. Note that the causal relationship from input to neurons in the main network is not shown, which caused some spikes with no or few incoming lines in the figure. (D–E) Similar to panels B and C, but activated at another time.

Figure 6:

(A) The spike trains from the main network and the first readout layer from 1619.23 s to 1619.68 s of the simulation. The input spike trains (neuron 1–100) were distinguished from the other neurons. The two cyan-highlighted readout spikes indicate two activations of the same readout neuron, where the learned PNG was highlighted in orange/purple, that is, excitatory/inhibitory neurons (exclusive of inputs). The three patterns in the input spike trains were marked with red/green/blue. (B) The receptive fields of readout neuron 42 and their actual received spikes. A receptive field is said to receive a spike if there is a spike within around its expected timing (here ). (C) Causal firing relationship inside the actually received spikes in panel B. The causal relationship was obtained if a spike of a presynaptic neuron arrived less than before a postsynaptic spike. Note that the causal relationship from input to neurons in the main network is not shown, which caused some spikes with no or few incoming lines in the figure. (D–E) Similar to panels B and C, but activated at another time.

4.2  The Learned PNGs

In Figure 6 we show a learned PNG from the first readout layer. When we compare Figures 6B and 6D, it turned out that for the same readout neuron, different parts of its receptive field received spikes when activated for different times. Nevertheless, the whole receptive field still captured a relatively complete group. This was because after multiple activations of the same readout neuron, on average, every neuron inside its receptive field had a larger chance of contributing to the readout spike due to their complementary delays.

The first readout layer usually learns small pieces of structural PNGs. This is because the size of learned PNGs is positively related to the firing threshold of the readout neurons, but neurons with high are less likely to fire and update their weights and delays. A natural way of learning larger PNGs is to use two layers of readout neurons, so that the first layer learns PNGs in the main network, while the second layer learns PNGs in the first layer (high-level PNG). In this way, a readout neuron in the second layer is selective to a set of learned PNGs in the main network, which can be conceptually considered as learning larger PNGs. We refer to these larger learned PNGs from the second layer as assembled PNGs (see Figure 7). An examination of the size and time span of the assembled PNGs revealed that they were indeed larger and lasted longer than the learned PNGs from the first layer (see Figure 8). In addition to learning larger PNGs, this hierarchical readout structure has other benefits, such as a higher RFHR (see Figure 9) and a higher percentage of learned PNGs consistent with network anatomy (see Figure 10).

Figure 7:

Similar to Figure 6. Only the differences are described here. (A) The spike trains from the main network and two readout layers. The two cyan-highlighted readout spikes in the second readout layer indicate two activations of the same readout neuron, where the spikes of the corresponding high-level PNG are also shown in cyan in the first readout layer. In the main network, we highlight spikes in the assembled PNG. (B) The receptive fields of the assembled PNG. (C–E) Similar to Figures 6C–E. Note again that the causal relationship from input to neurons in the main network was not shown, which caused some spikes with no or few incoming lines in the figure.

Figure 7:

Similar to Figure 6. Only the differences are described here. (A) The spike trains from the main network and two readout layers. The two cyan-highlighted readout spikes in the second readout layer indicate two activations of the same readout neuron, where the spikes of the corresponding high-level PNG are also shown in cyan in the first readout layer. In the main network, we highlight spikes in the assembled PNG. (B) The receptive fields of the assembled PNG. (C–E) Similar to Figures 6C–E. Note again that the causal relationship from input to neurons in the main network was not shown, which caused some spikes with no or few incoming lines in the figure.

Figure 8:

(A) The distribution of the PNG sizes (#neuron) from the first readout layer and the assembled PNGs. Here, only the sufficiently specific readout neurons were considered—readout neurons that fired within 2 s before the end of the simulation and had a receptive field hit rate (RFHR) not less than the median of all readout neurons in the same layer. For the second layer, the RFHRs of the assembled PNGs were used to determine sufficiently specific readout neurons. Note that the sizes of the learned PNGs were on average larger (“wider”) than those in Figure 8 in Izhikevich (2006); in our case, at least three synchronous spikes were needed to trigger postsynaptic firing, while it is two in Izhikevich (2006). We list some statistics here. Layer 1: max: 94; min: 5; median: 44; standard deviation: 17.4. Assembled: max: 106; min: 39; median: 73.5; standard deviation: 16. (B) The distribution of the learned PNG time spans. Layer 1: max: 18 ms; min: 7 ms; median: 13 ms; standard deviation: 3.5 ms. Assembled: max: 27 ms; min: 19 ms; median: 23 ms; standard deviation: 2.2 ms. Note that the time span of the learned PNGs from the first layer cannot exceed , because they were recovered from the incoming delays of the readout neurons.

Figure 8:

(A) The distribution of the PNG sizes (#neuron) from the first readout layer and the assembled PNGs. Here, only the sufficiently specific readout neurons were considered—readout neurons that fired within 2 s before the end of the simulation and had a receptive field hit rate (RFHR) not less than the median of all readout neurons in the same layer. For the second layer, the RFHRs of the assembled PNGs were used to determine sufficiently specific readout neurons. Note that the sizes of the learned PNGs were on average larger (“wider”) than those in Figure 8 in Izhikevich (2006); in our case, at least three synchronous spikes were needed to trigger postsynaptic firing, while it is two in Izhikevich (2006). We list some statistics here. Layer 1: max: 94; min: 5; median: 44; standard deviation: 17.4. Assembled: max: 106; min: 39; median: 73.5; standard deviation: 16. (B) The distribution of the learned PNG time spans. Layer 1: max: 18 ms; min: 7 ms; median: 13 ms; standard deviation: 3.5 ms. Assembled: max: 27 ms; min: 19 ms; median: 23 ms; standard deviation: 2.2 ms. Note that the time span of the learned PNGs from the first layer cannot exceed , because they were recovered from the incoming delays of the readout neurons.

Figure 9:

(A) RFHRs from readout layer 1, 2 and the assembled PNGs, computed along with the simulation every 40.5 s. (B) The box plot of the overall RFHRs over the whole simulation.

Figure 9:

(A) RFHRs from readout layer 1, 2 and the assembled PNGs, computed along with the simulation every 40.5 s. (B) The box plot of the overall RFHRs over the whole simulation.

Figure 10:

(A) The mapping between structural PNGs enumerated from network anatomy (top rectangle) and learned PNGs from the first readout layer (bottom rectangle). During the last 2 s of simulation, 249 (23.3%) of the structural PNGs were activated. (B) The mapping between structural PNGs (top rectangle) and the assembled PNGs from the second readout layer (bottom rectangle). (C) Comparing the number of activated PNGs mapped and not mapped to learned PNGs from the first readout layer during the last 2 s of simulation. Mann-Whitney U test revealed a significant difference at p-value (sample size: 79 mapped, 170 not mapped). (D) Similar to panel C but for the second layer. Mann-Whitney U test revealed a significant difference at p-value (sample size: 65 mapped, 184 not mapped).

Figure 10:

(A) The mapping between structural PNGs enumerated from network anatomy (top rectangle) and learned PNGs from the first readout layer (bottom rectangle). During the last 2 s of simulation, 249 (23.3%) of the structural PNGs were activated. (B) The mapping between structural PNGs (top rectangle) and the assembled PNGs from the second readout layer (bottom rectangle). (C) Comparing the number of activated PNGs mapped and not mapped to learned PNGs from the first readout layer during the last 2 s of simulation. Mann-Whitney U test revealed a significant difference at p-value (sample size: 79 mapped, 170 not mapped). (D) Similar to panel C but for the second layer. Mann-Whitney U test revealed a significant difference at p-value (sample size: 65 mapped, 184 not mapped).

4.3  Learning Performance

The learning performance was quantitatively measured by receptive field hit rate (RFHR), the fraction of receptive field where spikes were actually received, as shown in Figure 9. The RFHR was computed in the following way. First, for each readout spike, its receptive field was generated according to equation 3.1 with . And its RFHR was computed as the fraction of the receptive fields, which actually received a spike within around the expected timing. Then for each readout neuron, we computed the average RFHR of all readout spikes fired within the last 2 s in each 40.5 s interval. And each point in Figure 9A was computed as the maximum average RFHR of all readout neurons from the same layer. The reason for using only a short period (2 s) before the end of each interval was to avoid incorrect RFHRs due to the constantly updating readout weights and delays under the joint STDP.

For the first readout layer, at the beginning of the simulation (0–100 s), none of the readout neurons had more than five strong incoming weights (), so that we did not compute the RFHR at that period. After the beginning period, the RFHR rapidly grew and fluctuated around median value 0.56, which meant that the readout neurons could update their incoming weights and delays along with the dynamic development of the structural PNGs, so that on average, more than half of the receptive fields actually received spikes. The RFHR was not full because of the noisy spikes in the input, which generated transiently activated PNGs that interacted with the persistent PNGs, resulting in missing spikes in the receptive fields.

For the second readout layer, its RFHR grew rapidly before 200 s, fluctuated at 200 to 1000 s, and then increased to near 1. The second layer had higher RFHR than the first layer, with median value 0.86. The higher RFHR meant that the second layer received a cleaner firing pattern from the first readout layer, compared to the main network. In fact, each readout layer acted as a noise filter, where the useful information was carried by the readout spikes while the noise was reduced. A Wilcoxon test (nonparametric version of paired t-test because of the nongaussianity) confirmed the statistically significant difference between the RFHRs of two layers at p-value .

We also computed the RFHR of assembled PNGs, which fluctuated around 0.3 at 100 s to 1000 s and grew to around 0.4 after 1000 s, resulting in the median value 0.37. Its value was approximately the product of the RFHRs of the two layers; therefore, it was lower than them on average.

4.4  Validation of the Learned PNGs

Having examined the performance of learning self-organized structural PNGs, we next validated the learned PNGs by (1) finding the mapping between structural PNGs enumerated from network anatomy (see Figure 10) and (2) checking the statistical significance of the readout spikes to confirm that they did not appear by chance with the varying firing rate (see Figure 11).

Figure 11:

(A) The distribution of readout spike numbers in the first layer during the last 2 s of the simulation when using actual network and surrogate spike trains as the input. The Mann-Whitney U test revealed a statistically significant difference between them. (B) The distribution for the second readout layer. The Mann-Whitney U test confirmed a statistically more significant difference than the first layer.

Figure 11:

(A) The distribution of readout spike numbers in the first layer during the last 2 s of the simulation when using actual network and surrogate spike trains as the input. The Mann-Whitney U test revealed a statistically significant difference between them. (B) The distribution for the second readout layer. The Mann-Whitney U test confirmed a statistically more significant difference than the first layer.

To find the structural PNGs in the main network, we enumerated all three trigger neurons in the 100 input neurons, which had combinations (Izhikevich, 2006; Martinez & Paugam-Moisy, 2009). For each combination, we enumerated all 900 noninput neurons, which were connected to all three trigger neurons, as the fourth neuron. The fourth neuron determined the firing pattern of the three trigger neurons according to the complement of their delays, so that the three trigger spikes could arrive at the fourth neuron at the same time to ensure postsynaptic firing. The pattern was fed to the network and was run for 150 ms. For each combination, network responses containing at least 10 spikes were counted as structural PNGs. In total, 1069 structural PNGs were found. Learned PNGs from the first readout layer were obtained from the incoming weights and delays of sufficiently specific readout neurons—that is, readout neurons fired within 2 s before the end of the simulation and had an RFHR not less than the median of all readout neurons in the same layer. For the second readout layer, we used the assembled PNGs and applied the same procedure to identify sufficiently specific readout neurons. The reason of mapping only sufficiently specific readout neurons traced back to the lateral inhibition between readout neurons, so inevitably some readout neurons were inhibited by others and not able to learn. These readout neurons had incorrect receptive fields indicated by low RFHR, which were meaningless to be mapped to structural PNGs. The detailed mapping method is described in the appendix.

The results in Figure 10 indicated that only a small portion of structural PNGs from network anatomy was mapped. An obvious reason was that the maximum number was limited by the number of readout neurons in both layers (150 and 50). Another reason was that the structural PNGs were obtained by enumeration from network anatomy, with no guarantee of being actually activated. It turned out that 249 (23.3%) of them were activated in the last 2 s of simulation.

For the learned PNGs, all of them were mapped to activated structural PNGs. Further, in Figures 10C and 10D, the mapped activated PNGs were activated more frequently than those not mapped, confirmed by Mann-Whitney U test. It suggested that the readout neurons tended to learn the frequently activated structural PNGs, which were informative since they represented the inputs. In other words, the readout neuron became selective to the information represented by structural PNGs under ongoing activity.

For the first readout layer, 40 sufficiently specific readout neurons were mapped to 79 structural PNGs. There were 15 not mapped learned PNGs, because the structural PNGs were obtained from enumerating all 3 trigger neurons in the 100 input neurons, which did not cover all structural PNGs, such as those triggered by 4 neurons. A complete enumeration led to combination explosion, while the joint STDP approach did not face this problem.

For the second readout layer, 21 sufficiently specific readout neurons were mapped to 65 structural PNGs, and 3 were not mapped. Comparing two readout layers, two main differences were found. First, the second layer had a higher percentage of mapped neurons (layer 1: , layer 2: ). This was because of the cleaner firing pattern in the first readout layer compared to the main network, so that more readout neurons fired early without being laterally inhibited and learned the correct receptive fields. Second, less structural PNGs were mapped to the assembled PNGs in the second layer (layer 1: 79, layer 2: 65). The reason was that the learned PNGs from the first layer were smaller and did not last as long as the assembled ones (see Figure 8), so that the mapping condition was more relaxed, which required only small matched portions in the structural PNGs. The mapping condition for assembled PNGs was stronger, which required larger matched portions, leading to less mapped structural PNGs.

Finally, we checked the statistical significance of the readout spikes to confirm that they did not appear by chance with the varying firing rate. In Figure 11, we compared the readout spike numbers from both layers during the last 2 s of the simulation, when using actual and surrogate spike trains as the input, where the latter were generated by inverting the time of the actual network spike trains (Izhikevich, 2006). In this way, the important spike train statistics were preserved, while only the relative spike timing was changed. The reason for not using a longer period was to avoid incorrect results due to the constantly updating readout weights and delays under joint STDP. Using the Mann-Whitney U test (nonparametric version of t-test because of the nongaussianity), the two-sided p-value showed statistical significance for both layers. The second layer had more significant difference, again due to the noise reduction of the first layer, so that the receptive fields fit the actual spike trains better and the surrogate spike trains less well. Also note that some readout neurons fired many times even for the surrogate data because of the relatively low firing threshold (see Table 3 in the appendix).

5  Discussion

5.1  A PNG-Based Neural Encoding and Decoding Framework

As a cell assembly model, PNG is considered to be one of the powerful, but elusively difficult to detect, neural encoding and decoding methods. We have shown that the joint STDP allowed learning the spike timing inside PNGs at millisecond scale under two realistic settings: artificially generated continuous spike stream and self-organized network activity. In line with Buzsaki (2010), that “the reader is not necessarily an independent, isolated unit,” the detection in our work becomes an intrinsic part of the network and the readout neurons become differentiated members in the PNG to indicate whether subsets of the PNG have been activated according to their spike timing. It opens the possibilities of studying neural encoding and decoding based on PNG, such as analyzing the interaction between multiple PNGs and studying the representation ability of PNGs as a mechanism of memory (Szatmáry & Izhikevich, 2010). Further, the PNG readout neurons act as a noise filter and decode useful information into spikes, which propagate to downstream networks for higher-level processing or to motor cortex to close the loop (Tetzlaff, Dasgupta, Kulvicius, & Wörgötter, 2015). These possibilities form a PNG-based neural encoding and decoding framework where the joint STDP plays an important role (see Figure 12).

Figure 12:

A PNG-based neural encoding and decoding framework.

Figure 12:

A PNG-based neural encoding and decoding framework.

5.2  Related Work in the Literature

Hüning, Glünder, and Palm (1998) also maximized the arrival synchrony by delay plasticity. Our work differs with them from two aspects. First, in their work, the delay learning happens continuously during or at the end of an interval when the sum of PSPs is higher than a learning threshold. In our work, the delay learning happens at the beginning of that interval, when the PSP rises across the threshold and triggers a postsynaptic spike. Second, in their work, the learning rule is to maximize the peak of the sum of PSPs by gradient descent, so that the spikes are pulled into synchrony toward the mean arrival time. In our work, the joint STDP maximizes arrival synchrony by minimizing the relative difference between spike arrival times and pulls them forward since for .

Eurich et al. (2000) also suggested a joint use of weight and delay STDP. However, their learning windows are different from ours. They stated that synchronous spike arrival (i.e., complementary delays) is a stable solution when for and for . In contrast, our delay learning window is reversed with for and for . In fact, without the weight-related gains, our delay learning is indeed unstable. With the help of weight-related gains, it actually goes into quasi-equilibrium states where the delays become smaller slowly with speed proportional to . Before reaching , the relative differences between delays are kept (see the appendix). In fact, is small enough compared to , and it is not possible for biological brains to encounter exactly the same pattern for many times since STDP constantly shapes the network. As a result, for realistic cases, the dynamics is at quasi-equilibrium. In addition, for leads to reduced latency for strengthened synapses, a functional consequence consistent with the weight STDP.

In the literature, there are some arguments on the necessity of the delay STDP according to Senn, Schneider, and Ruf (2002), where they suggested that slow, unbiased fluctuations in the delay and unreliable synapse can implicitly drift the delays to maximize the arrival synchrony. To compare, first, in their work, the arrival times are pulled to the point , which maximizes the bi-alpha-shaped weight learning window (). However, we use an exponential-shaped window where . Second, their approach is stochastic, where the drift to optimality is a result of slow evolutionary advantage rather than an explicit delay plasticity that adapts to the development of PNGs faster.

There is some biological evidence to support delay STDP. Research in the calyx of Held in the mammalian auditory central nervous system reveals that activity-dependent modification of the synaptic transmission delay exists together with weight plasticity (Tolnai, Englitz, Scholbach, Jost, & Rübsamen, 2009; Bakkum, Chao, & Potter, 2008), which is found to depend on the strength of local transients (Bollmann & Sakmann, 2005), although on a smaller scale compared to upstream factors (Lorteije & Borst, 2011).

5.3  Limitations

There are two major limitations in our approach. First, the structural PNGs to be learned must be activated frequently so that the distortion on the readout weights and delays caused by the random noise between two successive activations can be overcome. In the experiment in section 3, there were 4461 instantiations of the three patterns (known structural PNGs) in the 405 s spike trains, implying that each same pattern appeared about every 272 ms on average. In the experiment in section 4, since the input spike trains were taken from the experiment in section 3, which contained both patterns and random noises, the network activity contained both PNGs activated by the patterns and the random noise. For PNGs activated by the patterns, their activations had the same period, 272 ms, which enabled the readout neurons to develop receptive fields with stable RFHRs. For PNGs activated by the random noise, their activations were less persistent and less frequent, which could be activated every 1 s or longer according to Izhikevich (2006). The PNGs activated by the random noise transiently interacted with PNGs activated by the patterns and created fluctuations in the RFHRs.

Second, the readout neurons have a preference to learn larger structural PNGs when multiple PNGs with different sizes are presented (see Figures 6A and 7A). And the readout neurons could fire unreliably or become silent for the small structural PNGs. This is because larger structural PNGs bring more spikes to the readout neuron, making it more possible to fire and learn. Therefore, the readout firing threshold plays an important role in the joint STDP, which should be fine-tuned. In fact, in biologically detailed neuron models such as the Hodgkin-Huxley model, there is no unique threshold and the transition to firing is smooth, which depends not only on the amplitude of inputs but also on the speed of the increase (Gerstner et al., 2014). An adaptive firing threshold should be used to address this limitation.

Appendix

A.1  Parameter Values

Table 1:
Joint STDP Parameters.
SymbolDescriptionValueUnit
,  Weight learning window time constant 16.8, 33.7 ms 
,  Weight learning window amplitude Section 3: ,  – 
  Section 4: ,   
,  Effective weight window length ,  ms 
,  Non-Hebbian homeostatic term Section 3: , 0 – 
  Section 4: ,   
,  Delay learning window time constant 16.8, 16.8 ms 
,  Delay learning window amplitude ,  – 
,  Weight-related gain offset ,  – 
 Maximum weight – 
,  Minimum/maximum delay 1, 35 ms 
SymbolDescriptionValueUnit
,  Weight learning window time constant 16.8, 33.7 ms 
,  Weight learning window amplitude Section 3: ,  – 
  Section 4: ,   
,  Effective weight window length ,  ms 
,  Non-Hebbian homeostatic term Section 3: , 0 – 
  Section 4: ,   
,  Delay learning window time constant 16.8, 16.8 ms 
,  Delay learning window amplitude ,  – 
,  Weight-related gain offset ,  – 
 Maximum weight – 
,  Minimum/maximum delay 1, 35 ms 
Table 2:
Weight STDP Parameters in the Network in Section 4.
SymbolDescriptionValueUnit
,  Weight learning window time constant 16.8, 33.7 ms 
,  Weight learning window amplitude Section 3: ,  – 
  Section 4: ,   
,  Effective weight window length ,  ms 
,  Non-Hebbian homeostatic term ,  – 
 Maximum weight – 
,  Minimum/maximum delay 1, 20 ms 
SymbolDescriptionValueUnit
,  Weight learning window time constant 16.8, 33.7 ms 
,  Weight learning window amplitude Section 3: ,  – 
  Section 4: ,   
,  Effective weight window length ,  ms 
,  Non-Hebbian homeostatic term ,  – 
 Maximum weight – 
,  Minimum/maximum delay 1, 20 ms 
Table 3:
SRM Neuron Model Parameters.
SymbolDescriptionValueUnit
,  PSP time constant 4.6, 4.6 ms 
 Refractory time constant 12 ms 
 Absolute refractory period ms 
,  Resting and reset potential 0,  mV 
 Firing threshold Section 3 mV 
     10 (weight STDP)  
     15 (joint STDP)  
  Section 4  
     2.6 (network)  
     9 (readout layer 1)  
     3 (readout layer 2)  
Readout lateral inhibition weight Section 3 – 
      (weight STDP)  
      (joint STDP)  
  Section 4  
      (readout layer 1)  
      (readout layer 2)  
SymbolDescriptionValueUnit
,  PSP time constant 4.6, 4.6 ms 
 Refractory time constant 12 ms 
 Absolute refractory period ms 
,  Resting and reset potential 0,  mV 
 Firing threshold Section 3 mV 
     10 (weight STDP)  
     15 (joint STDP)  
  Section 4  
     2.6 (network)  
     9 (readout layer 1)  
     3 (readout layer 2)  
Readout lateral inhibition weight Section 3 – 
      (weight STDP)  
      (joint STDP)  
  Section 4  
      (readout layer 1)  
      (readout layer 2)  

Note: The same model for readout and network neurons.

A.2  Fine-Tuning Parameters

In section 3, the readout firing threshold was selected from a list [0.1, 0.15]() (100 is the number of spike sources); readout lateral inhibition weight from [0.1, 0.15, 0.2, 0.25, 0.3](); non-Hebbian homeostatic term from [, , ]; and ratio between delay and weight learning rates from [, , , , ]. For each combination of parameter values, the pattern hit rate H and false alarm frequency F were computed over the last 60 s of the simulation for all 10 readout neurons and for three patterns. For each pattern, the readout neuron with the highest was selected; therefore, we had three values for three patterns, respectively. The mean of these three values was used as the final score to find the best parameter combination.

In section 4, the parameters of the first readout layer were fine-tuned; then we kept the best value and fine-tuned the parameters of the second layer. For the first layer, was selected from a list [0.01, 0.02]() (900 is the number of noninput neurons in the network); readout lateral inhibition weight from [0.1, 0.2](); and from [, , ]. For the second layer, was selected from a list [0.02, 0.03]() (150 is the number of neurons in the first readout layer) and other parameters from the same range as the first layer. The criterion for determining the best parameter combination is the RFHR for the assembled PNGs from the second layer because it measures the overall performance of both layers.

A.3  Input Spike Trains

The way to generate the input spike trains in section 3 was similar to Masquelier, Guyonneau, et al. (2009). The spike trains had 100 spike sources and were 405 s long. Because of memory constraints, for each spike source, a 135 s sequence was first generated and repeated three times. This repetition had no impact on the results.

Each spike source independently emitted spikes according to an inhomogeneous Poisson process with the instantaneous firing rate r within [0, 20]Hz and maximum change rate, allowing an increase from 0 to 20 Hz in 30 ms and vice versa. For each spike source at each time step , it had a probability of to emit a spike at a random time within this time bin; then r was updated by and clipped in [0, 20]Hz, where was updated by a random number from . The second-order update was to ensure the smoothness of r with time. In addition, to make sure the patterns defined in the following paragraphs contained at least one spike from each spike source, a spike was added whenever a source had been silent for more than 30 ms, which resulted in a minimum firing rate of 33.3 Hz.

Three patterns were embedded in the spike trains. The pattern occurrence times were randomly placed throughout the sequence, so that each pattern had equal occurrence possibility, the total duration of patterns was one-third of the whole simulation, and the minimum interval was 60 ms. Finally, the number of instantiations of the first pattern was 1434, second 1569, and third 1458—altogether 4461. The spike sources in each pattern were arbitrarily chosen and grouped together for better visualization.

The major difference from Masquelier, Guyonneau et al. (2009) was that each pattern contained only one spike from each spike source. In our case, the patterns were defined by first randomly selecting a 30 ms period from the spike train and then randomly taking one spike of each source within this period (“wavefront”). The reason for using the 30 ms time span was that it was shorter than so that the complementary delays were possible to cover the whole pattern. The reason for using wavefront patterns was that the delays complementary to the spike timing can be obtained without ambiguity.

Once the patterns were defined, they were copy-pasted according to the pattern occurrence times. A gaussian jitter with 1 ms standard deviation was added before pasting. To keep the firing rate stable, the closest spike in the pasted location was removed before pasting. Finally, the spike times were rounded to the precision of with duplicate spikes removed.

The input spike trains in section 4 were generated by repeating the input spike trains in section 3 for four times, which were . They were then fed to the 100 input neurons to generate the network activity.

A.4  Latency, Pattern Hit Rate, and False Alarm Frequency in Section 3

Before computing the latency, we first defined a response interval for each pattern so that the readout spike falling into it was treated as a valid response. For readout neurons under weight STDP, since the delays from spike sources to readout neurons were constant at , the response interval was after each pattern. For readout neurons under joint STDP, since the delays were variant, the response interval started at the th early spike of pattern i plus after each pattern, where was the minimum number of maximally strengthened spikes to trigger the postsynaptic spike. The response interval ended at the half of pattern duration () plus after each pattern.

The latency was computed for each readout spike, by first finding its pattern response interval and then subtracting the pattern occurrence time. If it was outside any response interval, it was counted as a false alarm, with the latency denoted as 0. The false alarm frequency (Hz) was computed using the number of false alarms divided by a period T. The pattern hit rate was computed as the fraction of patterns with readout spike within its response interval over all patterns during the period T. Figure 4 was generated using a sliding window of and overlap.

A.5  Mapping between Learned PNGs and Structural PNGs in Section 4

The mapping was done using static readout neurons (Sun, Yang, Sourina, & Huang, 2015). Specifically, for each structural PNG, we created a readout neuron with static connections (no STDP) to all its member neurons, where the incoming weights were and the delays were . Here, is the spike timing of the ith neuron inside structural PNG. In this way, the static readout neuron receives maximum PSP when the arriving spikes follow the spike timing of its corresponding PNG. The network was connected with 1069 static readout neurons (1069 structural PNGs); then we fed each learned PNG to the network. If a static readout neuron fired for that learned PNG, it was mapped since it shared common spike timing with some jitter, where is the firing threshold of the static readout neuron. The amount of tolerable jitter for mapping was controlled by both PSP time constant and of the static readout neuron. Here, and , where N1 is the size of the learned PNG; N2 is the size of the structural PNG; and 30 mV is to avoid an overly high firing threshold when , such as cyclic PNGs. The static readout neurons were also used to detect activated PNGs in Figure 10, where and .

A.6  Proof for Equations 2.5 and 2.6

In this section, we give a proof for the quasi-equilibrium weights and delays in equations 2.5 and 2.6 under the dynamics of equations 2.3 and 2.4 and the following conditions.

  • There is only one structural PNG involving neuron set with spike timing at .

  • There is only one readout neuron i, so that there is no lateral inhibition.

  • We analyze only the case where the readout spike is caused by the structural PNG, not by the noise.

The proof also holds for multiple structural PNGs and readout neurons. This is because with the lateral inhibition, some readout neurons have the chance to learn only one particular structural PNG, thus reducing to the case of single structural PNG and readout neuron. The proof also holds for mediate level of random noise between multiple structural PNG activations. The random noise can make the weights slowly drift to 0 with speed proportional to , and the delays slowly drift to with speed proportional to . If the noise density and duration are mediate so that the update caused by the structural PNGs can overcome these drifts, the readout dynamics is mainly PNG driven (also mentioned in section 5.3). These cases are validated in sections 3 and 4.

The proof is divided into part A and part B:

  1. is a quasi-equilibrium state.

  2. All quasi-equilibrium states are of the form of .

A.6.1  Part A of the Proof

is a quasi-equilibrium state.

Suppose is distorted by a small amount,
formula
A.1
formula
A.2
where and are small compared to and , respectively. The postsynaptic spike time is changed to under the distorted weights and delays .

We show that under a small distortion, recover to , but with T becoming smaller slowly with speed proportional to (quasi-equilibrium, not fixed point) unless (fixed point).

(A1). When .

(A1.1). For the weight, we have
formula
A.3

Under the assumption that the distortion is small, we have ; therefore, and recovers to for . We note that holds only for small distortion.

(A1.2). For the delay, we note that the following four statements are equivalent:

  • dij is of the form .

  • The delay dij is complementary to the spike timing in the structural PNG.

  • The spike arrival times of any synapses j and k are the same—i‥e. , for , and .

  • The relative difference between spike arrival times of any synapses j and k is 0— for , and .

Therefore we equivalently prove the last statement. Formally, after the distortion, recovers to 0.

Without loss of generality, assuming , we have
formula
A.4
where
formula
A.5
formula
A.6
formula
A.7
formula
A.8
formula
A.9
formula
A.10

Therefore .

Continuing equation A.4:
formula
A.11
formula
A.12
formula
A.13

After the update, becomes smaller and finally recovers to 0; therefore, recovers to the form . The recovery speed is faster if is smaller, that is, the distortion is larger or is larger. Once again, if the distortion becomes large enough so that in equations A.8 and A.9 does not hold, the dynamics goes to another state. Note that in equation A.13, the equality holds only when holds for all , which corresponds to no distortion and no weight-related gain offset .

However, is not a fixed point because
formula
A.14
where .

In this case, becomes smaller slowly with speed proportional to and keeps the form , or equivalently T becomes smaller slowly. After learning many activations of the same structural PNG, reaches , at which the form cannot be kept anymore. However, is small enough compared to , and it is not possible for biological brains to encounter exactly the same input many times since STDP constantly shapes the network. On the other hand, if , it shuts down the dynamics with , therefore leading to a fixed point. A similar conclusion that dij is at quasi-equilibrium unless also holds for and . We note that without the weight-related gains, , so becomes smaller faster and reaches soon after a few updates, so the form cannot be kept.

To conclude for , under a small distortion, recovers to and recovers to the form . However, becomes smaller slowly with speed proportional to (quasi-equilibrium, not fixed point) unless (fixed point).

(A2). When and .

(A2.1). For the weight, under the assumption that the distortion is small, we have . For and , we have , so , and
formula
A.15
Therefore and recovers to for , and . Note that holds only for small distortion; otherwise, the dynamics goes to another state.

(A2.2). For the delay, we do not check whether it is of the form because these synapses do not contribute to the postsynaptic spike ().

Similar to (A2.1), under the assumption that the distortion is small, we have ; therefore,
formula
A.16
where .

Similar to (A1.2), becomes larger slowly with speed proportional to (quasi-equilibrium, not fixed point) unless (fixed point). The exact value is indeterminate, which depends on the ratio between delay-weight learning rates . But we can still give a range: is always larger than because , and upper bounded by .

(A3) When .

For the case of no spike for , recovers to if or keeps its initial value if ; the delays are kept at the initial value. In the following, we consider only the case with random noise spikes for .

(A3.1) For the weight, the spikes arrive randomly before or after the postsynaptic spike at . In this case,
formula
A.17

Therefore, and recovers to for .

(A3.2) For the delay, we have
formula
A.18
Meanwhile the variables satisfy the following conditions: , , , and . Therefore equation A.18 continues as
formula
A.19
where B is a variable on the same scale of and and is a variable on the same scale of and . Therefore, and recovers to for .

In conclusion, combining (A1, A2, and A3), under a small distortion, recover to , with T becoming smaller slowly with speed proportional to (quasi-equilibrium, not fixed point) unless (fixed point).

A.6.2  Part B of the Proof

All quasi-equilibrium states are of the form of .

We equivalently prove its contrapositive: if the weights and delays are not of the form of , it is not a quasi-equilibrium state.

First, there is no equilibrium state when . This is because the weight learning window is half negative and half positive with no zero point within its effective length. When the weight is larger, the latency is reduced, so the weight is further strengthened, which forms a positive feedback and finally leads to . The weakened weight also follows a positive feedback and leads to .

Next, for the delays, we show that is not a quasi-equilibrium state and converges to the form of .

(B1) For synapses in the set , their weights are maximized. The only reason is that the spikes arrive regularly before the postsynaptic spike: (assume negligible ). Similar to equation A.11, the relative difference between spike arrival times of any synapses j and k is diminishing to zero when , and , so dij cannot be kept and evolves to the form , where .

(B2) For synapses in the set , their weights are minimized. There are two possible reasons (assume negligible ): the spikes arrive regularly after the postsynaptic spike or arrive randomly.

In case of the spikes arriving regularly after the postsynaptic spike, , the synapse satisfies . Meanwhile dij becomes larger slowly due to and satisfies .

In case of the spikes arriving randomly, the synapse is not involved in the structural PNG: . Since the integral of is negative, dij cannot be kept and evolves to .

In conclusion, if the weights and delays are not of the form of , it is not a quasi-equilibrium state. Equivalently, all quasi-equilibrium states are of the form of , where . Since Q is not unique, the quasi-equilibrium states can have different Qs.

A.7  The Roles of Noise in the Network Dynamics

Noise in the neuron set P can make large deviations so that , leading to the switching behavior between multiple quasi-equilibrium states with different Qs (see Figure S1 in the online supplement). The switching behavior potentially leads to robust learning, since the readout neuron can switch to a different Q to avoid malfunctioned neurons.

Noise in the strengthened neuron set Q accelerates the speed of recovering to the quasi-equilibrium state, since larger ( for ) in equation A.5 makes M larger and thus smaller.

Noise outside P leads to for . This helps the readout neuron learn new PNGs when the previously learned one disappears, because with the short delay, spikes from a previously weakened synapse are more likely to arrive before the postsynaptic spike to get strengthened.

Acknowledgments

We thank Eugene M. Izhikevich from Brain Corporation and Jose C. Principe from the University of Florida for the initiative discussions on polychronous neuronal group. We thank Sui Guorong from the University of Shanghai for Science and Technology for the insightful discussions on the principles of neural networks.

References

Bakkum
,
D. J.
,
Chao
,
Z. C.
, &
Potter
,
S. M.
(
2008
).
Long-term activity-dependent plasticity of action potential propagation delay and amplitude in cortical networks
.
PLoS One
,
3
(
5
),
e2088
.
Bollmann
,
J. H.
, &
Sakmann
,
B.
(
2005
).
Control of synaptic strength and timing by the release-site CA2+ signal
.
Nature Neuroscience
,
8
(
4
),
426
434
.
Buonomano
,
D. V.
, &
Maass
,
W.
(
2009
).
State-dependent computations: Spatiotemporal processing in cortical networks
.
Nature Reviews Neuroscience
,
10
(
2
),
113
125
.
Buzsáki
,
G.
(
2010
).
Neural syntax: Cell assemblies, synapsembles, and readers
.
Neuron
,
68
(
3
),
362
385
.
Eurich
,
C. W.
,
Pawelzik
,
K.
,
Ernst
,
U.
,
Thiel
,
A.
,
Cowan
,
J. D.
, &
Milton
,
J. G.
(
2000
).
Delay adaptation in the nervous system
.
Neurocomputing
,
32
,
741
748
.
Gerstner
,
W.
,
Kistler
,
W. M.
,
Naud
,
R.
, &
Paninski
,
L.
(
2014
).
Neuronal dynamics: From single neurons to networks and models of cognition
.
Cambridge
:
Cambridge University Press
.
Gewaltig
,
M. O.
, &
Diesmann
,
M.
(
2007
).
NEST (NEural Simulation Tool)
.
Scholarpedia
,
2
(
4
),
1430
.
Gilson
,
M.
,
Masquelier
,
T.
, &
Hugues
,
E.
(
2011
).
STDP allows fast rate-modulated coding with Poisson-like spike trains
.
PLoS Computational Biology
,
7
(
10
),
e100223
.
Guise
,
M.
,
Knott
,
A.
, &
Benuskova
,
L.
(
2014
).
A Bayesian model of polychronicity
.
Neural Computation
,
26
(
9
),
2052
2073
.
Guyonneau
,
R.
,
VanRullen
,
R.
, &
Thorpe
,
S. J.
(
2005
).
Neurons tune to the earliest spikes through STDP
.
Neural Computation
,
17
(
4
),
859
879
.
Hüning
,
H.
,
Glünder
,
H.
, &
Palm
,
G.
(
1998
).
Synaptic delay learning in pulse-coupled neurons
.
Neural Computation
,
10
(
3
),
555
565
.
Izhikevich
,
E. M.
(
2006
).
Polychronization: Computation with spikes
.
Neural Computation
,
18
(
2
),
245
282
.
Lorteije
,
J. A.
, &
Borst
,
J. G. G.
(
2011
).
Contribution of the mouse calyx of held synapse to tone adaptation
.
European Journal of Neuroscience
,
33
(
2
),
251
258
.
Martinez
,
R.
, &
Paugam-Moisy
,
H.
(
2009
).
Algorithms for structural and dynamical polychronous groups detection
. In
Proceedings of Artificial Neural Networks–ICANN 2009
(pp.
75
84
).
New York
:
Springer
.
Masquelier
,
T.
,
Guyonneau
,
R.
, &
Thorpe
,
S. J.
(
2009
).
Competitive STDP-based spike pattern learning
.
Neural Computation
,
21
(
5
),
1259
1276
.
Masquelier
,
T.
,
Hugues
,
E.
,
Deco
,
G.
, &
Thorpe
,
S. J.
(
2009
).
Oscillations, phase-of-firing coding, and spike timing-dependent plasticity: An efficient learning scheme
.
Journal of Neuroscience
,
29
(
43
),
13484
13493
.
Senn
,
W.
,
Schneider
,
M.
, &
Ruf
,
B.
(
2002
).
Activity-dependent development of axonal and dendritic delays, or, why synaptic transmission should be unreliable
.
Neural Computation
,
14
(
3
),
583
619
.
Sun
,
H.
,
Yang
,
Y.
,
Sourina
,
O.
, &
Huang
,
G. B.
(
2015
).
Runtime detection of activated polychronous neuronal group towards its spatiotemporal analysis
. In
Proceedings of the International Joint Conference on Neural Networks
(pp.
1
8
).
Piscataway, NJ
:
IEEE
.
Sutherland
,
G. R.
, &
McNaughton
,
B.
(
2000
).
Memory trace reactivation in hippocampal and neocortical neuronal ensembles
.
Current Opinion in Neurobiology
,
10
(
2
),
180
186
.
Szatmáry
,
B.
, &
Izhikevich
,
E. M.
(
2010
).
Spike-timing theory of working memory
.
PLoS Comput Biol
,
6
(
8
),
e1000879
.
Tetzlaff
,
C.
,
Dasgupta
,
S.
,
Kulvicius
,
T.
, &
Wörgötter
,
F.
(
2015
).
The use of Hebbian cell assemblies for nonlinear computation
.
Scientific Reports
,
5
,
art. 12866
.
Tolnai
,
S.
,
Englitz
,
B.
,
Scholbach
,
J.
,
Jost
,
J.
, &
Rübsamen
,
R.
(
2009
).
Spike transmission delay at the calyx of held in vivo: Rate dependence, phenomenological modeling, and relevance for sound localization
.
Journal of Neurophysiology
,
102
(
2
),
1206
1217
.
Zenke
,
F.
,
Agnes
,
E. J.
, &
Gerstner
,
W.
(
2015
).
Diverse synaptic plasticity mechanisms orchestrated to form and retrieve memories in spiking neural networks
.
Nature Communications
,
6
,
art. 6922
.