Designing spiking neural networks for robust and reconfigurable computation

How can spiking neural networks compute in ways that are robust against external disruptions such as the loss of individual components and how can they be reconfigured to compensate for disruptions? We here study these questions for basic computational tasks known as $k$-winner-take-all computations that are generalizations of rank ordering and that may be combined to yield universal forms of computations.

Rank ordering of signals by their (average) strength is a fundamental computational operation, particularly useful in attention-related tasks in both natural and artificial systems.^1–3 It provides an effective way for extracting the most important information from high-dimensional input spaces by simply ordering in terms of relevance (or strength). However, such operations are computationally costly for high-dimensional inputs, and may retain, at times, a large amount of irrelevant information (depending on the application). A less costly, but closely related, operation is called partial rank ordering, in which a subset of the $k$ strongest out of $N>k$ signals is identified. Often such operations are referred to as $k$-winner-takes-all (⁠ $k$-WTA) computations. Still, already for fixed $k$⁠, combinatorially many computational outputs (results of the task) need to be accounted for by any system that performs $k$-WTA computations and the number $k$ of “winners” might be variable in addition.

Besides the mathematical analysis of $k$-WTA computations, a variety of systems have been proposed for their implementation. Particularly, partial rank ordering can be performed by neuron networks.^4,5 For example, to date, there exist multiple bio-inspired implementations of one-winner-takes-all (⁠ $1$-WTA) functionality.^6–8 More general $k$-WTA operations can be performed by exploiting complex periodic orbits in symmetrical oscillator networks, similar to heteroclinic computing.^5,9–15 Furthermore, simple hardware implementation of a neural-circuit performing WTA-calculations has been recently suggested,¹⁶ exhibiting a mixture of excitatory and inhibitory couplings. However, with all these approaches, either the number of winners $k$ is fixed or not easily reconfigurable, or computations typically take a long time.^{5,10–12,14,15,17,18}

Recently, a fast and re-configurable $k$-WTA implementation via a symmetric neural network with inhibitory pulse coupling was proposed.¹⁹ Via adjusting a single parameter, the global coupling strength, the number of winners $k$ can be chosen freely and hence be adapted to different scenarios. In contrast to most existing $k$-WTA implementations, the system is made up from very simple, identical parts and usually converges only within a few spikes. The network presented in Ref. 19 is a computational application of a specific multiplicative coupling scheme in which the effect of an incoming pulse depends linearly on the state of the receiving oscillator. Although closely connected to more common coupling types where the strength of the inhibition does not depend on the oscillator voltage,²⁰ this multiplicative coupling is the most natural way to achieve the specific type of phase compression which is necessary for the $k$-WTA functionality.

In this article, we analytically characterize the mechanisms of $k$-WTA computations based on proportional inhibition. We adapt the original Mirollo–Strogatz formalism²¹ to multiplicative pulse-coupling to analytically describe the dynamics of the $k$-WTA network proposed in Ref. 19 and explore under which conditions it exhibits the desired computational features. In particular, we identify parameter regions in which the system may perform specific computations and provide analytic expressions for the appropriate coupling strength as a function of external variables. In contrast to the seminal work,¹⁹ which uses “brute-force” parameter scans and numerical integration, the current work addresses the question of how to set the global coupling strength in order to obtain collective dynamics with a specified number of $k$ firing neurons with a more analytical approach, hence potentially allowing for a more direct (re-)calibration and thus reconfiguration of the computational network.

It has been shown that symmetrical networks of such neuronal units may be used to perform $k$-WTA calculations.¹⁹ An idealized diagram of a complete computing unit is shown in Fig. 1(a). In this context, the variables $I i$ [see Eq. (2)] are taken as input signals, which determine the free dynamics and period lengths $T i$ of the corresponding neurons, see Fig. 1(b). The reset pulses of the system not only mediate the global coupling, but also provide the output of the computational system: After setting the input signals $I i$⁠, the network converges to a collective dynamic in which only the $k$ fastest neurons reset, leaving the other outputs silent. The underlying mechanism can be summarized as follows. Whenever some neuron $i$ resets, the down-concavity of $f i IF( x i)$⁠, together with the proportional coupling as described by Eq. (5) leads to a relative compression of the voltages $x i$ of all neurons $i≠j$⁠, see Fig. 1(c) (for details, see Ref. 19). If the differences in input signals $I i$ and the global coupling strength $ϵ$ are both sufficiently large, neurons with shorter free period $T i$ may overtake other neurons repeatedly, keeping them from reaching the threshold value altogether by recurrently inhibiting them, see Fig. 1(d). After a (typically short) transient, the system enters a periodic orbit, in which only the $k$ fastest out of all $N$ neurons reset and send pulses, while the neurons which correspond to the $N$– $k$ smallest input signals stay silent. Within any periodic orbit, the spike patterns of the individual neurons are interpreted as output vector by discriminating only between spiking and non-spiking neurons, so that the $N$-dimensional, real-valued input vector $( I 1 , … , I N ) T$ is mapped to a binary $N$-dimensional vector consisting of $k$ ones and $N$– $k$ zeros, identifying the $k$ highest and $N$– $k$ lowest input signals. The exact number of spiking neurons $k$ depends on the particular choice of input signals $I i$⁠, as well as on the global coupling strength $ϵ$⁠.

As outlined in Ref. 19, the considered spiking neural networks can be reconfigured to different required numbers of winners $k$ and input signal differences by adapting one global parameter (the coupling strength) after, for example, the loss of a neuron (and its connections). Under some conditions on the input signal vector and the system parameters, the disrupted and reconfigured network is then capable of performing the same (⁠ $k$-WTA) computation as the original non-disrupted network.

Figure 2(a) shows the number of winners $k$ as a function of the input spacing parameter $ΔI$ and the coupling strength $ϵ$⁠. For each choice of $ϵ$ and $ΔI$⁠, the collective dynamics is evaluated, starting from 100 different random initial conditions. Throughout, we find the number of winners $k$ to be independent from the specific initial conditions, so that the parameter space consists of cohesive regions with unique $k$⁠. In particular, we find that a wide range of parameter configurations leads to period-one dynamics. In these, every neuron spikes at most once, so that the sequence length $p$ of the resulting periodic orbit (i.e., the number of reset events per orbit) is equal to the number of spiking neurons, $p=k$⁠. However, period-one dynamics with different $k$ are usually divided by transitional regions in which neurons spike more than once per period, $p>k$⁠, see the close-up in Fig. 2(b).

In general, the exact structure of the dynamical space might depend strongly on system parameters and on the range of input signals $I i$⁠. Meanwhile, “brute-force” parameter scans and direct numerical integration as done for Fig. 2 and in Ref. 19 can be computationally expensive.

Hence, in the following, we develop an analytical framework which allows to find the system parameters which exhibit specific types of orbits in a more direct manner, without simulating the underlying system. For instance, given fixed input signals, as, e.g., given by Eq. (6), which coupling strengths enable computations for which $k$⁠.

As the transitions between different $k$ are closely connected to period-one orbits (or, more precisely, their absence), we put a particular emphasis on the description of period-one orbits, before generalizing our description to arbitrary dynamics.

In the following sections, we derive sets of analytical conditions to determine whether a specific periodic orbit is consistent, given specific neuron models, neuron frequencies ${ ω i}$⁠, and coupling strength $ϵ$⁠. In doing so, we follow the general approach proposed in Ref. 24. Solving these conditions for the global coupling strength $ϵ$ yields the adequate choices for the desired dynamics. In particular, we identify coupling strengths which lead to collective dynamics with a specific number of winners $k$⁠, and hence calibrate the $k$-WTA system more directly, with no need for numerical simulations of the collective dynamics.

Every periodic orbit can be characterized by its sequence of neuron resets $s i$⁠, $i∈{1,…p}$⁠, where the sequence length $p$ is the total number of reset events in the given orbit. For mathematical analysis, it is convenient to encode the reset order of neurons in terms of a mapping $σ:(1,…,p)↦σ(i)≡ σ i∈(1,…,N)$⁠, such that $( s 1,…, s p)=( σ 1,…, σ p)$⁠. We also use $σ$ to refer to the orbit itself. The subset of spiking neurons within a specific orbit $σ$ is $S≡ S σ={ σ i |i∈{1,…,p}}$⁠, which contains $k$ neurons. The complement of $S σ$⁠, $S ¯≡ S ¯ σ={1,…,N}∖ S σ$ contains the $N$– $k$ neurons which do not reset in the considered orbit and conclusively remain silent.

Given a full description of the system parameters, the event sequence $s 1,…, s p$⁠, or equivalently, the mapping $σ$⁠, uniquely defines a physical orbit, that is, determines phases $ϕ m∈{1,…,N}$ of all neurons at any times, as well as the event times $t i$⁠, $i∈{1,…,p}$ at which the $i$th neuron in the considered orbit, neuron $σ i$⁠, spikes, relative to the time $t 1=0$ of the first reset event. We denote the time step between two successive reset events at times $t i$ and $t i + 1$ of neurons $σ i$ and $σ i + 1$⁠, respectively, as $Δ t i= t i + 1− t i$⁠. As the orbits are periodic by definition, we treat $σ$ as cyclical, so that $σ i + p≡ σ i$ and $Δ t i + p=Δ t i$⁠. Also note that the description of a physical periodic orbit in terms of a mapping $σ$ are unique only up to cyclical permutations of $σ$⁠.

The approach for finding system parameters which allow for a specific orbit as given by $σ$ is as follows. First, we derive a set of $p$ periodicity conditions, each one guaranteeing the periodicity of a single reset event $s i$⁠, $i∈{1,…,p}$ at time $t i$⁠. Second, we have to ensure that the assumed spike ordering is consistent in so far that all time steps $Δ t i$ are positive. Third, we have to guarantee that no neuron $l∈ S ¯$ reaches threshold when the system is in the orbit $σ$⁠, leading to another set of inequalities. Finally, we have to ensure that neurons $j∈ S$ reach the threshold only at reset times $t i$ with $σ i=j$⁠, leading to additional inequalities.

In Secs. V A–V C, we exemplarily demonstrate the approach for period-one dynamics, that is, orbits where each of the $k$ spiking neurons resets exactly once, so that the number of reset events $p$ is equal to the number of winners $k$⁠. In doing so, we illustrate that typically violations of the inequalities correspond to transitions between different dynamical regions or numbers of winners $k$⁠. For a discussion of how to treat arbitrary orbits, we refer to  Appendix A.

In the following, we provide an analytical description of period-one orbits with fixed $k$ and derive expressions for the transitions to other dynamical regimes. In doing so, we explore the microscopic mechanisms mitigating these transitions and demonstrate how to characterize parameter choices that lead to a desired computational functionality.

Note that Eq. (22) closely resembles the periodicity Eq. (18), but without a reset before the first (or after the last) line. Also note that because of $ϕ ℓ ( t i − )>0$⁠, the phases at later times are always larger than they would be with a reset at time $t i$⁠, given fixed ${Δ t j}$ and free frequency $ω ℓ$⁠. This implies that if a neuron $m$ with free frequency $ω m$ is reset in a given orbit, $m∈ S$⁠, also all other neurons $n$ with $ω n≥ ω m$ must reset in the given orbit, as otherwise condition (20) is violated for $l=n$ at the time where neuron $m$ spikes. Indeed this is a central property of all self-consistent periodic orbits for the considered system type: The $k$ spiking neurons in a given orbit $k$ are always the $k$ neurons with the largest input signal and therefore the largest free frequency. This also implies that if two neurons $m$ and $n$ have the same frequency $ω m= ω n$⁠, it is not possible to have one neuron spike and the other one be silent, regardless of the coupling strength, restricting the possible options for the number of winners $k$⁠.

Note that the critical value $ϵ min ( 2 )$ does not depend on the order in which the neurons spike in the considered orbit so that for $ϵ≤ ϵ min ( 2 )$⁠, no period-one dynamics with given $k$ are consistent, regardless of the order of spiking neurons. In fact, also non-period-one orbits with sequence length $p>k$ with given $k$ are not possible below the critical value $ϵ min ( 2 )$⁠, as these require even higher coupling strengths than the period-one dynamics with the same $k$⁠, compare Fig. 2. As typically $ϵ min ( 2 )≫ ϵ min ( 1 )$ (for $k<N$⁠) and $ϵ min ( 2 ) ≫ ϵ min ( 3 )$ (see the next section), the critical value $ϵ min ( 2 )$ usually is the lower boundary of the coupling strength for a fixed number of winners $k$⁠, regardless of the orbit type, see Fig. 3.

What happens on a mechanistic level, if one starts with a consistent orbit satisfying condition (37) and continuously decreases the coupling strength, $ϵ→ ϵ min ( 2 ) +$? According to Eq. (39), the maximum value of the phase $ϕ ℓ(t)$ within the orbit $σ$ continuously increases, and formally, at coupling strength $ϵ= ϵ min ( 2 )$⁠, takes on $1$ just at as the $k$-fastest neuron (neuron $s$⁠) spikes. However, as the reset pattern $σ$ does not account for a reset of neuron $ℓ$⁠, at $ϵ= ϵ min ( 2 )$ the periodic orbit $σ$ becomes inconsistent and a discontinuous change to an orbit with $k+1$ spiking neurons, in which also neuron $ℓ$ spikes, occurs. Typically, the coupling strength $ϵ= ϵ min ( 2 )$ is still too high for period-one dynamics with the new number of winners $k+1$⁠, so that we observe a more complex orbit with $p>k$⁠, in which faster neurons spike multiple times, while two or more slower neurons take turns. For a typical example, see Fig. 3.

As we are interested in values of the coupling strength $ϵ$ for which condition (45) is satisfied for all tuples $i,j∈{1,…,k},j≠i$⁠, for given $k$⁠, ${ ω i}$⁠, and spiking order $σ$ we have to check condition (45) only for the tuple $(i,j)$ that maximizes $ϕ j ( t i − )$⁠.

Which ordering of neurons represents the most stable period-one orbit fulfilling condition (45), that is, gives the smallest maximum phase $ϕ j ( t i − )$? Considering Eq. (44) again, we note that $ϕ j ( t i − )$ is larger for larger $μ$ as well as for larger ratios $ω σ i/ ω σ j$⁠. We find the most stable reset pattern permutation to be the one in which all neurons reset in ascending order of their free frequencies and then start with the slowest neuron again (standard ordering), see Fig. 3. This is because for the ascending patterns, the smaller the number of reset events $μ$ from $i$ to $j$⁠, the larger the ratio $ω σ i/ ω σ j$ becomes such that both contributions compete. Hence, in order to identify values of $ϵ$ where any type of period-one dynamics with given $k$ is possible, we consider only orbits with ascending free frequencies.

In  Appendix C, we show that, for any period-one orbits with standard ascending ordering, if Eq. (45) is satisfied for $μ=k−1$⁠, it also holds for any other $μ∈{1,…,k−2}$⁠, so that we have to consider only $μ=k−1$ for evaluating condition (45). The interpretation is that the first violation of the condition (45) happens in terms of some neuron $i$ reaching threshold exactly one event “too early.”

What happens at the point where condition (46) is violated? For $ϵ→ ϵ max ( 3 ) −$⁠, the phase of at least one neuron $σ i∈ S$ starts to reach the threshold $1$ at some time $t j$ which does not correspond to the proper reset event $s i$ as defined by the mapping $σ$⁠. Depending on the input configuration and system details, at $ϵ= ϵ max ( 3 )$ either a more complicated orbit with $p>k$⁠, but the same number of winners $k$⁠, occurs, or the systems dynamics changes directly to period-one dynamics with a new number of winners $k+1$⁠. However, for $k>2$⁠, we usually observe the first case. For a typical example, consider Fig. 3.

Figure 4 illustrates the identification of parameter regions with period-one dynamics for linear neurons with equally spaced input signals $I i=(i−1)ΔI$⁠, for different choices of the input spacing parameter $ΔI$⁠. Figure 4(a) shows both simulation results as well as the derived boundaries for period-one dynamics for different numbers of winners $k$⁠. Indeed, the analytical boundaries are consistent with the simulation results. Figure 4(b) focuses on a smaller parameter region with $k∈{5,6,7}$⁠, illustrating the transition between different $k$⁠. While the lower boundary of the coupling strength $ϵ$ for period-one dynamics with given $k$ (solid black line) corresponds also to a transition to a larger value of $k$⁠, the upper boundary of the coupling strength (dashed black line) leads to more complicated orbits with $p>k$⁠, but with the same number of winners. For $k<N$⁠, we find $ϵ min= ϵ min ( 2 )$ throughout.

We point out that for linear neuron dynamics a closed-form solution for the time steps $Δ t i$⁠, $i∈{1,…,p}$ is not only possible for period-one dynamics, but also for arbitrary complex dynamics with $p>k$⁠, see  Appendix B.

For nonlinear free (uncoupled) neuron dynamics, or more specifically, free neuron models which lead to nonlinear transfer functions, the constraints can generally not be solved in closed form. However, the corresponding equations can be solved for the critical coupling strengths $ϵ$ easily by using standard methods for root finding. As such nonlinear neuron dynamics are commonly used for modeling biological systems, as well as in engineering application, we briefly describe the procedure, which treats each set of constraints (1), (2), and (3) [Eqs. (19), (20), and (23)] separately.