Spiking-bursting activity allows the simultaneous propagation of different units of information throughout the network (multicoding). Therefore, different spatio-temporal spaces can be simultaneously used to globally encode and store information. In the network discussed in this paper, we consider two coexisting units of information in each neural signal: the bursting frequency and the neural fingerprints included within the burst. In the first dimension, the network must generate and coordinate spatiotemporal patterns of propagating transient bursting activity (rhythmic encoding modality). To achieve this, we impose two constraints (Wiedemann and Lüthi, 2003; Tabak et al., 2010): (i) predominance of excitatory synapses and (ii) a refractory period in each neuron following hyperexcitation. Information processing in the second dimension is based on the emission and recognition of specific neural signatures (Tristán et al., 2004; Carrillo-Medina and Latorre, 2015), i.e., information in this dimension propagates encoded in a spike-timing modality. An intra-unit contextualization mechanism drives the signature emission and recognition processes. This does not only allow us to illustrate a novel information processing strategy in the context of spiking neural networks, but also the dynamical richness that subcellular plasticity can provide to these networks.

#### 2.1.1. Neuron Spontaneous Dynamics

Many spiking models generate output bursts depending on the parameter settings and/or the input stimuli (e.g., the models by Hindmarsh and Rose, 1984, Komendantov and Kononenko, 1996, or Liu et al., 1998 that we have previously used to investigate the functional meaning of neural signatures). However, simulations show that the neural signatures in these models mainly depends on the network connectivity (Latorre et al., 2002) and, to our knowledge, none of the existing spiking models displays an adaptive fingerprint as required by our study. A possible alternative to this issue is using the mechanism described in Marin et al. (2014) to tune neuron busting models and produce neural signatures equivalent to those observed in living cells. However, the generation of realistic signatures is out of the scope of this proof-of-concept.

To describe the individual behavior of each unit, we define a stochastic model operating in a discrete event framework. The neuron activity is considered as a discrete variable and characterized in time by V(t), its "membrane potential." **Figure 1A** illustrates schematically the neuron spontaneous dynamics. Our model neuron integrates and processes the information received through its different input channels (synaptic integration), adapts its firing pattern to the incoming information (intra-unit plasticity), and generates a coherent signed output signal. During subthreshold activity, the spontaneous evolution of the neuron activity is determined by the probability p—the transit probability of the internal state per time step. When the membrane potential of a neuron ni reaches the firing threshold (TH), this generates a sequence of spikes (not a single spike). The temporal distribution of spikes within the response burst is given by a firing sequence composed of concatenating the signature (Si) and the preferred output spike pattern (Pi) of the neuron. Then, the stochastic dynamics of a single neuron depends on the temporal evolution of the neuron activity and whether it is under (subthreshold activity) or over (spiking-bursting activity) the firing threshold. Formally:

• During subthreshold activity (Vi(t) < TH):

$$V\_i(t+1) = \begin{cases} V\_i(t) + I\_{syn} + 1 & \text{with probability } p \\ V\_i(t) + I\_{syn} & \text{otherwise} \end{cases} \tag{1}$$

where Isyn is the synaptic input (Equation 3) and p the transit probability of the internal state per time step.

FIGURE 1 | (A) Schematic representation of the stochastic neuron model (see main text for details). S = {ISI1, ISI2} and P = {P1, P2} denote the neuron signature and the preferred output pattern, respectively. Note that the intraburst firing pattern is different in the first and the second burst. This is because, as (B) illustrates, the neuron recognizes a signature at the time step pointed by the arrow and intra-unit plasticity changes the neuron response in the spike-timing encoding modality. (B) Example of signature recognition. For each incoming spike, the local informational context keeps track of the corresponding input channel and spike timing (e.g., C1-100 means that at time step 100 a spikes arrived to the neuron through channel C1). This transient memory provides an intra-unit contextualization mechanism to the single neuron. For example, if the arrow in (A) corresponds to time step 136 and an input spike arrives through channel C1, the neuron can contextualize this spike and determine that the signature {8, 3} have been received four times in the recent history. If this value is greater than the learning threshold (Li ), the neuron recognizes this signature and, consequently, modulates its output firing pattern as illustrated in the second burst of (A). (C) Network topology. Each neuron is directly connected to its eight nearest neighbors with periodic boundary conditions. Then, neighbors of the white unit are the blue neurons.

• During the generation of the burst (Vi(t) >= TH):

$$V\_i(t+1) = \begin{cases} AP & \text{if } t = t\_1 + t\_n \\ TH + 1 & \text{if } t = t\_1 + t\_n + 1 \text{, } \forall n \neq N \\ 0 & \text{if } t = t\_1 + t\_N + 1 \\ V\_i(t) + 1 & \text{otherwise with probability } p \\ V\_i(t) & \text{otherwise} \end{cases} \tag{2}$$

where N is the number of spikes in the firing sequence (Si + Pi), tn denotes the timing of the nth spike in this sequence, (i.e., t1 corresponds to the initial timing of the burst) and being AP the peak membrane potential to generate a spike. Note, that during the burst generation synaptic input (Isyn) is not taken into account (cf. Equation 1 and 2). After generating a burst, neurons have a refractory period of RP time steps during which Vi(t) = 0. Then, subthreshold dynamics starts again.

#### 2.1.2. Synaptic Input

Synaptic input arrives to a neuron through two kind of input channels: connections with other neurons and an external channel to introduce external stimulation into the network. Each neuron in the network is connected to its eight nearest neighbors (**Figure 1C**) with periodic boundary conditions. As in every spiking neural network, neurons communicate with each other through the generation and propagation of spikes. Then, the interchange rule is defined by:

$$I\_{\rm sym} = \mathcal{g}\_{\rm e} \cdot pulse\_{\rm e} + \sum\_{j} \mathcal{g}\_{ji} \cdot pulse\_{j} \tag{3}$$

where ge defines the weight of the external stimulus, pulsee is 1 when an action potential is delivered through the external channel and 0 otherwise; and, similarly, gji is the weight of the connection between neurons nj and ni and pulsej is 1 when Vj(t − 1) = AP and 0 otherwise. Note that Equation 3 does not apply neither during the generation of a burst (Equation 2) nor during the refractory period, i.e., in these situations synaptic input is not considered.

It is important to highlight that in this paper we do not discuss synaptic learning (see Section 4). This implies that gji is constant for all the synapses and, consequently, the neighborhood of every neuron does not change.

#### 2.1.3. Intra-Unit Plasticity

Incorporating subcellular plasticity to a neuron model implies that a mechanism inside the cell allows tuning the neuron dynamics to incoming signals and/or to particular processing states. We consider here a history-dependent contextualization mechanism driving the spike-timing encoding modality. This intra-unit contextualization modulates the preferred output pattern as a function of previous incoming spike patterns.

As in the non-spiking signature neural network paradigm, to implement local contextualization, each individual neuron uses a transient memory, called local informational context. The local informational context keeps track of the information received during a time window of Mi time units, providing a historydependent contextualization mechanism to the single neuron processing. In the case of our spiking network, for each incoming spike, the neuron stores in its local context the joint information about the input channel and the spike timing (**Figure 1B**). In this way, different intra-unit plasticity rules can be defined to take into consideration the input spike timings. In particular, the following rule can be used to recognize specific neural signatures:

• when a spike arrives to a target unit, this checks whether the spike pattern received though the corresponding input channel appears in its local informational context so many times as a given learning threshold, Li . If so, the receptor recognizes this fingerprint, which implies that the preferred output pattern is overwritten with the recognized fingerprint.

**Figures 1A,B** illustrate how intra-unit plasticity tunes the output firing pattern in response to the fingerprint recognition. During the generation of the first burst in the time series, the neuron does not recognize any signature. Therefore, there is not a preferred output pattern and the burst only contains the signature of the neuron (S = {ISI1,ISI2}). At time step 136 (pointed by the arrow), an spike arrives through channel C1. The neuron can use its local informational context (**Figure 1B**) to contextualize this spike. In our case, this means to identify the incoming pattern through this channel (in this case {8, 3}) and to determine that this fingerprint has been received four times from time step 100. Then, assuming that the learning threshold is Li = 4, the neuron's preferred output spike pattern changes due to the recognition of the signature S ′ = {8, 3}. As a consequence, the intraburst firing pattern of the second burst in the time series varies to encode additional information (in the example, the sequence P = {P1 = 8, P2 = 3}). The neuron emits the new preferred output pattern until a new fingerprint is recognized or until the recognized fingerprint appears less than Li times in the local informational context (keep in mind that this is transient memory). Note, that intra-unit plasticity can be used to compute different aspect of the output signal as a function of the local contextualization, not only the spiking firing pattern. For instance, a particular cell could increase/decrease its level of activity or generate an output spike in response to specific incoming patterns independently of the synaptic weight.

During the input processing, channels are checked randomly in each iteration. In this way, when the target neuron recognizes multiple signatures in the same iteration, the last processed prevails over the others. Plasticity rule does not apply during the generation of a burst—i.e., once the neuron starts firing, the output spike pattern cannot change.

## 2.2. Analysis Methods

#### 2.2.1. Rhythmic Encoding Modality

To illustrate the spatio-temporal patterns generated in the bursting informational dimension, we generate activity movies representing the membrane potential evolving dynamics. In these movies, the evolution in time of the activity of a given unit (Vi(t)) is represented with a color scale. Regions with the same color have synchronous behavior. Red corresponds to neurons with a membrane potential over the firing threshold (Vi(t) > TH), i.e., they are generating a burst. Intermediate colors between blue and red, represent subthreshold activity. The cooler the color, the lower the level of activity.

Spatio-temporal patterns of spiking or spiking-bursting activity in one dimensional signals are usually detected and analyzed by means of spectral methods. However, in higher dimensions, the coefficients produced by the multidimensional Fourier transform are hard to interpret. On the other hand, wavelet-based techniques have proven to be useful tools for signal analysis (Stollnitz et al., 1996; Mallat, 1999). Unlike the Fourier transform coefficients, the wavelet transform coefficients are determined both by a resolution component and a time (or space) component and, therefore, they represent the resolution content at a given portion of the original signal. Thus, to quantitatively characterize the bursting rhythmic activity in our network, we perform a wavelet-based analysis.

In particular, we use the same discrete wavelet transform (DWT) analysis employed in Latorre et al. (2013a) to characterize the global network dynamics of a model of the inferior olive. The method consists in considering the spiking-bursting spatiotemporal patterns produced by the network as sequences of images evolving in time. As a first step in the characterization, a two-dimensional basis is generated by direct Cartesian product of the one-dimensional Haar basis (Stollnitz et al., 1996). Then, the two dimensional non-standard DWT is calculated for each frame of network activity. The idea behind this characterization method is that the number of wavelet coefficients in a given frame, C(t), provides an estimation of the complexity of the image corresponding to the spatio-temporal pattern at time t. A low number of coefficients means that the image is smooth or is composed of smooth components. In contrast, a high number of coefficients corresponds to complex images. In this way, the DWT analysis transforms the multidimensional spiking-bursting activity in the network, Vi(t), into a one dimensional signal, C(t). This signal provides an useful characterization of the bursting dynamics in which both the frequency and the spatial complexity can be discussed. From the frequency perspective, a simple visual inspection of the evolution of C(t) allows to detect the presence of different rhythmic patterns in the network. Furthermore, these rhythms can now be studied by means of the one dimensional Fourier transform. From the spatial complexity of the patterns, very high values of C(t) correspond to almost random behavior of every neuron, with no patterns present; intermediate high values indicate the presence of complex spatial structures in the patterns; while completely synchronized networks produce a small number of coefficients. Note, that C(t) ranges between 0 and the number of neurons in the network.

#### 2.2.2. Spike-Timing Encoding Modality

The spike-timing encoding is related to the spreading of specific intraburst spike patterns through the network and the synchronization mechanisms that allow a group of neurons to recognize and emit the same signature at a given moment (Tristán et al., 2004; Carrillo-Medina and Latorre, 2015). To graphically illustrate the dynamic spatial organization of the spike patterns within the network, we generate activity movies representing the fingerprint-based evolving dynamics (e.g., see **Figure 5**). Each point in the 50 × 50 square represents with a color code the neural signature recognized by a given neuron within the network at a given moment. In this manner, neurons with the same color recognize the same signature. White color identifies the units that do not recognize any fingerprint.

To quantitatively analyze this encoding strategy, we compute the evolution of the number of neurons that recognize and emit each individual signature per time unit. This measure provides an estimation over time of the level of activity in the network related to each signature.

### 3. RESULTS

We have conducted experiments in which multiple datasets are presented to regular networks with different parameters. Independently of the network size and the number of neighbors per neuron, it is possible to find a broad range of synaptic weights and neuron parameters allowing the network to simultaneously encode information in the rhythmic and the spike-timing modality. However, the emerging phenomena that we describe here can be more easily illustrated in autonomous networks with a low level of bursting activity, since in these cases, the spatio-temporal activity in the different dimensions arises due to external stimulation. In autonomous networks, i.e., networks not receiving external input, the level of bursting activity depends on the transit probability of the internal state (p), the firing threshold (TH), and the duration of the refactory period (RP). These parameters modulate the ratio of bursts produced by an isolated neuron. The greater the value of the stochastic probability p, the higher the mean bursting frequency. Similarly, the bursting frequency also grows with low values of TH and RP.

Thus, in the following sections, we focus on neurons where p = 0.05, TH = 50, RP = 50, and AP = 200 (units are dimensionless). Note, that AP—the peak membrane potential to generate a spike—has no influence on information processing, the only requirement is being greater than TH. We discuss results of square-shaped networks of 50×50 of such units, with periodic boundary conditions and where each unit is connected through an excitatory synapse (gji = 1) to its eight nearest neighbors as shown **Figure 1C**. External stimuli consist of tonic spiking signals at a given frequency introduced into a randomly chosen cell during a give time period. The neural signature of every neuron has six spikes, with all the ISIs in the range 2–12 (dimensionless). These signatures are randomly generated and assigned at the beginning of the simulation. The rest of parameters are specified in the corresponding experiment description. Vi(0) is chosen randomly between 0 and 40 a.u. for all neurons in the network.

#### 3.1. Rhythmic Encoding Modality

The degree of synchrony among the membrane potential of the neurons constituting the network characterizes the global spiking-bursting activity in the network. For fixed values of p and TH, the degree of synchrony varies as a function of the synaptic transmission strength among neurons (i.e., gji in Equation 3). For small values, each neuron fires nearly independently. As synaptic weights grow, the degree of synchrony increases because the generation of a burst in a given unit sequentially propagates to its neighbors and so on (**Figure 2**). The higher synchrony occurs in networks with combinations of firing thresholds and excitatory synapses that allow a target neuron to reach the firing threshold between frames of 33 a.u.

when it receives a burst (gji · #spikes\_in\_burst ≥ TH). However, as we mention above, here we are interested in autonomous with a low level of bursting activity.

Depending on the synaptic parameters, burst propagation provides autonomous networks the ability to generate welldefined spatio-temporal patterns in the form of propagating wave fronts of transient spiking-bursting activity. Note that local contextualization modulates intraburst firing patterns, but it has not any influence on burst timings. To illustrate these spatiotemporal patterns, we generate activity movies representing the membrane potential evolving dynamics (see Section 2.2.1 for details). As representative example of the spontaneous collective bursting rhythms generated by the network, bottom panel in **Figure 2** displays snapshots of the activity movie of the network shown in the top panel.

The spontaneous generation of transient spatio-temporal patterns of spiking or spiking-bursting activity is a feature with relevant functional implications observed in different living neural media. However, we are interested in the network response to stimuli. Therefore, from the encoding perspective, the most interesting feature of the network, appearing even in networks with a small synaptic transmission among neurons, is its ability to develop dynamical patterns of spiking-bursting activity in response to data onset. These patterns allow the network to encode information using the frequency of different bursting rhythms induced by stimuli. To illustrate how the network of **Figure 2** encodes a single input using this spatio-temporal space, **Figure 3A** shows snapshots of its collective spiking-bursting dynamics when a unit in the left-top corner receives an external tonic spiking signal. When the stimulus is introduced into the neuron, its firing frequency increases. Then, the spiking-bursting

FIGURE 3 | (A) Snapshots of an activity movie illustrating the spiking-bursting spatio-temporal patterns generated by the network of Figure 2 (Mi = 400 a.u. and Li = 4) when a tonic input with a period of 100 time units between consecutive spikes is introduced into a single neuron (arrow in the first frame points to the neuron that receives incoming stimulus). Sequences develop in time from left to right and from top to bottom. The time interval between frames is 33 a.u. The stimulated unit increases its bursting frequency due to the external stimulation, and this generates new spatio-temporal patterns of transient spiking-bursting activity from this unit. (B) Characterization with the DWT coefficients of the activity of the network in the top panel: first without stimuli (snapshots in Figure 2B belong to this period), then when the selected neuron receives the incoming data (grayed area identifies the period while the input is active), and finally without any input again. (C) Normalized power spectra of the wavelet analysis for the three periods. Left: without stimuli. Middle: during the stimulation. Right: when the stimulation is over and after the reverberation period. Power spectra are calculated using time series of 500,000 time units. The DWT analysis demonstrates that the global network dynamics changes when data are introduced into the network. It also shows that the network is a dynamical working memory of spiking-bursting rhythms, since the network dynamics generated in response to data onset reverberates after the input is retired.

activity originated in the stimulated neuron propagates to the surrounding units because of excitation. Thus, this neuron becomes the origin of a new rhythm that coexists with those generated spontaneously by the network (if any).

The DWT analysis (Section 2.2.1) corroborates the rhythm encoding in the network transient spiking-bursting dynamics. Changes in the collective spiking-bursting dynamics in response to data onset are reflected in a change in the evolution of the DWT coefficients whose shape characterizes the spiking-bursting activity of the network. As an example, **Figure 3B** illustrates how the collective dynamics of the network in **Figure 3A** changes when data are introduced into the stimulated unit. Initially, no data is present and the network spontaneously generates spatiotemporal patterns as the ones shown in the bottom panel of **Figure 2**. In this situation, the DWT coefficients oscillate with a nearly homogeneous frequency capturing the spontaneous spiking-bursting rhythm. The spontaneous rhythm frequency depends on the stochastic probability p and can be estimated by means of the Fourier transform of the wavelet analysis of the network activity. For instance, in the network of **Figures 2**, **3**, the spontaneous rhythm frequency is around 6.2·10−3 (see frequency peak in **Figure 3C**, left). On the other hand, the oscillation of the DWT coefficients between a high and an intermediate value indicates, respectively, the nearly independent neuron behavior during subthreshold activity and a high transient synchronization in the network during the spreading of the spiking-bursting wave fronts. Then, the external stimulus is introduced into the network during a given time interval (grayed area). At this point, the network collective dynamics stepwise changes. A first remarkable change in the evolution of DWT coefficients is observed in the oscillation amplitude. Now, the DWT coefficients tend to oscillate around two high values. This change points out the complex spatial structure of the new emerging dynamics. Not obtaining low or intermediate values in the DWT analysis during the stimulation period indicates that, in this network, the propagation of the wave fronts originated in the stimulated unit does not imply a complete transient synchronization in the whole ensemble. Another relevant change in the DWT coefficients during the stimulation is a frequency increase (cf. left and middle power spectra in **Figure 3C**), pointing out that the rhythm evoked by the stimulus prevails over the spontaneous rhythm (6.2 · 10−3 vs. 10 · 10−3 a.u. −1 ). The frequency of the spiking-bursting rhythms evoked by external stimulation depends on the frequency of the input, since the stimulated neuron follows the stimulus. These changes indicate that the network has encoded the incoming information in a characteristic spiking-bursting rhythm. Finally, no input is present again and the network recovers the spiking-bursting autonomous activity (cf. **Figure 3C**, right). The DWT analysis indicates that the stimuli-evoked rhythms can reverberate for long periods after data onset. This implies that the network behaves as a working memory in the spiking-bursting spatiotemporal space. For each network configuration, the mean reverberation period of the rhythms encoding different inputs is nearly the same, i.e., the memory capability of the network in this information dimension is independent of the data and only depends on the synaptic parameters.

The emerging collective dynamics analysis in networks that receive multiple tonic stimuli with different frequencies indicate that spatio-temporal patterns of spiking-bursting activity allow the network to encode information using several coexisting and coordinated rhythms. Top panel in **Figure 4** displays an example of the complex spatial organization of the patterns generated by a network receiving 10 different inputs. The snapshots clearly show the increased complexity of the patterns, since, now, the network organizes clusters of neurons oscillating at different frequencies (cf. top panel in **Figure 3**). Each of the unit receiving external data becomes the source of a rhythm that propagates through the network competing with the rhythms encoding other inputs.

FIGURE 4 | Figure equivalent to Figure 3 but when the network (Mi = 400 a.u. and Li = 4) receives 10 inputs. (A) Sequences develop in time from left to right and from top to bottom. The time interval between frames is 33 a.u. In response to data onset, the network starts generating 10 different coexisting rhythms encoding incoming information. (B) The different spiking-bursting rhythms encoded within the network are captured by the DWT analysis. While external inputs are present, the oscillation frequency of the DWT coefficients is not homogeneous (see inset), which reveals the coexistence of the different rhythms. Inset shows the normalized power spectrum of the wavelet analysis of a time series of 500,000 time units while the 10 external stimuli are present. The number of coefficients increase (cf. Figure 3B) denotes the increase in the spatial complexity of the patterns.

As we show above, while an input is active, the corresponding rhythm survives in the network. Therefore, when more than one stimulus is present, the competition among the inputevoked spiking-bursting rhythms is a winnerless competition. Note that there is no inhibition in the network nor subcellular plasticity rules limiting the spiking-bursting activity. Winnerless competition allows the encoding of multiple coexisting spikingbursting rhythms. This competition dynamics is captured by the DWT analysis (bottom panel in **Figure 4**). When multiple data are introduced into the network, the number of DWT coefficients remains high with a non-homogeneous oscillation frequency. This reveals the complex spatial structure of the patterns and, on the other hand, the coexistence of multiple spiking-bursting rhythms within the network.

We have previously shown that the spiking-bursting rhythms evoked by a single stimulus reverberate for a while when the stimulation is over. The reverberation period drastically increases when the network receives multiple stimuli. The greater the number of external inputs, the greater the number of sources of spiking-bursting activity. This translates into a higher spikingbursting activity in the network and explains the increasing reverberation period. Depending on the synaptic strength and the value of p, in this situation, the network even becomes a long-term memory of spiking-bursting rhythms. We would like to emphasize that the rhythms that survive for longer periods in

short-term memories or the ones that persistently reverberate in long-term memories are not always the higher frequency stimulievoked rhythms nor the rhythms encoding the last data presented to the network.

## 3.2. Spike-Timing Encoding Modality

One of the major characteristics of the proposed network is the intra-unit contextualization of input signals, responsible of the spike-timing encoding modality. In this section, we study the complex collective dynamics induced by this intraunit information processing strategy. These emerging collective dynamics can give us important clues about the underlying computational properties of the network.

As expected, the ability of an individual unit to recognize specific fingerprints varies as a function of the intra-unit parameters shaping local contextualization, i.e., the maximum size of the local informational context (Mi) and the fingerprint learning threshold (Li). Depending on the value of these parameters, specific intraburst firing patterns can propagate through autonomous networks. However, the more interesting phenomena from the information processing viewpoint are related to the mechanisms that allow the network to generate and organize spatio-temporal patterns in response to data onset. Therefore, we focus our attention on networks in which the signature recognition does not occur without external stimuli. When these networks receive incoming data, they aid the study of the information encoding in the fingerprint-based spatiotemporal space by analyzing how the signatures of the stimulated units propagate throughout the network.

Again, we first address the analysis of networks receiving a single stimulus. When a neuron receives an external tonic input, this unit increases its bursting frequency (see Section 3.1). This increase can make the neighbor units recognize the neural signature of the stimulated neuron and propagate the corresponding intraburst firing pattern. In this situation, new intriguing collective dynamics arise in the network. To illustrate the dynamic spatial organization of the neural signatures traveling through the network, we generate activity movies representing the fingerprint-based evolving dynamics (see Section 2.2.2 for details). These activity movies point out that the network generates in this dimension well-defined transient patterns of activity in response to data onset. The emerging spatio-temporal patterns are related to the spatial organization and clusterization of the signatures traveling through the network. To give insight into the generation and propagation of these complex spatio-temporal structures, **Figure 5** shows snapshots of the activity movies of two representative networks in which the same unit receives an input. Note, that the only signature traveling through the network corresponds to the stimulated unit. If we consider that at a given moment two neurons that recognize the same signature belong to the same cluster; we can study the specific properties of the dynamic organization of the patterns by calculating the clustering coefficient and the average shortest path between neurons belonging to the same cluster. This analysis indicates that the fingerprint-based spatio-temporal patterns are initially originated in the stimulated unit (see initial frames in the

FIGURE 5 | Snapshots of two representative activity movies illustrating the fingerprint-based encoding mechanism. (A) Mi = 500 a.u. and Li = 5. (B) Mi = 400 a.u. and Li = 4. Sequences develop in time from left to right and from top to bottom. The time interval between frames is 1000 a.u. Note that the propagation of the fingerprint-based spatio-temporal patterns is slower than the corresponding spiking-bursting rhythms (cf. bottom panel and Figure 3). The color code identifies neurons recognizing the same signature, being white color used for neurons that do not recognize any signature. The first frame in each sequence indicates that, in the absence of stimuli, neural signatures do not propagate in these networks. When the external stimulus is introduced into a neuron located in the left-top corner (second frame in both panels), new collective dynamics emerge and the network organizes transient spatio-temporal patterns of activity related to the propagation of the signature of the stimulated unit (blue regions). Note that this is the only signature that travels throughout the network. These localized patterns of activity encode the who of incoming data.

sequences of **Figure 5**). Then, depending on the parameters Mi and Li , they can propagate locally or globally as transient wave fronts; or as localized clusters with a fixed spatial organization that occasionally become the source of new transient patterns. The generation of localized transient patterns of activity in the fingerprint spatio-temporal space suggests a collective coding strategy based on the emission and recognition of specific neural fingerprints. This mechanism allows the network to encode information regarding the origin of incoming data (input source) in a distributed network form.

The information encoded in the spike-timing modality and the encoded in the rhythmic modality coexist in the network. A relevant property observed in the simulations is that a neural fingerprint does not necessarily travel over the propagating wave fronts encoding the corresponding spiking-bursting rhythm. The spreading velocity of the fingerprint-based spatio-temporal patterns is always slower than the corresponding spikingbursting spatio-temporal patterns velocity (cf. time interval between frames in **Figures 3**, **5**; 33 vs. 1000 a.u). Likewise, the spatial organization of the patterns in the different spatiotemporal spaces is not correlated. If we consider that at a given moment two neurons over the firing threshold belong to the same cluster, we can calculate the clustering coefficient and the average shortest path for the spiking-bursting patterns and compare the self-organizing properties of the patterns encoded in both information dimension. This analysis points out that the spiking-bursting patterns always consist of propagating transient wave fronts from the stimulated unit traveling through the whole network. Meanwhile, the fingerprint-based patterns can also be originated in the stimulated unit, but they can propagate locally or globally as transient wave fronts or remain bounded in specific regions of the network.

A simple way to characterize the fingerprint-based dynamics is computing the number of neurons that recognize and emit a given firing pattern. This allows us to identify the signatures encoded in the network. **Figure 6A** depicts the characteristic evolution of the level of activity related to the fingerprint of the data source in three representative networks receiving the same single input during three different stimulation periods. This figure corroborates the results derived from the snapshots shown in **Figure 5**. When the stimulation begins, the signature of the stimulated unit starts propagating through the network. The number of neurons recognizing this fingerprint grows until reaching a stationary level that depends on the value of Mi and Li . Then, the network dynamics consists of a fluctuation around the steady level (e.g., see blue traces in **Figure 6A**). This dynamic is kept while the stimulation is sustained. When the stimulation ends, the stimulus-evoked activity does not immediately disappear from the network (cf. red and green traces in **Figure 6A**). This is an interesting result that demonstrates that intra-unit contextualization can be a mechanism to implement intrinsic memory in the network, giving rise to both shortterm and long-term memories. In short-term memories (bottom and middle panel in **Figure 6A**), the stimuli-evoked dynamics reverberate for a while. This reverberation effect constitutes a mechanism providing the network the ability of acting as a dynamical working memory that transiently stores incoming data. In contrast, in long-term memory networks (top panel in **Figure 6A**), the information survives in the network in a permanent manner (maybe until a new input is received).

The collective dynamics in the fingerprint-based dimension is mainly driven by the intra-unit parameters Mi and Li . On the one hand, reducing the size of the local informational context of every neuron (Mi) decreases the number of neurons that recognize a given firing pattern. On the other hand, decreasing the learning threshold (Li) facilitates the recognition of the propagating fingerprints and, therefore, the level of activity in the network grows. The trade-off among the effect of these parameters determines if the network encodes information in the spike-timing modality and the mode of behavior in this dimension. To illustrate this, **Figure 6B** depicts a phase diagram locating the different behaviors in the space of intra-unit parameters.

With the experiments described so far, we investigate the ability of the proposed network to encode and process a single stimulus using an information processing strategy driven by local contextualization. If we repeat the same experiments but now introducing multiple inputs simultaneously, we observe that the presence of multiple stimuli makes the network generate coexisting transient spatio-temporal patterns of activity encoding the origin of the different inputs (**Figure 7**). These experiments reveal additional relevant computational properties that subcellular plasticity can provide to spiking neural networks. When multiple intraburst firing patterns spread through the network, a competition dynamics arises between them. A simple visual inspection of the snapshots shown in **Figure 7** reveals that the self-organizing properties of the patterns drastically change depending on the intra-unit parameters shaping the intra-unit plasticity rules. These define different modes of competition among the spreading fingerprints. This competition affects the global level of activity of each signature in the network and determines the spatial organization of the patterns. The competition dynamics among the different intraburst firing

FIGURE 6 | (A) Evolution of the mean number of neurons that recognize and emit the fingerprint of a unit receiving the same data in three different networks during three different periods. Each trace is calculated as the average of 10 experiments with different random seeds and location of the stimulated unit. These plots characterize the stimuli-evoked fingerprint-based dynamics. Top panel: Mi = 500 and Li = 5. Middle panel: Mi = 400 and Li = 4. Bottom panel: Mi = 350 and Li = 4. Units are dimensionless. In red traces, the stimulation period corresponds to the red region. In green traces, to the green region. And in blue traces, data are continuously present. In this spatio-temporal space, the network may act as a long-term memory (top panel) or as a short-term memory (middle and bottom panels) depending on the value of Mi and Li , i.e., the parameters associated to intra-unit contextualization. (B) Phase diagram illustrating the relationship between Mi and Li in networks where p = 0.05, TH = 50, RP = 50, and AP = 200 (units are dimensionless).

As result of the competition, only the patterns associated to a limited group of data (the winners) propagate. Note that the different competition regimes arise

patterns determines the coherence and coordination of the

depending on the values Mi and Li which shape the intra-unit contextualization mechanism.

coexisting patterns. We would like to highlight that the competition regimes observed in the activity movies arise in the absence of inhibitory connections, which hints at intra-unit contextualization as an effective mechanism to restrict the activity in networks without inhibition. Note that each neuron can only transmit one recognized firing pattern per burst. This limitation produces somehow a local competition among the patterns received by the neuron where only the "winner" is transmitted. This local competition is the basis of the global competition in the whole network.

The different dynamical modes observed in the activity movies are better characterized by the evolution of the number of neurons that recognize and emit each signature (**Figure 8**). Regardless the number of active inputs, the type of competition depends on the value of the parameters Mi and Li and may vary from a winnerless (WLC) to a winner-take-all (WTA) competition. In WLC networks, none of the signatures becomes a "winner," and therefore, none of them persistently prevails over the others. Depending on the intra-unit parameters, the network can display different winnerless regimes. **Figure 8A** illustrates a winnerless competition in which the level of activity related to every fingerprint is similar and remains fluctuating nearly a stationary level. This defines a collective dynamics where several coherent spatio-temporal patterns coexist within the network encoding simultaneously a great amount of data (e.g., in **Figure 8A** all the inputs introduced into the network). **Figures 8B,C** show winnerless regimes with alternating periods where some fingerprint has a higher level of activity.

An interesting phenomenon observed with some network settings is that some regions within the network specialize in the emission of firing patterns encoding the origin of different stimuli although they do not receive any external input. This phenomenon occurs without any kind of supervised synaptic nor intra-cellular learning, i.e., it is a self-organizing property of the network. These emitter areas are usually related to winnerless competitions where the prevailing fingerprints change

accordingly to the patterns originated in these areas (**Figure 9**). Conversely, when a winner-take-all competition occurs, only the signature or signatures that win the competition propagates through the network (e.g., see **Figure 8D**). In the WTA network shown in **Figure 7D**, all the neurons tend to recognize and emit simultaneously the prevailing fingerprint. However, depending on Mi and Li , this can also spread as evolving transient patterns equivalent to the shown in **Figure 7C** when the dark green input prevails over the others. Note that, in some sense, the winnerless competitions displayed in **Figures 7B,C** consist of sequences of transient winner-take-all competitions.

The reverberating spatio-temporal patterns encoding the origin of incoming data continue competing even when they are not sustained by an active input. Short-term memory networks have a limited ability to retain previously stored data when new information is introduced into the network. In these cases, the reverberation period drops as compared to networks receiving a single input, and the stored data are almost instantaneously forgotten, i.e., the corresponding patterns disappear because the patterns encoding the last incoming stimulus win the competition. However, in long-term memory networks, coexisting coherent spatio-temporal patterns related to multiple fingerprints can be observed even when the corresponding input is not active.

## 4. DISCUSSION

The present work introduces a spiking neural network that makes use of multicoding strategies for information propagation and subcellular plasticity to locally contextualize or discriminate data received by a unit. Furthermore, each neuron in the network has a neural signature that allows its unequivocal identification by the rest of the cells. This network is an encoder and generator of spatio-temporal patterns that take advantage of the multiple simultaneous encoding modalities present in the network to transform dynamic inputs into different spatio-temporal spaces, and organize and coordinate coexisting patterns of transient activity in response to data onset.

The discussed experiments are aimed at analyzing the emerging collective dynamics in two information dimensions.

and from top to bottom. When neurons in the pointed area start generating patterns encoding a given input, the collective behavior changes accordingly to these patterns and the corresponding fingerprint prevails over the others. Note that the existence of these emitter areas is a self-organizing property of the network.

On one hand, a spiking-bursting spatio-temporal space, where information processing is driven by synaptic transmission.  The specific properties of the dynamic organization of the patterns are different in each information dimension, so that, the life cycle of the information encoded in both encoding schemes is independent. When multiple patterns in the same dimension coexist in the network, a competition emerges between them. We show that various forms of competition can arise without inhibitory connections in the network. Depending on the parameters shaping simple intra-unit plasticity rules, the competition regime may vary from a winnerless (i.e., the network stores multiple data simultaneously) to a winner-takeall competition (i.e., one datum or a group of them prevails over the others). The stimuli-evoked spatio-temporal patterns and the corresponding competing dynamics can survive for long periods after data onset. This reverberation effect allows the network to memorize incoming data. This can display shortterm or long-term memory capabilities in the different spatiotemporal spaces. When the network behaves as a short-term memory, the spatio-temporal patterns encoding incoming data in the corresponding scheme transiently reverberate after the stimulation ending. Conversely, in long-term memories, the stimulus leads the network to a new stable state and the patterns persistently survive. The memory ability of the network in each dimension varies as a function of the synaptic and/or intraunit parameters. Therefore, different simultaneous processing strategies can be implemented within the network.

These results illustrate the dynamical richness and large flexibility of the proposed network to encode and process information in different spatio-temporal spaces. We argue that plasticity mechanisms inside individual cells and multicoding strategies can provide additional computational properties to spiking neural networks, which could enhance their capacity and performance. In particular, local contextualization mechanisms allow individual neurons to process the multiple simultaneous codes in their input signals selectively or globally in order to completely decide or weight the decision about their output in the different encoding schemes. This information processing provides a framework to model complex highdimensional processes that can be applied to different real-world computational problems. The ideas relating multicoding with local information discrimination have a direct application in problems that benefit from multifunctionality and parallelism. These are desirable features for many technical applications of ANNs, representing a potential advantage when processing large amounts of data or multiple decision-making criteria must be developed, for instance, in multiobjective optimization problems (Saini and Saraswat, 2013; Wang et al., 2014) or in control systems [e.g., multifunctional prosthesis controllers that must quickly detect and classify multiple characteristic simultaneous myoelectric signals (Saridis and Gootee, 1982; Hudgins et al., 1993; Karlik et al., 2003; Li et al., 2010)]. Another straightforward application of these concepts is in problems where a global task is solved by means of solving independent partial tasks. An example is the wide scope of multidimensional sorting problems, specifically when the order in a particular dimension can be independent of the order in other dimensions, or when there is no global sorting criteria in any dimension. Non-spiking signatures neural networks have been successfully applied to this type of problems (Latorre et al., 2011). Areas of application for multidimensional sorting are scheduling, planning and optimization, between others (Catoni, 1998; Aref and Kamel, 2000). On the other hand, the different dynamical modes observed in the network are relevant in the context of multiple technical applications. Winnerless competition is usually associated to sequential information processing (Seliger et al., 2003; Rabinovich et al., 2006a; Arena et al., 2009; Kiebel et al., 2009; Latorre et al., 2013b), which has a wide application in many artificial intelligent systems in tasks such as inference, planning, reasoning, natural language processing, and others (Sun and Giles, 2001; Wörgötter and Porr, 2005). Similarly, pattern recognition in different spiking ANNs is based on winner-take-all dynamics (Bohte et al., 2002a; Gütig and Sompolinsky, 2006; Schmuker et al., 2014).

In this paper, we have imposed some constraints and assumptions in order to facilitate the presentation of our results. Results obtained with larger regular networks (up to 1000 × 1000); higher levels of bursting activity; and different number and/or distribution of spikes in the neural signatures are equivalent to the results presented in Section 3. In experiments with signatures with an arbitrary number of spikes, new interesting fingerprint-based dynamics emerges in the network and results are not exactly the same. In these simulations, not only the fingerprints belonging to a neuron propagate, but also specific firing sequences built with combinations of these signatures propagate throughout the network. In some sense, these networks do not only encode information regarding the input source, but they also generate new information. It is also important to note that, for simplicity, we only consider two encoding schemes in the network. However, bursting activity allows easily including additional units of information (e.g., the burst duration or the number of spikes in the burst). In this line, and regarding a selective processing of input messages, experimental evidence indicates that some neural systems exhibit functional or behavioral neural signatures representing different states or associated to the task performed at a given moment (Klausberger et al., 2003; Somogyi and Klausberger, 2005; Kaping et al., 2011). The concept of neural fingerprint that underlies the strategy of the discussed network can be extended to consider the emission and recognition of multiple fingerprints with a different meaning within the same signal. In this situation, subcellular plasticity in the form of intraunit information contextualization mechanisms would allow individual neurons to perform a distinct processing of incoming signals, for example, as a function of specific emitters and/or functional states.

Although not addressed in this paper, subcellular plasticity and multicoding mechanisms for information processing can be combined with the features that underlie information processing in the existing spiking neural network paradigms. In this line, for example, plenty of work has been done on synaptic plasticity in spiking neural networks, since modifications of the synaptic connections are traditionally considered the physiological basis of learning in the nervous system. These works are mostly related to unsupervised synaptic learning methods, such as Spike-Timing Dependent Plasticity (STDP) (Song et al., 2000; Bohte et al., 2002b; Kube et al., 2008; Meftah et al., 2010), with an increasing interest into supervised synaptic learning (Bohte et al., 2002a; Belatreche et al., 2007; Yu et al., 2013). The combination of learning rules including not only the modification of the synaptic weights, but also the parameters that affect the local discrimination of input signals can greatly contribute to enhance the spiking ANNs' computational power. In this vein, our results can be of particular interest in the context of the generation and recognition of spatio-temporal information. Different spiking neural networks have been proposed to process, classify, and store spatio-temporal patterns (Laje and Buonomano, 2013; Yu et al., 2013). We speculate that incorporating multicoding strategies and different types of subcellular plasticity to other successful spiking ANN paradigms can potentially allow these networks to process, classify and store more complex data. For example, a highly relevant application of the referred spiking networks is the analysis of EEG spatio-temporal data.

#### REFERENCES

To consider a multicoding mechanism that incorporates the neural fingerprint-based dimension to these networks could permit an analysis of coexisting brain rhythms from multiple simultaneous perspectives. In particular, the fingerprint-based spatio-temporal patterns could facilitate the analysis of the propagation trajectories and the identification of possible information sources and sinks in different cognitive processes.

Because of their functional similarity to biological neurons, spiking neural networks have been extensively used by the computational neuroscience community as a powerful tool for studying neural information processing (e.g., see Izhikevich, 2003; Deco et al., 2008; Izhikevich and Edelman, 2008). Results obtained with our simple model could also be relevant from this perspective Information storage in the nervous system has been typically studied considering the adaptation of the synaptic connection strengths (e.g., see Zipser et al., 1993). Our simulations suggest that mechanisms inside individual cells modulating their intrinsic dynamics could also be an effective mechanism to implement intrinsic memory, both in shortand long-term memory networks. On the other hand, many biological neural systems (including many areas of the human brain) continuously receive a great amount of inputs from many different sources and, nevertheless, they exhibit a low level of activity and only respond to specific inputs (Shoham et al., 2006; Sato et al., 2007; O'Connor et al., 2010; Barth and Poulet, 2012). We hypothesize that neural dynamics based on the propagation of specific neural fingerprints and a contextualization mechanisms like the one studied here could explain why these system are so sparsely active. Target neurons would only fire when they recognize a characteristic firing pattern in their incoming stimuli; while signal not recognized would be simply ignored. Obviously, to test this hypothesis more realistic spiking models for the activity of the neurons must be developed.

#### AUTHOR CONTRIBUTIONS

RL conceived and designed the study. JC and RL conducted the experiments. JC and RL analyzed the data. JC and RL wrote, read, and approved the manuscript.

#### FUNDING

This work was supported by UAM-Banco Santander (CEAL-AL/2015-16) and MINECO/FEDER DPI2015-65833-P.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Carrillo-Medina and Latorre. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Mechanisms of Winner-Take-All and Group Selection in Neuronal Spiking Networks

Yanqing Chen\*

The Neurosciences Institute, La Jolla, CA, USA

A major function of central nervous systems is to discriminate different categories or types of sensory input. Neuronal networks accomplish such tasks by learning different sensory maps at several stages of neural hierarchy, such that different neurons fire selectively to reflect different internal or external patterns and states. The exact mechanisms of such map formation processes in the brain are not completely understood. Here we study the mechanism by which a simple recurrent/reentrant neuronal network accomplish group selection and discrimination to different inputs in order to generate sensory maps. We describe the conditions and mechanism of transition from a rhythmic epileptic state (in which all neurons fire synchronized and indiscriminately to any input) to a winner-take-all state in which only a subset of neurons fire for a specific input. We prove an analytic condition under which a stable bump solution and a winner-take-all state can emerge from the local recurrent excitation-inhibition interactions in a three-layer spiking network with distinct excitatory and inhibitory populations, and demonstrate the importance of surround inhibitory connection topology on the stability of dynamic patterns in spiking neural network.

#### Edited by:

Sander Bohte, Centrum Wiskunde & Informatica, Netherlands

#### Reviewed by:

Da-Hui Wang, Beijing Normal University, China Victor de Lafuente, Universidad Nacional Autónoma de México, Mexico

#### \*Correspondence:

Yanqing Chen<EMAIL_ADDRESS>
Received: 12 November 2016 Accepted: 20 March 2017 Published: 21 April 2017

#### Citation:

Chen Y (2017) Mechanisms of Winner-Take-All and Group Selection in Neuronal Spiking Networks. Front. Comput. Neurosci. 11:20. doi: 10.3389/fncom.2017.00020 Keywords: neuronal spiking network, phase transition, learning and memory, Winner-take-all (WTA), neural computation, Robotics

## 1. INTRODUCTION

Facing with vast amount of multi-sensory information, Central Nervous System (CNS) seems to process only a small subset of those inputs at any given time, no matter whether they come from external or internal sources.  At the level of neuronal networks, a network in which all neurons respond the same to all stimuli would convey no information about the stimulus. In order to be useful, neurons must come to respond differentially to variety of incoming signals. Many neural models and theories have been proposed to account for such ability. Winner-Take-All (WTA) network is one of such proposed mechanisms for developing feature selectivity through competition in simple recurrent networks, and it has received much attention on both theoretical and experimental grounds. The primary theoretical justification is the ability of such networks to explain how the maps, which are ubiquitous in the cerebral cortex, can arise (Kohonen, 1982; Goodhill, 2007). WTA networks can also explain how a network can come to make useful distinctions between its inputs. WTA networks coupled with synaptic learning rules and homoestatic plasticity can explain how this takes place in a self-organized fashion from an initially undifferentiated state. Finally, WTA models are often employed at the behavioral level in

theoretical models of higher-level cognitive phenomenon such as action-selection, attention (Itti and Koch, 2001; Walther and Koch, 2001) and decision making (Wang, 2002; Furman and Wang, 2008).

Another mechanism proposed for feature selectivity is the phenomenon of spatially localized bumps in neuronal networks (Somers et al., 1995; Laing and Chow, 2001; Wei et al., 2012). If we view multiple neurons within a bump as mulitiple-winners of excitatory and inhibitory competition, bump activity in spiking networks can be treated as a soft WTA or k-Winner-Take-All phenomenon (see Maass, 2000 for their definition). In this paper we use Winner-Take-All (WTA) and "bump activity" inter-changeably to describe the same stable group activity that arises from inter-connected excitatory-inhibitory neuronal networks. On a more general level, both bump activity and WTA phenomenon can be viewed as a type of pattern formation process in networks of excitatory and inhibitory neurons (for example, patterns of stable grid in Wilson and Cowan, 1973; Ermentrout and Cowan, 1979), and an example of activity dependent neuronal group selection process (Edelman, 1987).

Population rate-based WTA models have been extensively studied and are well understood (Dayan and Abbot, 2001). But, the connections between rate models to the real biological neural systems are not direct, because they are different from the real nervous systems whose neurons are spiking. So it is necessary to study the networks of spiking neurons, such that the biological interpretation of spike models can be more directly linked to real nervous systems. Modeling and understanding spiking networks is not simple because spiking neurons are highly nonlinear and their action potentials are discrete. As a result, it is always more difficult to obtain analytical solutions for spiking firing properties than rate models.

Analysis has shown conductance-based spiking models can be approximated by simple rate models under certain conditions (such as in an asynchronous state in Shriki et al., 2003). This approach has been applied to the study of hyper-column in a spiking model of visual cortex (Shriki et al., 2003). The orientation selectivity in their study, is modeled as the appearance of a unimodal "bump"-like spiking activity in a ringconnected spiking network, similar to an earlier study (Laing and Chow, 2001). Both approaches applied approximations from the rate models and used Fourier analysis to calculate the conditions for the appearances of bump activity. Recent work specifically studied recurrent spiking WTA networks, which are closer to real biological systems than previous rate models (Rutishauser and Douglas, 2009; Rutishauser et al., 2011). Even though these newer network models can receive spike input and generate spike output, their network structures are still very simplified. For example, excitatory and inhibitory neurons are modeled into one single population (Laing and Chow, 2001), and inhibitory population are reduced into one unit (Rutishauser et al., 2011), or removed altogether and modeled as direct inhibitory connections among excitatory neurons (Oster et al., 2009).

In a recent report we presented a robust and more biologically-realistic WTA network structure with distinct excitatory and inhibitory populations with arbitrary number of units (Chen et al., 2013). This WTA network has been implemented into a robot that accomplished a sequence learning and mental rotation task (McKinstry et al., 2016). In our spiking models each neuron type has very detailed biological parameters to model different neuronal transmitters and receptor types similar to previous work (Izhikevich and Edelman, 2008). We showed that surround inhibition and longer time constants from NMDA and GABAb conductances are sufficient to achieve stable "bump" spiking activity in a selected winner neuronal group while all the other neurons are inhibited and quiet. However, detailed biological properties, such as STSP (short-term synaptic plasticity), NMDA voltage gating etc., prevented a formal analytical analysis of the whole model. Also, it is not clear any of those biological details or a specific type of synaptic connections are crucial for the emergence of bump activity.

To identify the most important mechanistic factors for the spiking WTA networks, here we study a simplified spiking network after some biological details are removed. For example, based upon what we have noticed previously, turning off STSP, NMDA voltage-gating and excitatory-to-excitatory connections does not change the overall properties of WTA phenomenon. On the other hand, we preserve some important biological features such as the four different synaptic connections and conductance types (AMPA, NMDA, GABAa, and GABAb), because we found that these four individual conductance types contribute to different aspects of the "bump" stability. By examining functions of these individual conductances and the topologies of excitatoryinhibitory connectivity, we provide a detailed analysis of the conditions on which a stable bump activity can emerge from this recurrent spiking network. Our analysis thus provide a mechanistic analysis on how a neuronal group selection process can occur in an activity dependent manner in neural systems.

## 2. METHODS

## 2.1. Network Structure

Here we analyze a basic 3-layer spiking neuronal network with different neuron types with realistic biological parameters. The first layer of excitatory neurons (IN – input cells) takes input signals (e.g., arbitrary analog patterns) and translates them into spiking activity. The input signal we considere here in this paper is a type of unstructured random currents evenly distributed within a certain range and injected into the 100 input neurons (IN). IN cells are randomly connected to the next excitatory layer (E) with initial weights evenly distributed between 0 and a maximal value. The random input currents and random connections to the excitatory layer we analyzed here provide a baseline condition in which we test how the recurrent/reentrant connectivity between excitatory and inhibitory neurons by themselves can accomplish winner-take-all competition to random but unstructured input patterns (without obvious firingrate differences among input neurons) and without synaptic modifications. The successful WTA network structure then can be trained to discriminate more complex and structured patterns through spike-timing dependent learning rules such as STDP. Such learning process will modify the synapses between these two excitatory types so that a selected E and I neurons (the WTA group) will respond to preferred input patterns more quickly for practical applications. We have demonstrated these in previous reports (Chen et al., 2013; McKinstry et al., 2016) where the same WTA network structures were implemented in a humanoid robot to process real world complex visual inputs, to learn visual-motor association and sequencing, and to accomplish a "mental rotation" and delayed-match-to-sample task.

We also implemented the above network using adaptive exponential spiking models and obtained similar results. For simplicity the analysis below uses the Izhikevich model (Izhikevich and Edelman, 2008), and excitatory (E) and inhibitory (I) neurons use the same parameters in the following equation:

$$C\dot{\nu} = k(\nu - \nu\_r)(\nu - \nu\_t) - \mu - I\_{sym} \tag{1a}$$

$$
\dot{u} = a \{ b(\nu - \nu\_r) - u \} \tag{1b}
$$

Parameters in these equations are the same as explained before (Izhikevich and Edelman, 2008). That is, v is the membrane voltage in millivolts (mV), C is the membrane capacitance, vr is the neuron's resting potential, vt refers to its threshold potential, u represents the recovery variable defined as the difference of all inward and outward voltage-gate currents. Isyn is the synaptic current (in pA) originated from spike input from other neurons. a and b are different constants. When the membrane voltage reaches a threshold, i.e., v > vpeak, the model is said to generate a spike, and two variables in Equations (1a, 1b) are reset according to v ← c and u ← u + d while c and d are parameters for different cell type.

We use a simplified synaptic current form with four basic conductances from AMPA, NMDA, GABAA, and GABAB channels. For simplification, voltage-gating of NMDA channel is reduced to a constant factor. This is done through calculating an average number for the voltage-gating term for the NMDA conductance (i.e., [(v + 80)/60]2 /[1 + ((v + 80)/60)2 ] on Page 11 of the Supplementary Information in (Izhikevich and Edelman, 2008)) for the normal range of voltages: v = [−60, 60] , and the result is equivalent to a voltage-independent NMDA channel with smaller gain factor than AMPA channels (see Appendix for the individual conductance gain factors we used). So synaptic current Isyn is composed of four different current types originated from those four conductances multiplied with the voltage differences between their individual reversal potentials:

$$I\_{\rm syn} = \mathcal{g}\_{\rm AMPA} \left( \nu - 0 \right) + \mathcal{g}\_{\rm NMDA} \left( \nu - 0 \right) + \mathcal{g}\_{\rm GABA\_A} \left( \nu + 70 \right)$$

$$+ \mathcal{g}\_{\rm GABA\_B} \left( \nu + 90 \right). \tag{2}$$

Shown in the above equation, reversal potentials of AMPA and NMDA channels are 0 and reversal potentials for GABAA and GABAB channels are −70 and −90 mV respectively.

As described before, each conductance has exponential decay with different time constants in millisecond (ms):

$$
\dot{\mathbf{g}} = -\mathbf{g}/\mathbf{r},\tag{3}
$$

while τ = 5, 150, 6, and 150 for the AMPA, NMDA, GABAA, and GABAB channels respectively.

To simplify the analysis, there are equal numbers (400 in all the subsequent analysis) of excitatory (E) and inhibitory (I) neurons in our basic network model in **Figure 1**, although their numbers can be in any ratio. In fact, in our previous published full models (Chen et al., 2013; McKinstry et al., 2016) the ratio of E and I neurons were set at 4:1 to more closely resemble the real cortex. We also explored different types of connection topologies in the connections from excitatory to inhibitory neurons (E to I), the reentrant inhibition from basket cells to pyramidal neurons (I to E) and the inhibitory connections within basket cells themselves (I to I). In our study, Inhibitory to Excitatory and Inhibitory to Inhibitory connections are kept the same topological type and total weights are kept equal. Throughout the simulation the total connection weights to each neuron are normalized to be a constant for each connection type. The total weights for each connection type (E to I and I to E) are two parameters we explored systematically. As a first step, we firstly only consider one type of inhibitory conductance (GABAA) to obtain analytical solutions for the conditions of Winner-Take-All state. GABAB conductances are added after an analytical solution is found, a comparison of the transition plots can be found in the Appendix.

peaks with a flat gap (zero value connectivity) in between. We call this specific network connectivity as surround inhibition type for the one dimensional case and it is a simplified version of the two-dimensional Central-Annual-Surround

(CAS) type of topology we described before (Chen et al., 2013).

## 3. RESULTS

To classify different types of spike dynamics for the surround inhibition network in **Figure 1**, for each neuron, we record the number of spikes between 2 and 3 s after the simulation had reached steady state without synaptic plasticity (STDP off). We then characterize the behavior of the network by the maximum number of spikes generated by any excitatory neuron. **Figure 2E** plots this maximal firing rate for every combination of the E to I weights vs. the I to E weights. The analysis is repeated and plotted in **Figure 2F** for the inhibitory population.

**Figure 2** shows different types of dynamic firing patterns in the 2-dimensional parameter space. When only one type of connection weight (excitatory or inhibitory) is high but the other weight is low, either excitatory or inhibitory neurons are in a quasi-random/rhythmic state in which one group of neurons fires in high Gamma frequency range (>40 Hz, see **Figures 2A,G**). When both connection weights are relatively high (see **Figure 2C**), both excitatory and inhibitory neurons have high maximal firing rates where excitatory neurons have a maximal firing rate larger than 35 Hz and inhibitory neurons have a maximal firing rate of larger than 100 Hz. If we look at the corresponding spike raster plot in **Figure 2C**, only a subset of excitatory and inhibitory neurons maintain such high firing rates while majority of other neurons are silent. We call this Winner-Take-All (WTA) state in which only a small subset of neuronal groups persistently fire high frequency and keep the rest of neurons from firing using surround inhibition. The region of the

FIGURE 2 | Classification of different dynamic spiking patterns in the surround inhibition network we defined in Figure 1 and the phase diagram for the transition between different firing states. (A–C,G–I) Are raster plots which show all spikes within a half second interval for each neuron in the network. (E,F) are maximal firing rates of excitatory and inhibitory neurons in the 2-d parameter space (total excitatory weight in y-axis and total inhibitory weight in the x-axis). Red Points in (D) are transition curves constructed from (E) where maximal firing rate of excitatory neurons changed from below to above 40 Hz. (A–C,G–I) Are example spiking patterns under their specific parameter combinations which are marked on the 2-d plot in (D). Subplot (C) represents a winner-take-all (WTA) state where only a small group of excitatory and its corresponding inhibitory neurons fire persistently while others are silent. Black curves on the middle row–(D,E,F) are the same analytic transition condition for the WTA region based upon the analysis described in the Appendix (Equation A.14). Basically it is a curve where SE · SI equals to a constant defined by Equation (A.14), and it delineates the Winner-Take-All region (i.e., max firing rate large than 40 Hz for only a selected group of excitatory neurons) in the parameter space very well for both the excitatory and inhibitory neurons.

parameter space with such WTA states is delineated by the right curve in **Figure 2D** where the maximal firing rate of excitatory neurons increased to greater than 35 Hz from lower firing rates in the middle region (from the blue area in **Figure 2E** transition to the red area on the top right), and roughly corresponds to a similar increase of maximal firing rate to above 100 Hz for inhibitory neurons in **Figure 2F**.

Subplots **Figures 2A,B,H,I** all belong to an intermediate region in the parameter space in **Figure 2B** between two curves where maximal firing rates for both excitatory and inhibitory neurons are relatively low. Within this parameter range, excitatory and inhibitory neurons are either quasi-synchronized (**Figure 2A**) or precisely synchronized and firing rhythmically (**Figure 2H**), or exhibit moving bump activity (**Figure 2I**) or as combinations of rhythmic and moving bump activity. In all these cases, single excitatory neuron cannot maintain a stable high gamma frequency spiking activity unless connection weights are changed, moving to the WTA region on the top-right of the second curve in **Figure 2D**. **Figure 2D** thus provides a phase diagram for the neuronal network defined in **Figure 1**.

Notice that this maximal firing rate is not the neuron's instantaneous firing rate, but is the total number of spikes within a 1 s window. This definition is useful to discriminate a stable high firing rate neuron vs. a neuron firing a short burst less than 1 s and then becoming quiet (especially for stable vs. traveling activity, see **Figure 2C** vs. **Figure 2I**).

**Figure 3** summarizes patterns of spike dynamics with different connection topologies. Compared to the surround inhibition type analyzed above, all the other connection types do not support a Winner-Take-All state manifested as stable bump activity shown in **Figure 2C**. This is because under those connection types, excitatory and inhibitory neurons cannot maintain high maximal firing rates when both excitatory and inhibitory weights are high and did not have a red area on the upper-right region shown in **Figures 2E,F**. The most common firing patterns for those connection types are quasi-rhythmic firings in the 10–20 Hz range for excitatory neurons resembling an epileptic state while some short burst of unstable bump activity in inhibitory neurons. Our results suggest that, among different types of connectivity topologies we analyzed, only

surround inhibition can generate a stable bump spiking activity and maintain a WTA state.

## 3.1. Mechanism of Winner-Take-All Neuronal Group Selection and Emergence of Bump Activity

Our above analysis suggests that surround inhibitory topology supports emergence of bump activity. To explore the mechanism of WTA and which neuronal properties are essential for such behavior, we applied the same analysis as in **Figure 2** to neuronal network in **Figure 1** when Short-Term-Synaptic-Plasticity (STSP) or NMDA voltage-gating is on, or change excitatory and inhibitory neurons' parameters to different type. In all cases, a similar WTA region was found for every conditions, even though the transition curves that delineate the emergence of stable bump activity are shifted to different positions in the parameter space (see results in Chen et al., 2013). We also analyzed the same neuronal network with a different set of individual spiking models, i.e., the adaptive exponential models and found the similar WTA region as long as the topology of the inhibitory connections are surround type. These analyses suggest that detailed neuronal properties such as exact models of the spiking neuron, STSP or NMDA gating etc., are likely not fundamental for the existence of stable bump activity, but the type of connectivity topology (i.e., surround inhibition) is more important for such behavior.

Both the Izhikevich neural model and the adaptive exponential model we used are conductance based with models of inhibitory and excitatory currents of different time scales. So we suspect that different time constants of NMDA, AMPA, GABAa, and GABAb channel conductance might play some role for the emergence of bump activity. To demonstrate this, **Figure 4** shows the time evolutions of AMPA, NMDA, and GABAa currents along with the spiking activity in the simplified network in **Figure 1** starting from a zero conductance initial condition. It demonstrates the detailed transition from a rhythmic synchronized firing state into a stable bump activity. Looking at detailed dynamic changes of the individual excitatory and inhibitory currents should shed light on how the transition is occurred.

From **Figure 4** we can see, differences in time constants will determine how fast a specific channel conductance returns back to zero after a burst of spiking activity. For excitatory neurons specifically, because of its short time constants (6 ms), AMPA current fluctuates around a similar level with large variances. NMDA currents, on the other hand, are accumulating to higher levels because of longer time constant (150 ms) even though both currents are generated by the same spiking input from the input neurons. Similar phenomenon can be seen for inhibitory neurons. When those neurons fire rhythmically before about 300 millisecond, AMPA conductance jumps to high level (from 7 to 9 nS) after each spike then drops down to zero very fast (red curve in **Figure 4D**), while NMDA conductance only drops a small amount each cycle and overall level still increases to much higher value (red curve in **Figure 4E**). Initially inhibitory neurons fire after excitatory neurons in each rhythmic cycle and they synchronize to each other with a time delay. If excitatory to inhibitory weights (E to I) are larger than a certain value, such that NMDA currents for inhibitory neurons increase faster than excitatory neurons (**Figure 4E**), the delay between inhibitory and excitatory neurons diminishes and GABAa currents become effective within the same cycle to inhibit other neurons. As a result some inhibitory and excitatory neurons stop firing in the rhythmic cycle, eventually lead to a winner group that persistently fires and shuts off their surrounding neighbors. Notice that in the simplified model in **Figure 4**, GABAb (with longer time constant of 150 ms) currents are omitted and set to zero, which lead to a moving bump activity for this specific parameter set. If GABAb conductance is restored to the original level as in the full model, bump activity becomes stable. It implies that time constant of GABAa and GABAb channel conductance is related to the stability of the bump activity.

As a summary, we think the combinations of long and short time constants from excitatory and inhibitory conductance plus the surround inhibitory connectivity support a mechanism for emergence of bump activity and winner-take-all phenomenon in this basic spiking neuronal network. This neuronal group selection mechanism provides a basis for modeling learning and map-formation process for sensory motor integration and other higher cognitive processes.

## 4. ANALYTICAL ANALYSIS OF THE TRANSITION CURVE FOR WTA PHENOMENON

## 4.1. Differentiation in Inhibitory Conductances Lead to Spiking Activity Pattern Transition and Neuronal Group Selection

To identify the mechanism of bump activity in spiking networks, based upon the transition plots shown in **Figure 4**, we look at voltages, conductances of all 400 excitatory neurons at one specific time point (t = 980 ms. in **Figure 4**). **Figure 5** shows AMPA, NMDA, and GABAa conductances along with voltages for all those excitatory neurons at this specific time point. From the network structure defined in **Figure 1**, we know that the excitatory conductances (AMPA and NMDA) are determined by excitatory synapses from the input layer (because we omitted self-excitation), where inhibitory conductances (GABAa and GABAb) are only determined (triggered) by spikes from inhibitory neurons onto those excitatory neurons. Because of the uniform random connection from input layer, AMPA and NMDA conductances are around the same level and undifferentiated for all excitatory neurons. From **Figure 5**, voltages are above threshold and neurons fire only at locations where GABAa conductances below a certain value. So in order for a bump to emerge and a subgroup of neurons selected to be active, GABAa or GABAb conductances have to be differentiated, i.e., they have to be small for some neurons and remain high for all other neurons. We suspect this condition can only be met by surround type of inhibitory connectivity. It is obvious that lowest GABA conductances lead to highest firing rate for excitatory

neurons. Since we have local feedforward excitation to inhibitory neurons in our network, the bump area in inhibitory neurons with highest firing rate should also have lowest inhibitory conductances. This difference in GABA conductances is true for both excitatory and inhibitory neurons because GABA conductances are determined by the same inhibitory spikes. This suggests that in order for a bump to emerge, local inhibition to the nearest neighbors should be lower than inhibition to neurons outside of the bump. Notice the three other inhibition topology in **Figure 3** all have peak (or flat) inhibition locally, so even if a neuronal group emerge with highest firing rate, the strong local inhibition will force their firing rates to decrease, and let the other sub-threshold neurons to fire. So this is likely the reason why we did not obtain stable bump activity using those inhibition connectivities. On the other hand, surround inhibition

winner-take-all state early and then bias the transition in excitatory cells.

type defined in **Figure 1** might be the most simple form of inhibition topology that could let a bump emerge and stabilize. Below analysis will further prove this point.

In contrast with models using negative weights to represent inhibitory connections, our spiking models' synaptic weights and excitatory/inhibitory conductances are all positive (which obviously is more biologically realistic). From Equation (2) we can see, it is only because of the differences in reversal potentials between excitatory and inhibitory channels, the current generated by excitatory and inhibitory conductances could have different signs (excitatory current coming into the neuron and inhibitory current coming out of neuron). In order for a neuron to fire, synaptic current has to be below a threshold Isyn < −Ith where Ith is about 100pA (Isyn has to be negative for a neuron to fire because it was defined as an outward 60

25

current in Equation 1a). Firing threshold Ith can be found from calculating F-I curve for the specific spiking neuron model we used. Figure A.3 (in the Appendix) plots the firing rates vs. amount of injected current (equivalent to −Isyn) for the Izhikevich neuron model result from numerical simulations. It shows that neurons start to fire when absolute value of the injected current is above 100 pA and then increase their firing rate approximately linearly until above 100 Hz. we can use this information to simplify the spiking activity into a rate model. As indicated on the last paragraph, AMPA and NMDA conductances are approximately uniform for excitatory neurons and they can not contribute to the differences in firing rates, so in order for the excitatory population to fire differentially, the difference between highest and lowest GABA conductances for individual neurons has to be larger than a certain value. This value can be estimated using Equation (2). If min(gGABAA ) is 0, for a resting potential of vr = −60mV, GABAA conductance has to be larger than the following value so injected synaptic current −Isyn will be below the firing threshold Ith:

$$\lg\_{\text{GABA}\_A} > \left(-I\_{l\mathfrak{h}} + 60 \* (\lg\_{\text{AMPA},E} + \lg\_{\text{NMDA},E})\right) / 10. \tag{4}$$

In **Figure 5**, gAMPA,E + gNMDA,E for excitatory neurons is around 4nS (Appendix will show how this value can be estimated analytically), so gGABAA has to be larger than 14 nS to keep subthreshold neurons from firing. This number is consistent with the result plotted on **Figure 4**, **5** that neurons fire and form a bump area where GABA conductances are below 14 nS and areas with GABA conductance larger than 14 nS are completely quiet.

If we consider both the GABAA and GABAB conductances based upon Equation (2) and using the same idea as above, conditions for inhibitory conductances will be the following for the winner-take-all state:

$$10 \cdot \text{g}\_{\text{G}ABA} + 30 \cdot \text{g}\_{\text{G}ABA} > \left( -I\_{th} + 60 \ast (\text{g}\_{\text{A}\text{MPa},E} + \text{g}\_{\text{N}\text{MDA},E}) \right) \tag{5}$$

Equations (4, 5) can be used further to identify the exact condition for the WTA state and to locate the transition curve in **Figure 2**. Using two cosine bumps as surround inhibitory connection weights, Appendix gives the analytic form of GABA distribution of a bump solution for neurons connected one dimensionally and uses it to obtain analytic conditions for the WTA state in the parameter space (see Equations A.13, A.14). Such analytic conditions are expressed as formulas combining single neuron property and the conductance parameters (such as time constants, gain factors for different inhibitory and excitatory conductances). Based upon these formulas we can locate the Winner-Take-All and bump activity in the parameter space fairly precisely (see black curves in **Figure 2** and the white curve in Figure A.4 in the Appendix), thus provide a mechanistic explanation for the emergence of winner-take-all state and stationary bump activity in this 3-layer spiking network we analyzed here.

### 4.2. Origin of Traveling Wave and Instability of Bump Activity – Driven by AMPA Conductances

Parameters in **Figure 4** are located very near the transition curve in the parameter space (see **Figure 2**), so the bump is not spatially stable and moves across different neurons. To identify the origin of such instability, we selectively set AMPA gain of excitatory or inhibitory neurons to 0 in order to see their effects on the bump stability. This is equivalent to selectively block AMPA conductances in either excitatory or inhibitory neurons in real biological neural systems. We found that if AMPA conductances in inhibitory neurons are set to 0 but not in excitatory neurons, bumps become more or less stable. On the other hand, when AMPA conductances are blocked and set to 0 in excitatory neurons but not in inhibitory neurons, we can have a moving bump with a constant spatial speed (see **Figure 6**). In fact, we can estimate the moving speed of bumps based upon the parameters we defined in **Figure 2**. So we believe the source of the bump instability is from the AMPA conductances in inhibitory neurons. Previous studies have associated stabilities of bump activity with dynamic synapses (Fung et al., 2009). Notice that AMPA conductances have much shorter time constants than others thus more associated with faster synapses (similar to GABAa), so in this sense there might be a connection and some agreement between our observation on **Figure 6** and dynamic synapses analyzed by Fung et al. (2009).

#### 5. CONCLUSION AND DISCUSSION

In this paper we derive global properties of spiking neuronal networks related to bump activity and Winner-Take-All state mainly through analysis of the dynamics of excitatory and

inhibitory conductances. To achieve the analysis of collective behavior, individual spiking properties are approximated by its firing rate property such as the conductance/firing rate curve (Figure A.3 in the Appendix) or F-I curves. In this regard, detailed properties of individual spiking model might not be crucial for global activity such as the emergence of bump and WTA state. For example, if we use different parameter sets for excitatory and inhibitory neurons such as changing the inhibitory neurons to be basket cell type, we can still found the WTA region in the parameter space in **Figure 2D**, but the exact location of the transition curve is shifted to a different place because basket neurons have different conductance/firing rate curve and different Ith, gth, k values in Equation (A.13). This could explain the transition curve in our full model with more detailed biological properties has the same functional form, but in different location in the parameter space (see Chen et al., 2013) because it included more detailed single neuronal spiking properties such as NMDA voltage gating and STSP etc. In fact, we used adaptive exponential spiking model to substitute the Izhikevich neuron models and obtained similar phase plot and transition curve for the bump activity and the WTA state.

We suggest that all conductance-based spiking models with distinct excitatory and inhibitory populations could have the similar collective Winner-Take-All behavior as analyzed here. Detailed spiking model properties such as F-I curve and firing threshold (Ith and gth) would determine the exact location of the transition curve in **Figure 2**. Global connectivity topology and different time constants (dynamics) of excitatory and inhibitory conductances are likely to be the determinant of system-wide spiking activity patterns.

#### 5.1. Importance of the Inhibitory Topology

The most important feature of our winner-take-all network is its surround inhibition topology. The reason we chose two sine function peaks as surround inhibitory connection weights is just because of its mathematical convenience, since convolutions of sine/cosine functions are much easier to solve than other types of functions. In fact, connection topologies using Gaussian peaks or torus (for two-dimensional neuronal arrays) were used in our previous model (Chen et al., 2013) and similar stable WTA results were obtained. We believe using other type of function for inhibitory connection would also work, as long as there is a low inhibitory weight locally. Comparing four different connection topologies from **Figures 2**, **3**, the reason why only surroundinhibition supports stable bump activities is because its maximal inhibitory connection weight is not to the nearest neighbors, but to slightly distant neurons. This gap of zero inhibitory weight can be very small (e.g., w down to 0 and equivalent to a no-self-inhibition case), and we can still find solutions for stable bump or bumps (in fact, w determines how many bumps can emerge and we will have a 2-bump solution when w is close to 0, see Appendix and Figures A.1, A.2). So as long as there is a local valley of inhibitory weight, stationary bumps could emerge because only decreasing inhibition could allow a bump to sustain.

Mechanistically it appears that the most important requirement for a bump solution is the stable differentiation in inhibitory (GABAA or GABAB) conductance distributions across the neuronal population. That is, for some neuronal groups, GABA conductances should be low to allow bumps to emerge and for the other neurons, they need to be high enough to keep the rest of neuronal population from firing spikes. As long as this condition is met, more detailed biological properties such as local self-excitation, short-term synaptic plasticity (STSP), voltage-gating of NMDA channels etc. can be added to the model without destroying the overall bump stability.

## 5.2. Why Traditional Center-Surround Topology Might Not Lead to Stable Bumps in Models with Distinct Excitatory and Inhibitory Populations

Previous rate-based population models (Dayan and Abbot, 2001) had most often used center-surround type of connection topology as shown in the middle row of **Figure 3** (Narrow excitation/Wide inhibition). Similarly, many spiking models with excitatory/inhibitory conductances on the same units used the same topology (Laing and Chow, 2001). By simple subtraction, narrow excitation and wide inhibition can lead to a "Mexican Hat" type of effective connectivity which supports winner-takeall in previous firing rate models. But, as we see from the analysis above, in biologically more realistic spiking models with distinct excitatory and inhibitory neuronal populations, multiple types of conductances cannot cancel each other easily because they are generated by precise spike timing and have different time constants. The "classical" center-surround topology can not guarantee a stable "Mexican Hat" type of net connectivity because sensitive spike timing differences between different neurons prevent easy subtraction of excitatory and inhibitory weights at every time point. In fact, as shown in **Figure 4**, the emergence of winner-take-all in spiking networks is a direct

#### REFERENCES

result of precise spike-timing–the coincide of excitatory and inhibitory population firing spikes lead to a sub-population of inhibitory neurons fire earlier than the rest of populations which then let them suppress and shut off the other neurons in the network (see **Figure 4G**). This is the reason that we believe a surround-type of inhibitory topology is essential for a stable spiking WTA network because it can support the emergence of a winner-group without shutting off themselves too early.

In summary, WTA network analyzed here demonstrates how variability and randomness in spiking time of individual neurons can lead to global pattern changes and phase transition in collective neuronal groups. Analytic solutions for the phase transition curve provided in this paper will help to increase our understandings of different functional roles of excitatory and inhibitory neural connections on the emergence and stability of firing patterns in the brain.

### AUTHOR CONTRIBUTIONS

YC developed the original idea of this paper, performed computational modeling and mathematical derivation for the analytical solution and wrote the paper.

#### ACKNOWLEDGMENTS

This work was supported in part by DARPA through ONR Grant N00014-08-1-0728 and by AFRL Cooperative Agreement FA8750-11-2-0255 to Neurosciences Research Foundation and by a grant from The G. Harold & Leila Y. Mathers Charitable Foundation. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon. The views, opinions, and/or findings contained in this article are those of the authors and should not be interpreted as representing the official views or policies, either expressed or implied, of the Defense Advanced Research Projects Agency, the Air Force Research Laboratory, the Department of Defense, or the U.S. Government.

## APPENDIX

The Appendix for this article can be found online at: http:// journal.frontiersin.org/article/10.3389/fncom.2017.00020/full# supplementary-material

**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Role of Architectural and Learning Constraints in Neural Network Models: A Case Study on Visual Space Coding

#### Alberto Testolin1 \*, Michele De Filippo De Grazia1 and Marco Zorzi 1, 2 \*

*1 Department of General Psychology and Padova Neuroscience Center, University of Padova, Padova, Italy, 2 San Camillo Hospital IRCCS, Venice, Italy*

#### Edited by:

*Marcel van Gerven, Radboud University Nijmegen, Netherlands*

#### Reviewed by:

*Michael W. Spratling, King's College London, UK Kandan Ramakrishnan, University of Amsterdam, Netherlands*

#### \*Correspondence:

*Alberto Testolin<EMAIL_ADDRESS>Marco Zorzi<EMAIL_ADDRESS>
Received: *30 November 2016* Accepted: *27 February 2017* Published: *21 March 2017*

#### Citation:

*Testolin A, De Filippo De Grazia M and Zorzi M (2017) The Role of Architectural and Learning Constraints in Neural Network Models: A Case Study on Visual Space Coding. Front. Comput. Neurosci. 11:13. doi: 10.3389/fncom.2017.00013* The recent "deep learning revolution" in artificial neural networks had strong impact and widespread deployment for engineering applications, but the use of deep learning for neurocomputational modeling has been so far limited. In this article we argue that unsupervised deep learning represents an important step forward for improving neurocomputational models of perception and cognition, because it emphasizes the role of generative learning as opposed to discriminative (supervised) learning. As a case study, we present a series of simulations investigating the emergence of neural coding of visual space for sensorimotor transformations. We compare different network architectures commonly used as building blocks for unsupervised deep learning by systematically testing the type of receptive fields and gain modulation developed by the hidden neurons. In particular, we compare Restricted Boltzmann Machines (RBMs), which are stochastic, generative networks with bidirectional connections trained using contrastive divergence, with autoencoders, which are deterministic networks trained using error backpropagation. For both learning architectures we also explore the role of sparse coding, which has been identified as a fundamental principle of neural computation. The unsupervised models are then compared with supervised, feed-forward networks that learn an explicit mapping between different spatial reference frames. Our simulations show that both architectural and learning constraints strongly influenced the emergent coding of visual space in terms of distribution of tuning functions at the level of single neurons. Unsupervised models, and particularly RBMs, were found to more closely adhere to neurophysiological data from single-cell recordings in the primate parietal cortex. These results provide new insights into how basic properties of artificial neural networks might be relevant for modeling neural information processing in biological systems.

Keywords: connectionist modeling, unsupervised deep learning, restricted Boltzmann machines, autoencoders, sparseness, space coding, gain modulation, sensorimotor transformations

## INTRODUCTION

Artificial neural network models aim at explaining human cognition and behavior in terms of the emergent consequences of a large number of simple, subcognitive processes (McClelland et al., 2010). Within this framework, the pattern seen in overt behavior (macroscopic dynamics of the system) reflects the coordinated operations of simple biophysical mechanisms (microscopic dynamics of the system), such as the propagation of activation and inhibition among elementary processing units. Though this general tenet is shared by all connectionist models, there is large variability in processing architectures and learning algorithms, which turns into varying degrees of psychological and biological realism (e.g., Thorpe and Imbert, 1989; O'Reilly, 1998). When the aim is to investigate high-level cognitive functions, simplification is essential (McClelland, 2009) and the underlying processing mechanisms do not need to faithfully implement the neuronal circuits supposed to carry out such functions in the brain. However, modelers should strive to consider biological plausibility if this can bridge different levels of description (Testolin and Zorzi, 2016).

Recent theoretical and technical progress in artificial neural networks has significantly expanded the range of tasks that can be solved by machine intelligence. In particular, the advent of powerful parallel computing architectures based on Graphic Processing Units (GPUs), coupled with the availability of "big data," has allowed to create and train large-scale, hierarchical neural networks known as deep neural networks (LeCun et al., 2015, for review). These powerful learning systems achieve impressive performance in many challenging cognitive tasks, such as visual object recognition (Krizhevsky et al., 2012), speech processing (Mohamed et al., 2012) and natural language understanding (Collobert et al., 2011). However, while the impact of deep learning for engineering applications is undisputed, its relevance for modeling neural information processing in biological systems still needs to be fully evaluated (for seminal attempts, see Stoianov and Zorzi, 2012; Khaligh-Razavi and Kriegeskorte, 2014; Güçlü and van Gerven, 2015).

One critical aspect of most deep learning systems is the reliance on a feed-forward architecture trained with error backpropagation (Rumelhart et al., 1986), which has been repeatedly shown to yield state-of-the-art performance in a variety of problems (LeCun et al., 2015). However, the assumptions that learning is largely discriminative (e.g., classification or function learning) and that an external teaching signal is always available at each learning event (i.e., all training data is "labeled") are clearly implausible from both a cognitive and a biological perspective (Zorzi et al., 2013; Cox and Dean, 2014). Reinforcement learning is a valuable alternative and it has already shown promising results when combined with deep learning (Mnih et al., 2015; Silver et al., 2016), but there is a broad range of situations where learning seems to be fully unsupervised and its only objective is that of discovering the latent structure of the input data in order to build rich, internal representations of the environment (Hinton and Sejnowski, 1999). We argue that more realistic neurocognitive models should therefore also exploit unsupervised forms of deep learning, where the objective is not to explicitly classify the input patterns but rather to discover internal representations by fitting a hierarchical generative model to the sensory data (Hinton, 2007, 2013; Zorzi et al., 2013). Compared to its supervised counterpart, this modeling approach emphasizes the role of feedback, recurrent connections (Sillito et al., 2006), which carry top-down expectations that are gradually adjusted to better reflect the observed data (Hinton and Ghahramani, 1997; Friston, 2010) and which can be used to implement concurrent probabilistic inference along the whole cortical hierarchy (Lee and Mumford, 2003; Gilbert and Sigman, 2007). Notably, topdown processing is also relevant for understanding attentional mechanisms in terms of modulation of neural information processing (Kastner and Ungerleider, 2000).

A powerful class of stochastic neural networks that learn a generative model of the data is that of Restricted Boltzmann Machines (RBMs), which can efficiently discover internal representations (i.e., latent features) using Hebbian-like learning mechanisms (Hinton, 2002). RBMs constitute the building block of hierarchical generative models such as Deep Belief Networks (Hinton and Salakhutdinov, 2006) and Deep Boltzmann Machines (Salakhutdinov, 2015). These unsupervised deep learning models have been successfully used to simulate a variety of cognitive functions, such as numerosity perception (Stoianov and Zorzi, 2012), letter perception (Testolin et al., under review), location-invariant visual word recognition (Di Bono and Zorzi, 2013), and visual hallucinations in psychiatric syndromes (Reichert et al., 2013). A similar approach has been used to simulate how early visual cortical representations are adapted to statistical regularities in natural images, in order to predict single voxel responses to natural images and identify images from stimulus-evoked multiple voxel responses (Güçlü and van Gerven, 2014). A temporal extension of RBMs has also been recently used to model sequential orthographic processing and spontaneous pseudoword generation (Testolin et al., 2016).

Unsupervised deep learning can be implemented using an alternative architecture based on autoencoders (Bengio et al., 2007), which are deterministic, feed-forward networks whose learning goal is to accurately reconstruct the input data into a separate layer of output units. Single-layer autoencoders are trained using error backpropagation, and can be stacked in order to build more complex, multi-layer architectures. However, despite the common view that RBMs and autoencoders could be considered equivalent (Ranzato et al., 2007), we note that their underlying architectural and learning assumptions are significantly different. In this study we empirically compare RBMs and autoencoders in terms of the type of internal encoding emerging in the hidden neurons. Moreover, we investigate how additional learning constraints, such as sparsity and limitation of computational resources (i.e., hidden layer size), could influence the representations developed by the networks. As a case study, we focus on the problem of learning visuospatial coding for sensorimotor transformations, which is a prominent example of how the emergentist approach based on learning in artificial neural networks has offered important insights into the computations performed by biological neurons (Zipser and Andersen, 1988).

Sensorimotor transformations refer to the process by which sensory stimuli are converted into motor commands. For example, reaching requires to map visual information, represented in retinal coordinates, into a system of coordinates that is centered on the effector. Coordinate transformations can be accomplished by combining sensory information with extra-retinal information, such as postural signals representing the position of eyes, head, or hand, thereby obtaining abstract representations of the space interposed between the sensory input and the motor output (Pouget and Snyder, 2000). Single-neuron recordings from monkey posterior parietal cortex have shown that the response amplitude of many neurons indeed depends on the position of the eyes, thereby unveiling a fundamental coding principle used to perform this type of signal integration (Andersen et al., 1985). The term gain field was coined to describe this gaze-dependent response of parietal neurons, and since then the notion of gain modulation has been generalized to indicate the multiplicative control of one neuron's responses by the responses of another set of neurons (Salinas and Thier, 2000). Another fundamental property unveiled by neuronal recordings is that the encoding of space used for coordinate transformations involves a variety of different, complementary frames of reference. For example, although many parietal neurons are centered on retinal coordinates (Andersen et al., 1985; Duhamel et al., 1992), others represent space using body-centered (Snyder et al., 1998) or effector-centered (Sakata et al., 1995) coordinate systems. Moreover, some neurons exhibit multiple gain modulation (Chang et al., 2009), suggesting more complex forms of spatial coding. For example, postural information related to both eye and head positions can be combined in order to encode "gaze" direction (Brotchie et al., 1995; Stricanne et al., 1996; Duhamel et al., 1997).

From a computational perspective, the seminal work of Zipser and Andersen (1988) showed that gain modulation could spontaneously emerge in supervised, feed-forward neural networks trained to explicitly map visual targets into headcentered coordinates, giving as input any arbitrary pair of eye and retinal positions. Similar results have been observed using more biologically-plausible learning settings, such as reinforcement learning (Mazzoni et al., 1991) and predictive coding (De Meyer and Spratling, 2011). Note that these learning settings assume that gain modulation emerges because the task implies to establish a mapping between different reference frames. However, it is unclear whether the form of modulation and the distribution of neuronal tuning functions is influenced by the type of learning algorithm and/or by the nature of the learning task (i.e., learning input-output mappings vs. unsupervised learning of internal representations). We also note that a popular alternative framework for modeling sensorimotor transformations is not based on learning, but rather stipulates that parietal neurons represent a set of basis functions that combine visual and postural information (for review, see Pouget and Snyder, 2000).

In summary, space coding represents an interesting case study for testing the adequacy of different neural network architectures and learning algorithms, because it provides a wealth of neurophysiological data (both at the population and single-neuron levels), and it departs from the classic problem of visual object recognition investigated in the large majority of deep learning research.

## MATERIALS AND METHODS

In this section we describe the space coding tasks used in our simulations, including training and test stimuli, the different learning architectures, and the procedures for analyzing the emergent neural representations.

## Space Coding Tasks

In this study we consider a visual signal in retinotopic coordinates and two different postural signals, one for eye position and another for a generic "effector," which might represent, for example, the position of the hand. We do not consider the integration between different modalities (see Xing and Andersen, 2000, for a computational investigation of multimodal integration in several coordinate frames). We implemented three types of space coding tasks to test the different learning architectures.

#### Unsupervised Learning with No Coordinate Transformation

The first learning architecture is depicted in **Figure 1A**. Unsupervised learning is represented by undirected arrows, which connect the sensory input to a separate layer of hidden neurons. The input signal to the network consists of a visual map, which represents target location in retinotopic coordinates, and two postural maps, which represent eye and effector positions. The learning goal is only to build a compact representation of these input signals in the hidden layer, which is later read-out by a simple linear associator in order to establish a mapping with the corresponding motor program. Details of input and output representations are provided in Section Dataset and Stimuli. The unsupervised learning phase does not involve any coordinate transformation because information about the motor program is not available.

#### Unsupervised Learning with Coordinate Transformation

The second learning architecture is depicted in **Figure 1B**. The input signal to the network still consists of a visual map and two postural maps, but in this case we also provide as input the corresponding motor program. In this setting the unsupervised learning phase implicitly involves coordinate transformation (i.e., different coordinate systems become associated). In order to compare the mapping accuracy of different learning architectures using the same method, the motor program is still read-out from hidden neurons via a simple linear associator.

#### Supervised Learning with Coordinate Transformation

The third learning architecture is depicted in **Figure 1C**, and it corresponds to the model used by Zipser and Andersen (1988). The input is the same of the unsupervised architecture shown in **Figure 1A**, but in this case supervised learning (directed arrows) is used to establish an explicit mapping between input signals

and motor programs. As for the previous architectures, accuracy of the motor program is also tested by read-out from hidden neurons via linear association.

#### Dataset and Stimuli

The representation format adopted for the sensory stimuli was the same used in previous computational investigations (Zipser and Andersen, 1988; Pouget and Snyder, 2000; De Filippo De Grazia et al., 2012), which is broadly consistent with neurophysiological data recorded in animals performing tasks involving coordinate transformations (e.g., Andersen et al., 1985).

The visual input to the models consisted in a real-valued vector representing the position of the stimulus as a Gaussian peak of activity in a specific location. These visible neurons simulate the activity of the cortical areas supplying retinotopic sensory information to the posterior parietal cortex. The retinotopic map consisted in a square matrix of 17 × 17 neurons, which employed a population code with Gaussian tuning functions (standard deviation = 4 ◦ ). Visual receptive fields were uniformly spread between −9 ◦ and +9 ◦ with increments of 3◦ , both in the horizontal and vertical dimensions.

Four postural maps, each one consisting of 17 neurons, were used to represent the horizontal and vertical positions of the eye and the effector. These visible neurons used a sigmoid activation function (steepness parameter = 0.125) to represent postural information between −18 and +18◦ , with steps of 3◦ .

The motor program consisted in a real-valued vector representing the target position of the stimulus. Similarly to the retinotopic map, it was coded as a square matrix of 25 × 25 neurons, which employed a population code with Gaussian tuning functions to represent target position in coordinates centered on the effector (standard deviation = 6 ◦ ). Motor programs were uniformly spread between −9 ◦ and +9 ◦ with increments of 3◦ , both in the horizontal and vertical dimensions.

In order to create the stimuli dataset, all possible combinations of visual input and postural signals were first generated, and the corresponding motor program (target location) was computed. We then balanced the patterns to ensure that target locations were equally distributed across the motor map to avoid position biases when decoding the motor program. This resulted in a total of 28,880 patterns, which were randomly split into a training set (20,000 patterns) and an independent test set (8,880 patterns). The latter was used to assess the generalization performance of the models.

## Learning Architectures

Despite they differ in several aspects, Boltzmann machines and autoencoders can both be defined within the mathematical framework of energy-based models (Ranzato et al., 2007), where the learning objective is to carve the surface of an energy function so as to minimize the energies of training points and maximize the energies of unobserved points. A set of latent variables is used to learn an internal code that can efficiently represent the observed data points, and since the number of latent variables is usually smaller than that of the observed variables the encoding process can be interpreted as a form of dimensionality reduction (Hinton and Salakhutdinov, 2006). In this unsupervised setting, the model learns the statistical structure of the data without the need for any explicit, external label.

#### Restricted Boltzmann Machines (RBMs)

Boltzmann machines are stochastic neural networks that use a set of hidden neurons to model the latent causes of the observed data vectors, which are presented to the network through a set of visible neurons (Ackley et al., 1985). In the "restricted" case, the network connectivity is constrained in order to obtain a bipartite graph (i.e., there are no connections within the same layer; see **Figure 2A** for a graphical representation). The behavior of the network is driven by an energy function E, which defines the joint distribution of the hidden and visible neurons by assigning a probability value to each of their possible configurations:

$$p(\nu, h) = \frac{e^{-E(\nu, h)}}{Z}$$

where v and h are the column vectors containing the values of visible and hidden neurons, respectively, and Z is the partition function. The energy function is defined as a linear combination of visible and hidden neurons' activation:

$$E(\nu, h) = -b^T \nu - c^T h - h^T W \nu$$

where W is the matrix of connections weights, b and c are two additional parameters known as unit biases and T denotes the transpose operator. Since there are no connections within the same layer, hidden neurons are conditionally independent given the state of visible neurons (and vice versa). In particular, the activation probability of the neurons in each layer conditioned on the activation of the neurons in the opposite layer can be efficiently computed in one parallel step:

$$\begin{aligned} P(h\_j = 1 | \nu) &= \sigma(c\_j + \sum\_i \boldsymbol{w}\_{ij} \boldsymbol{\nu}\_i) \\ P(\boldsymbol{\nu}\_i = 1 | h) &= \sigma(b\_i + \sum\_j \boldsymbol{w}\_{ij} h\_j) \end{aligned}$$

where σ is the sigmoid function, cj and bi are the biases of hidden and visible neurons (hj and vi respectively), and wij is the connection weight between hj and vi . Learning in RBMs can be performed through maximum-likelihood, where each weight should be changed at each step according to a Hebbian-like learning rule:

$$
\Delta W = \eta (\nu^+ h^+ - \nu^- h^-)
$$

where η represents the learning rate, v +h + are the visible-hidden correlations computed on the training data (positive phase), and v −h − are the visible-hidden correlations computed according to the model's expectations (negative phase). Model's expectations have been traditionally computed by running Gibbs sampling algorithms until the network reached equilibrium (Ackley et al., 1985). However, more efficient algorithms such as contrastive divergence (Hinton, 2002) speed-up learning by approximating the log-probability gradient. The reader is referred to Hinton (2010) and Zorzi et al. (2013) for more details about RBMs and for the discussion of hyper-parameters of the learning algorithm.

In our simulations, RBMs were trained using 1-step contrastive divergence with a learning rate of 0.03, a weight decay of 0.0002 and a momentum coefficient of 0.9, which was initialized to 0.5 for the first few epochs. Learning was performed using a mini-batch scheme, with a mini-batch size of 4 patterns, for a total of 100 learning epochs (reconstruction error always converged). Sparse representations were encouraged by forcing the network's internal representations to rely on a limited number of active hidden units, that is, by driving the probability q of a unit to be active to a certain desired (low) probability p (Lee et al., 2008). For logistic units, this can be practically implemented by first calculating the quantity q-p, which is then multiplied by a scaling factor and added to the biases of each hidden units at every weight update. When the sparsity constraint was applied, we always verified that the average activation of hidden units was indeed maintained below the desired level. All the simulations were performed using an efficient implementation of RBMs on graphic processors (Testolin et al., 2013). The complete source code is available for download1 .

#### Autoencoders

Similarly to RBMs, autoencoders rely on a single layer of nonlinear hidden units to compactly represent the statistical regularities of the training data. However, autoencoders are feed-forward, deterministic networks trained with error backpropagation (Bengio et al., 2007). The training data is presented to a layer of input units, and the learning goal is to accurately reconstruct such input vector into a separate, output layer. An autoencoder is therefore composed of a set of encoding weights W1 that are used to compute the activation of hidden h units given the activation of input units v, and a set of decoding weights W2 that are used to compute the network reconstructions v\_rec from the activations of hidden units:

$$\begin{aligned} h &= \sigma(W^1 \nu + c) \\ \nu\_- \text{rec} &= \sigma(W^2 h + b) \end{aligned}$$

where b and c are the vectors of output and hidden unit biases, and σ is the sigmoid function (see **Figure 2B** for a graphical representation). The error function E to be minimized corresponds to the average reconstruction error, which is quantified by the sum across all output units of the squared difference between the original and the reconstructed values:

$$E = \frac{1}{N} \sum\_{n=1}^{N} \sum\_{k=1}^{K} \left(\nu\_k - \nu\_{-} rec\_k\right)^2 + \beta^\* \Omega\_{sparsity}$$

where K is the number of output units and N is the number of training patterns. Similarly to RBMs, sparse representations can be induced by adding to the cost function a regularization term sparsity that takes a large value when the average activation value q of each hidden neuron diverges from a certain desired (low) value p. In particular, the sparsity constraint was implemented as the Kullback-Leibler divergence from q to p:

$$\Omega\_{\text{sparsity}} = \sum\_{i=1}^{H} \text{KL}(p \parallel q\_i),$$

where H is the number of hidden units. As for RBMs, when sparsity was applied we always verified that the average activation of hidden units was indeed maintained below the desired level.

1http://ccnl.psy.unipd.it/research/deeplearning

In our simulations, we used an efficient implementation of autoencoders provided by the MATLAB Neural Network toolbox (Demuth and Beale, 1993). Learning was performed using standard scaled conjugate gradient descent (Møller, 1993) with adaptive learning rate, using a weight decay factor of 0.0002 and a batch processing scheme, for a total of 150 learning epochs (reconstruction error always converged).

#### Feed-Forward, Supervised Networks

In order to better assess the impact of the learning regimen, we compared the unsupervised learning architectures described above with a standard, supervised architecture implemented as a feed-forward network with one hidden layer (Zipser and Andersen, 1988). Similarly to autoencoders, learning can be performed using error backpropagation (see **Figure 2C** for a graphical representation). We used an efficient implementation of feed-forward networks provided by the MATLAB Neural Network toolbox2 . Learning rate was set to 0.05 and training was performed for a total of 2500 learning epochs (output error always converged).

### Testing Procedure

For each experimental setting, we run 10 different networks in order to collect simulation statistics. In the results, we therefore always report mean values along with standard deviations.

#### Decoding Internal Representations by Linear Read-Out

Following unsupervised learning, a linear read-out was performed from the internal (hidden layer) distributed representations of the networks in order to assess how well they could support a supervised mapping to the target motor program through a simple linear projection (Pouget and Snyder, 2000). The read-out was implemented using a linear neural network trained with the delta rule (Widrow and Hoff, 1960). Learning was performed for 250 epochs using mini-batches of 20 patterns. Learning rate was set to 0.07, and weight decay of 0.000001 was used as a regularizer. Classifier performance was always measured on the separate test set. Test errors always matched those obtained on the training set, indicating that the read-out was robust to overfitting.

The output of the classifier was first compared with the target motor program by computing the Root Mean Squared Error (RMSE) between the two matrices. However, a more useful performance measure was obtained by first decoding the Center Of Mass (COM) of the output distribution, which was then compared with the actual coordinates of the motor program. This measure allows to quantify the read-out error in degrees: following Zipser and Andersen (1988), the mapping was considered to be successful if the error was below the distance between the centers of the Gaussian tuning functions in the retinotopic map (i.e., 3◦ ). If the latter mapping accuracy was not achieved, we did not consider the network for subsequent analyses. We found the RMSE and COM measures to be always consistent with each other, so we only report COM results.

Measuring Single-Neuron and Population Sparseness An index of single-neuron sparseness was computed using a well-established procedure employed in neurophysiological investigations (Rolls and Tovee, 1995; Vinje and Gallant, 2000), which describes the activity fraction a of each neuron across stimuli as:

$$a = \frac{\left(\sum r\_i/n\right)^2}{\sum \left(r\_i^2/n\right)}$$

where ri is the firing rate of the neuron to the i-th stimulus in the set of n stimuli. This is a useful measure of the extent of the tail of the distribution, in this case of the firing rates of the neuron to each stimulus. Mean single-neuron sparseness