Approximating Nonlinear Functions With Latent Boundaries in Low-Rank Excitatory-Inhibitory Spiking Networks

The neural circuits of the brain are unbelievably complex. Yet when it comes to studying how they compute, we often resort to highly simplified network models composed of neurons with graded activation functions, that is, rate neurons. The resulting rate networks have become the standard models of feedforward sensory processing (Yamins & DiCarlo, 2016; Lindsay, 2021) and recurrent task dynamics (Sussillo, 2014; Barak, 2017) and capture many aspects of neural circuits surprisingly well.

Of course, a mechanistic understanding of biological computation must eventually bridge back to the details of real circuits (Bernáez Timón et al., 2022). However, this task has proved surprisingly hard. Indeed, the more biologically detailed a network model is, the more difficult it tends to be to constrain and interpret (Eliasmith & Trujillo, 2014). We argue here that a fundamental part of this problem lies in two computational concepts of rate networks that are mismatched with biology.

The first concept is that of the feature detector (Martin, 1994) and pertains to the difference between rate-based and spike-based coding. In a nutshell, rate neurons operate in a regime of depolarized inputs; they are activated far beyond threshold whenever the neuron’s input pattern matches its pattern of synaptic weights. In the spiking domain, unless a large amount of external noise is added, a direct translation of this idea leads to regular spike trains (Softky & Koch, 1993; see Figures 1a–1c). While this may be an accurate description of some biological neurons (e.g., those at the sensory periphery), cortical neurons often fire spikes irregularly (Softky & Koch, 1993). The cause of this irregularity is that such neurons operate in a fluctuation-driven regime, in which excitatory and inhibitory input currents balance each other on average (Van Vreeswijk & Sompolinsky, 1996; Shadlen & Newsome, 1998; Haider et al., 2006), and not in a strongly depolarized regime (see Figures 1d–1f). However, despite much progress (see section 5 for more details), a general theory of computation in such a regime has remained elusive, especially if computations are limited to smaller networks with only tens of neurons.

The second mismatched concept of rate networks concerns function approximation and its relationship to Dale’s law, the common biologically observed distinction between excitatory and inhibitory neurons (Eccles, 1976). The flexibility of rate networks relies on the ability of each unit to linearly combine inputs in order to represent arbitrary input-output transformations on the network level (see Figure 1c). To be most effective, the neurons’ output rates must be combined with both positive and negative weights, resulting in mixed-sign connectivity on the level of individual neurons, and thus violating Dale’s law. Several studies have successfully incorporated Dale’s law (e.g., Parisien et al., 2008; Song et al., 2016; Miconi, 2017; Ingrosso & Abbott, 2019; Shao & Ostojic, 2023), and some have even suggested potential benefits to learning and robustness (Haber & Schneidman, 2022). However, these studies have primarily taken a bottom-up approach, incorporating Dale’s law for biological plausibility rather than for any computational necessity. It also remains to be seen how such sign constraints scale, as adding them to larger-scale machine learning benchmarks typically hinders (or at best matches) performance (Cornford et al., 2021; Li et al., 2023).

Neither of these mismatches is new, and both have puzzled the field for a while. One common line of reasoning suggests that irregular firing and Dale’s law are two among many manifestations of biological constraints on computation. From this view, one could justify the use of abstract rate networks as idealized versions of a constrained, noisy biological implementation. However, it is also possible that these mismatches reflect something fundamentally different about biological computation. Taking the latter perspective, here we offer a fresh approach to these puzzles.

Our study is based on a combination of two recent developments: low-rank connectivity and spike-threshold geometry. Networks with low-rank connectivities generate low-dimensional dynamics in a latent activity space (Seung, 1996; Eliasmith & Anderson, 2003; Landau & Sompolinsky, 2018; Mastrogiuseppe & Ostojic, 2018). They thereby allow the activity of individual neurons to be linked to the modes or patterns of population activities often observed in real neural circuits (Gallego et al., 2017; Saxena & Cunningham, 2019; Keemink & Machens, 2019; Jazayeri & Ostojic, 2021; Chung & Abbott, 2021; Langdon et al., 2023). Such insights can also be translated into spiking network models (Eliasmith & Anderson, 2003; Cimeša et al., 2023; Koren & Panzeri, 2022; DePasquale et al., 2023).

However, low-rank connectivity alone does not guarantee spike-based computation in the balanced operating regime (Eliasmith & Anderson, 2003). To achieve this, we take inspiration from spike-coding networks (SCNs; Boerlin et al., 2013), whose function can be geometrically understood by visualizing spike-threshold boundaries in a space of latent population modes (Calaim et al., 2022). Importantly, this perspective places the boundary between sub- and suprathreshold voltages at the center of computation, leading to a fundamentally different operating regime from rate networks (see Figures 1a and 1d). Such networks have been shown to exhibit irregular activity and have several other desirable biological properties like robustness and energy efficiency (Denève & Machens, 2016; Barrett et al., 2016). Based on a geometric reframing (Calaim et al., 2022; Mancoo et al., 2020), we generalize these networks here, removing several previous limitations of SCNs, while retaining their desirable properties.

We first show that a population of spiking neurons with rank-1 connectivity induces a latent variable readout and a spiking boundary in input-output space. We distinguish excitatory (E) and inhibitory (I) boundaries in this space and show that a combined rank-2 EI network is a universal function approximator of static input-output transformations. Importantly, as we demonstrate, flexible function approximation can already be achieved in small networks and can be theoretically understood without invoking large-network-size limits (e.g., mean field). Next, we illustrate a fundamental link between noise suppression and irregular firing and demonstrate how mistuned connectivity between the two populations in the rank-2 EI network can lead to amplification of noise. We then consider effects of slower synaptic dynamics and transmission delays on coding performance, and finally, we show when low-rank spiking networks can be approximated by equivalent low-rank rate networks.

While we limit ourselves here to static function approximation in rank-1 populations with few neurons, in section 5 we touch upon the implications for scaling up this framework to higher-rank networks with richer dynamical motifs, thereby providing a promising path to understand and construct spiking networks in biologically realistic regimes.

For simplicity, we focus on networks with rank-1 populations and use them to develop the central concepts of this article: latent variables and convex boundaries. Following an introduction of the general rank-1 case, we distinguish inhibitory (I) and excitatory (E) populations. Then, we combine the two populations into a rank-2 EI network and illustrate the resultant dynamics and input-output transformations.

The first term on the right-hand side describes a leak current; the second term a feedforward input, $c(t)$⁠, weighted by synaptic feedforward weights, $Fi$⁠; the third term the recurrent spike trains, $sj(t)$⁠, of presynaptic neurons, each weighted by recurrent synaptic weights, $Wij$ (without sign constraints for now); and the last term a constant background input. The constant background simply shifts the effective threshold of the neuron and can therefore be modeled equivalently by setting $bi=0$ and assuming an effective neural threshold $Ti$⁠. For simplicity, we also set the neuron’s membrane time constant $τ=1$⁠, so that the unit of time corresponds to $τ$ (set to 10 ms for all simulations).

We emphasize that linear, weighted sums of filtered spike trains are a common motif not only in the study of neural networks, but also in the analysis of population recordings (Fusi et al., 2016; Saxena & Cunningham, 2019; Keemink & Machens, 2019; Vyas et al., 2020; Jazayeri & Ostojic, 2021). Indeed, such “linear readouts” are used to extract information from neural recordings, either through explicit linear decoding (Dayan & Abbott, 2005) or the use of linear dimensionality-reduction methods such as principal components analysis or factor analysis (Cunningham & Yu, 2014; Pang et al., 2016). Consequently, we can also view $y$ as a component or mode of population activity.

To gain more intuition for these separate regimes, we start with a single neuron. For this neuron, the equation $V1=F1x+E1y=T1$ defines a line in $(x,y)$-space, that is, in the space of inputs and outputs (see Figure 2a, solid blue). On one side of the line, the neuron is subthreshold, and on the other side, it is suprathreshold. The output $y$ always leaks toward zero in the absence of firing. Once the neuron’s threshold is hit, it fires, and the readout bounces down according to the size of the neuron’s decoder weight to $y←y+D1$⁠. From a biophysical point of view, the neuron’s inhibitory self-connection acts as a reset, and the voltage moves from threshold (⁠$V1=T1$⁠) to a hyperpolarized value.

Now considering a network of several neurons, each threshold will trace a different line in the input-output space (see Figure 2b) and together they will form a single boundary (see Figure 2b, gray dashed line; Figure 2c, solid blue line). We will use subthreshold to denote the region where all neurons’ voltages are below threshold and suprathreshold to denote the region where at least one neuron is above threshold. Since $Ei≥0$⁠, the subthreshold region lies below each neuron’s threshold line. Just as in the single-neuron case, the latent variable, $y$⁠, leaks toward zero until it reaches the boundary, where a neuron will fire a spike, causing $y$ to jump back into the subthreshold region, and the dynamics continue. Note that this jump into the subthreshold region pushes $y$ away from all threshold boundaries and therefore corresponds to recurrent inhibition between all neurons. Since each part of the boundary is covered by a single neuron, only one neuron will fire spikes for each value of the input $x$ (see Figure 2d).

We emphasize that any integrate-and-fire network following equation 2.1 and assuming the sign and rank constraints, $Wij=EiDj≤0$⁠, will form a stable spike-threshold boundary. However, a rank-1 population with randomly distributed parameters, $(Fi,Ei,Ti)$⁠, will, more often than not, yield a boundary composed of only a subset of neurons (see appendix B for more details). In other words, some neurons’ boundaries will be well within the suprathreshold set and so will never fire any spikes. Furthermore, a random boundary may also cross the $x$-axis (⁠$y=0$⁠), resulting in some input values for which no neuron is active. Since such parameter regimes include silent neurons or fully silent activity regimes, we do not consider them. Instead, the networks shown in Figures 2b, 2c, and 2e were precisely tuned so that each neuron’s threshold lies along a smooth quadratic function (see appendix C).

We note that due to the positivity constraint, even the diagonal terms of the connectivity matrix, $Wii$⁠, are all positive. For the inhibitory network, in contrast, the diagonal terms were negative and thereby corresponded to the self-reset current after a spike. A positive diagonal term here means that the neuron will self-excite, and thereby immediately fire more spikes. While the absence of a reset term may seem unnatural, it provides a useful entry point for studying low-rank, all-excitatory networks in their idealized form. In practice, a self-reset can be introduced, which we describe below.

The dynamics of an example network are visualized in Figure 3. Once more, each neuron divides the input-output space into two halves, and the boundary of these half spaces is given by equating the voltage equation (see equation 2.10) with the threshold, that is, $Vi=Ti$ (see Figure 3a). The sub- and suprathreshold regions are similarly defined as in the all-inhibitory network and are shown for a single neuron or several neurons in Figures 3a and 3b. The output $y$ is now positive and still leaks toward zero in the subthreshold region, but now this leak is away from the boundary. If the boundary is breached, however, the respective neuron fires, which moves the latent variable further into the suprathreshold regime. The spike causes self-excitation of the firing neuron and, potentially, the crossing of other neurons’ thresholds. (Strictly speaking, this regime is mathematically ill defined. For now, we will be pragmatic, discretize time, and assume that only one neuron can fire per time step. In section 4, we relax this assumption.) In consequence, the dynamics beyond the boundary self-reinforces the growth of the output and becomes highly explosive (see Figure 3c). To illustrate the unstable dynamics on either side of the boundary, the boundary function in Figures 3b and 3c was chosen such that $y=0$ is contained within the subthreshold set for the displayed values of $x$⁠, meaning that in the absence of previous spiking, the network will remain silent. This choice is primarily illustrative, and we return to the spontaneously active case in the next section after reintroducing inhibition.

Before moving on, we make one additional comment about Dale’s law here. The concepts of stable and unstable spike-threshold boundaries form the main computational components of our framework. Intriguingly, sign-constrained inhibitory and excitatory populations are sufficient to generate such well-behaved boundaries, and we will see how the combination of the two leads to a unique computational regime. We note that while a stable boundary can be generated without adhering to Dale’s law, as is known from previous work on spike-coding networks (Calaim et al., 2022), the stability is only local and can collapse upon certain perturbations (see section 5).

Example boundaries for the inhibitory and excitatory populations are shown in Figures 4a and 4b, respectively. The blue boundary describes the inhibitory population and is a negative convex function of $x$ and $yE$⁠, denoted by $yI=-fcvxI(x,yE)$⁠. While we illustrate the boundary as a smooth, continuous surface, in reality it will be piecewise linear, made up of individual neurons’ thresholds as in Figure 2e (see appendix D). Since all decoders are negative, each spike of an inhibitory neuron drives the readout into the subthreshold regime, and the state space is restricted to negative values of $yI$ (Figure 4a). Similarly, the excitatory population is described by the red, negative convex boundary (Figure 4b), which is a function of $x$ and $yI$⁠, denoted $yE=-fcvxE(x,yI)$⁠. However, in this case, each excitatory spike drives the readout into the suprathreshold regime toward positive values of $yE$⁠.

We see that the two boundaries intersect (see Figure 4c), and we can distinguish four regions of the state-space around this intersection. If both populations are below threshold, the population activities will simply leak toward zero (gray area in Figure 4c). If one population is suprathreshold and the other subthreshold (red or blue areas in Figure 4c), then the suprathreshold population activity increases rapidly (through spiking) while the subthreshold activity decreases slowly (through the leak). If both populations are suprathreshold, both activities increase rapidly (purple areas in Figure 4c). In order to eliminate any ambiguities in the suprathreshold regime and to promote stability, we will prescribe that inhibitory neurons always fire before excitatory neurons when both are above threshold. This somewhat ad hoc rule implies that inhibition generally has faster dynamics than excitation and will be relaxed later in section 4.

To illustrate the stable EI dynamics, we simulate the two-neuron system in Figures 4d and 4e. Here, we have a single excitatory and a single inhibitory neuron. When the two outputs $yI$ and $yE$ leak toward zero, both neurons experience a decrease in recurrent inhibition and recurrent excitation, resulting in an overall depolarization of both neurons. The excitatory neuron fires first, and the output $yE$ moves to the right. At this point, both neurons experience an excitatory postsynaptic current. The inhibitory neuron then crosses threshold and, according to our rule above, fires before the excitatory neuron can fire another spike. The inhibitory spike inhibits the excitatory neuron and moves the output $(yE,|yI|)$ into the subthreshold set of both neurons. At this point, the dynamics repeat. While globally stable, the local dynamics of the limit cycle and repeatability of the spiking pattern not only depend on the shapes of the boundaries, but also the decoding weights (see appendix D). Principally the same pattern is observed when we consider a network with several neurons in each population (see Figure 4f). In the deterministic regime that we consider in this section, this oscillatory pattern only encompasses the two neurons that make up the crossing point for the given input $x$ (and local stability can still be assessed with equation 2.41), but the latent dynamics are reflected in all neurons’ voltages (not shown).

We note that in the general case, with multiple neurons per population, the E and I boundaries can in principle cross multiple times, leading to more interesting dynamics such as bistability (not shown, but see section 5).

As shown above, provided that stability conditions are satisfied, the dynamics of the rank-2 EI network will converge toward an oscillation around the crossing point of the two boundaries. Crucially, this crossing point depends on the input $x$⁠. If we add this third dimension back to the picture, we see that the crossing point between the two populations can vary if different neurons form the boundary at different values of $x$ (see Figure 5a). It turns out that such an approach can yield a rich set of possible input-output functions. Indeed, while the above equations cannot generally be solved for $yI$ or $yE$⁠, we can show that each latent variable on its own can, in principle, be any function of the input $x$⁠.

A toy example of function approximation is shown in Figures 5b to 5d, in which a nonconvex, saw-like function is approximated by the excitatory latent readout of a small network with $NE=3$ and $NI=3$ neurons. Such a network computation is intuitively straightforward to visualize and to construct (see appendix E). Following the notation above, we can simply design piecewise-linear convex functions $p(x)$ and $q(x)$ such that their difference approximates (or equals, in this case) the desired function (see Figure 5b). We note that more continuous functions can be approximated with more neurons, and similar nonconvex approximations can also be done with the inhibitory readout (not shown).

While the above derivation was done for a one-dimensional input $x$⁠, it also holds for $M$-dimensional inputs, so that the rank-2 EI network can approximate arbitrary mappings from $RM$ to $R$⁠. For simplicity, we limit ourselves here to rank-2 networks and 1D outputs, which also results in stereotypical, nonoverlapping spike trains (see Figure 5d). We contend that this is largely due to the simplicity of the 1D task, which should be remedied when considering higher-rank networks (see section 5).

Thus far, we have studied deterministic networks, and we have limited ourselves to rank-2 EI networks in which both populations share the same decoders (see equations 2.29 to 2.32; that is, excitatory neurons see the same inhibitory latent variable $yI$ as the inhibitory neurons themselves (and the same for the excitatory latent, $yE$⁠). Here, we will show how relaxing this constraint and adding small amounts of noise to each neuron leads to networks that can control or amplify input-independent noise.

The key consequence of the jittery boundary is that more than one neuron can fire for a given input $x$⁠. In the fully deterministic, inhibitory network, we saw that only one neuron becomes active for any given input $x$ (compare Figure 2d). However, once we add noise to the network, we observe instead that several neurons may become active for any input value $x$⁠, even though the latent variable is still faithfully represented (Figures 6b and 6c).

Geometrically, the individual thresholds start fluctuating around their default position, and if a subset of neural thresholds is close to each other, then the noise and self-reset can push any of these thresholds to the front. As a consequence, neurons within this subset take turns in enforcing the boundary for any particular input value $x$⁠. When the output $y$ decays toward the boundary, it reaches the threshold of the neuron that just happens to be more sensitive at that time point.

This type of redundancy can lead to neurons firing in a seemingly random fashion, simply because the latent variable $y$ can reach one of several thresholds, depending on the noise in the system. Figures 6d and 6e show two trials of the inhibitory network. Just as in previous SCN work (Boerlin et al., 2013; Denève & Machens, 2016; Calaim et al., 2022), we see that the spike trains change from trial to trial. Firing patterns become more irregular if either the number of neurons in the network or the amplitude of the injected noise increases, until they become close to Poisson. However, even in this limit, the amount of injected noise needed is relatively small compared to the synaptic inputs, and so the spiking irregularity does not simply reflect large amounts of external noise (see appendix F). Strictly speaking, the input-output function for such a network is no longer deterministic, but should rather be described by a distribution $p(y|x)$⁠, with a mean that is prescribed by a convex function from $x$ to $y$⁠.

The dynamics of the network can similarly be understood in this space. The firing rate vector, $r$⁠, will leak toward zero in the absence of spiking (compare equation 2.4). When it hits one of the thresholds, the respective neuron, say, neuron $k$⁠, fires, and updates the $k$th firing rate. In the absence of noise, the firing rate vector will always hit the same neuron (see Figure 7a), eventually reaching a fixed point (see Figures 7b and 7c). However, when the neurons’ voltages are contaminated by noise, the thresholds randomly move around their default positions, and other neurons can eventually be the ones that fire (see Figure 7d). As a consequence, the firing rate vector starts to diffuse along the surfaces spanned by the neural thresholds (see Figure 7e).

As a consequence of equation 3.4 the latent readout $y$ is exactly orthogonal to this random diffusion along the threshold boundaries and is therefore largely unaffected by it (see Figures 7b, 7e, and 7f). However, if we project the neural activities onto a readout positioned in the orthogonal subspace, then we can see the large fluctuations from the random diffusion, which is caused by the random, irregular spiking or the neurons (see Figure 7f, gray curve). Accordingly, the recurrent connectivity of the network ensures that the latent readout remains relatively independent of the noise generated by irregular activity, which will only appear in orthogonal directions. This result is similar to previous findings (Boerlin et al., 2013; Landau & Sompolinsky, 2021; DePasquale et al., 2023) and may also relate to the controllability of network dynamics (Kao & Hennequin, 2019).

We now return to the coupled excitatory-inhibitory network, where we have one latent variable for the excitatory and one for the inhibitory population, $yE$ and $yI$⁠, respectively. We first consider a noisy version of the simple two-neuron setup from Figure 4d, replacing each neuron by a population with $NE=50$ and $NI=50$ identical neurons, and adding a small amount of voltage noise (and self-reset) to each neuron, following equation 3.3. The boundary crossing in latent space remains unchanged (see Figure 8a), but is composed of two homogeneous populations. We then add one additional excitatory neuron (labeled $E2$⁠; see Figure 8a) with a different boundary; this neuron does not contribute to coding but is illustrative of how boundaries geometrically relate to voltage, as we discuss below. Finally, to illustrate the inhibitory stabilization of the network, we consider the effect of a positive current stimulation to all neurons of the inhibitory population, a protocol that typically results in a paradoxical effect in inhibition-stabilized networks (Sadeh & Clopath, 2021). Such a stimulation shifts the inhibitory boundary upward and should result in a new steady state of the latent variable activity in which both populations have reduced activity (see Figure 8b).

We simulate this larger, noisier network, applying the perturbation to the inhibitory population for the second half of the simulation time (see Figure 8c). We observe stable coding of the two latent variable readouts as predicted by the boundary crossing, with asynchronous irregular spiking activity. Furthermore, the stimulation of the inhibitory population causes a decrease in the activity of both populations, confirming the network as inhibition-stabilized. We note again that this relies on our assumption that the inhibitory population is faster than the excitatory population (see sections 4 and 5). Next, by separating the positive and negative contributions to the voltage, we see that the synaptic inputs are roughly balanced, such that neurons in the network fire irregularly due to positive input fluctuations (see Figure 8d).

Finally, we note that the additional excitatory neuron, $E2$⁠, does not fire any spikes; its voltage remains hyperpolarized (see Figure 8d, bottom), with the relative amount of hyperpolarization being roughly proportional to the distance between its threshold and the current latent dynamics (see Figures 8a and 8b). We thus see that neurons in the network will generally be in one of two dynamic states: they will either be contributing to coding by firing spikes, in which case their voltages will be balanced, or they will be silent with hyperpolarized voltages (see section 5).

Up to now we have studied an idealized spiking network model that ignores several basic properties of real neurons. Most notably, we have assumed that synaptic input currents are instantaneous; they arrive immediately and are infinitely short, and we have imposed that inhibitory neurons fire before excitatory neurons. As we will see, relaxing these assumptions will influence how sharp the boundary is and how accurately the input is mapped onto the output. In doing so, we arrive at a relationship between spiking dynamics and the smoothed, trial-averaged firing rate dynamics that are typically studied in rate networks.

We now consider a section of the boundary formed by one neuron, say, neuron $k$⁠, and visualize the dynamics in $(y,sk)$-space for the square pulse synapse model (see Figure 9b; shown for three values of $τs$—black, gray, and silver). The system relaxes into a steady-state oscillation at the boundary. However, the presence of the synaptic current, $sk$⁠, complicates the dynamical picture by adding another dimension, so that each neuron now depends on two dynamical variables.

As illustrated in Figure 9d, the infinitely steep boundary has become a finite-sized boundary. We can also see that the size of this boundary grows with the size of the decoder, $Dk$⁠, and is inversely related to the timescale of the synaptic dynamics, $τs$⁠. Moreover, when we add a second, identical neuron to the picture, both neurons will fire two spikes per oscillation cycle, so that the effective size of the boundary doubles (not shown). More generally, the size of the boundary therefore scales linearly with the number of neurons at the boundary as long as their thresholds are sufficiently closely spaced.

Let us next see how the boundary is affected when we add noise to the system. In section 3, we saw that the main effect of voltage or current noise was to randomly move the thresholds of the neurons. Similar effects can be obtained when synaptic weights or voltage resets slightly deviate from their idealized set points (Calaim et al., 2022).

To illustrate this softer boundary at the population level, we simulated the network from Figure 9f with the same setup from Figure 6, but now with finite square-pulse synaptic dynamics (⁠$τs=0.5τ$⁠, see Figures 10a to 10c). We see that in this case, the latent readout still closely follows the threshold boundary (see Figure 10a), and the irregular firing and trial-to-trial variability are still retained (see Figures 10b and 10c). This demonstrates that the finite, soft boundary can still be geometrically visualized in the same way as the idealistic, infinite boundary from before and that fragility issues previously associated with spike-coding networks (Chalk et al., 2016; Rullán Buxó & Pillow, 2020; Calaim et al., 2022) can be alleviated in this regime.

Formally, we have therefore shown that the dynamics of the latent variables of a rank-1 (inhibitory) firing rate network are equivalent to the trial-averaged dynamics of a rank-1 (inhibitory) spiking network, as long as parameters are set such that the networks generate a clear boundary (see sections 4.6 and 5).

We now return to the rank-2 E-I network and consider the effect of slower synapses here. First, we adopt the arguments made above for the stable inhibitory boundary and apply them to the unstable excitatory boundary. We then see that the two boundaries differ by a simple sign flip in the $(|y¯|,y¯˙)$-plot (compare (Figures 11a to 11b). Note that the output of the inhibitory population is negative in our convention (compare Figure 9c), and we here plot its absolute value to allow better comparison with classical plots of E-I networks.

We illustrate these soft boundaries in the two-dimensional latent space, $(|y¯I|,y¯E)$⁠, for two neurons (or two homogeneous populations) in Figures 11c and 11d. At a height of $(|dy¯I/dt|,dy¯E/dt)=0$⁠, we retrieve the nullclines of the dynamics, as illustrated by black lines in Figures 11c and 11d. We observe that the exact shapes and positions of the nullclines depend on the slope parameter $β$ of the sigmoids, which approach the spiking boundaries (see Figure 11e, dotted lines) for large $β$⁠.

For heterogeneous networks, we previously approximated each population’s boundary as lying along a smooth convex curve, similar to Figure 4f. When we soften the boundary, its overall shape is determined by both the individual neural thresholds and the overall smoothness of the sigmoidal boundary (see Figure 11f).

We illustrate the functionality of the rank-2 EI network with finite synapses by repeating the function approximation example from Figures 5b and 5d. As illustrated in Figure 12a, function approximation is made noisier with slowed-down synapses, but it is still reliable. While we no longer apply the ad hoc rule that inhibitory neurons fire first, the inhibitory synaptic dynamics should ideally be faster than the excitatory dynamics in order to ensure stability (in this case, $τE=τ/4$ and $τI=τ/10$⁠; see section 5).

Just as in the previous section, we can similarly derive the above equations by starting with a rank-2 firing rate network. When simulating such a network, we see that the latent readouts smoothly follow the intersection of the two boundaries, and rates will approximate the trial-averaged activities of the spiking network (see Figure 12b).

Beyond making the boundary softer, the combination of mistuned network parameters and realistic synaptic currents can also influence the spiking dynamics and thereby the accuracy of the network’s input-output mapping. Recall that the output in the original network fluctuates between the threshold boundary and the reset boundary (see Figure 2c). However, when synaptic currents are not instantaneous and the thresholds are not fixed, then the output can sometimes cross multiple thresholds, and multiple neurons fire. In this case, the output $y$ does not simply bounce back to the reset boundary but moves further into the subthreshold regime. This will be especially true if the synaptic time constants are very slow, or if there are significant delays in axonal conductance. These resulting oscillatory effects in $y$ decrease the accuracy with which the input is mapped onto the output. An example is shown in Figure 12c, in which $τsE=τsI=τ/4$⁠. Such deviations will also disrupt the efficiency of the code in terms of the number of spikes needed to represent a particular output level (Denève & Machens, 2016) (see section 5).

From a computational perspective, this loss of accuracy and efficiency is not desirable. These problems will likely be mitigated in higher-rank networks (Calaim et al., 2022) and thus may be less severe in scaled-up systems (see section 5). From a biological perspective, simultaneous crossing of multiple thresholds causes synchronous firing, which may relate to oscillatory dynamics in the brain (Chalk et al., 2016), suggesting that some residual effects may be present in real circuits.

We note that the above derivations were carried out for networks that have a block-wise rank-1 structure in which both self-connectivity and cross-connectivity of the excitatory and inhibitory subnetworks are rank-1. However, not all networks with this structure will permit a transition from spiking to rate networks as described here.

First, we have assumed that the cross-connections have the same decoders as the self-connections of each population (compare equations 2.29 and 2.32). As explained in section 3, when the cross-connections are not aligned, the spiking network will accumulate noise, and a match between trial-averaged spiking networks and rate networks can no longer be guaranteed. This restriction on cross-connection alignment is usually not followed in the design of firing rate networks.

Second, parameters in classical firing rate networks are usually chosen such that individual neurons operate in the full regime of their input function, from below threshold to large firing rates and even into saturation. The underlying reason is always that the neural activation functions serve as the basis of an underlying function approximation scheme (compare Figures 1a to 1c). Apart from feedforward networks, examples include line attractor networks in which neural thresholds are staggered in latent space so that the addition of all neural activation functions precisely counters the leak (Seung, 1996; Eliasmith & Anderson, 2003), or decision-making models that use saturation to obtain bistability (Mastrogiuseppe & Ostojic, 2018).

In contrast, in the spiking networks here, neural thresholds are aligned in latent space such that they join to form a boundary. As a consequence, each neuron touches the boundary at some point, and its activation function enforces the boundary. Indeed, parameters need to be chosen such that the boundary is so steep as to be basically unsurmountable by the latent variable. The boundary becomes particularly steep when the synaptic currents, $τs$⁠, become sufficiently short. Apart from that constraint, the precise shape of the activation function does not matter.

Importantly, this operating regime is crucial in order to obtain balanced synaptic current inputs into the neurons (compare Figures 1d and 1e). In the rate domain, our perspective bears a resemblance to other rate-based models featuring nonsaturating nonlinearities (Miller & Troyer, 2002; Hansel & van Vreeswijk, 2002; Ahmadian et al., 2013; Rubin et al., 2015; Hennequin et al., 2018). Indeed, we would reinterpret the respective neural activation functions as implementing soft boundaries.

In this article, we have proposed a new framework for constructing and interpreting a broad class of spiking neural network models based on spike-threshold boundaries. These boundaries come from the inherent inequality at the core of each neuron—the voltage threshold—which triggers a spike. Classic approaches usually seek to eliminate the threshold inequality by, for example, filtering or averaging spike trains, and to replace it with a continuous activation function (Dayan & Abbott, 2005; Gerstner et al., 2014). We take a different approach here. By defining a one-dimensional latent variable readout, we show that the spike-threshold boundaries of rank-1 networks can be visualized in input-output space, with each neuron’s threshold forming an affine boundary. A population of neurons delineates a convex boundary separating subthreshold from suprathreshold areas and for an inhibitory population, this boundary creates a stable input-output function. For an excitatory population, the boundary is instead unstable but can then be stabilized in a coupled EI network, where the dynamics sit at the crossing of the two boundaries. The boundary between sub- and supra-threshold regimes thus becomes the center of spike-based computation, and distinct from the “feature detector” perspective (see Figure 1).

The mathematics of systems with inequalities and hard boundaries has a rich history in constrained convex and nonconvex optimization (Boyd & Vandenberghe, 2004). Previous work on spike coding has demonstrated how the dynamics of spiking networks can be mapped onto convex optimization problems such as quadratic programming (Barrett et al., 2013; Mancoo et al., 2020), yielding stable, convex input-output transformations. We extended this work here by showing that coupled EI networks are capable of nonconvex function approximation, with links to difference of convex (DC) programming (Horst & Thoai, 1999; Lipp & Boyd, 2016). Similarly, previous work has linked EI networks to minimax optimization (Seung et al., 1997; Li & Pehlevan, 2020), but with a focus more on attractor dynamics rather than function approximation. Overall, the analogy between spiking thresholds and constraint boundaries offers a fruitful perspective that may lead to new insights and algorithms that harness the power of convex optimization (Boyd & Vandenberghe, 2004). Moreover, such techniques have recently caught on in deep learning in the form of deep implicit layers, offering potential avenues for learning and scaling up our framework (Amos & Kolter, 2017; Gould et al., 2021).

In the work presented here, we limited ourselves to small populations of neurons with rank-1 connectivity. Though increasing the number of neurons in rank-1 (or rank-2 EI) networks may enable closer and closer approximations to arbitrary continuous 1D functions, the $N→∞$ limit of rank-1 networks is likely too limiting to relate to biological spiking networks. Furthermore, though we do not demonstrate it here, this larger-scale limit requires the absence of noise and communication delays, in addition to having very large synaptic input currents. Instead, to stay within a biologically realistic and workable regime, scaling up the number of neurons should coincide with scaling up the rank of the connectivity (Calaim et al., 2022). Fortunately, many of the intuitions about spike thresholds and convex boundaries transfer to the higher-dimensional input-output spaces of arbitrary rank-K networks. Furthermore, many signatures of biological activity regimes discussed here, such as irregular firing and balance, become more robust in this higher-dimensional regime (Calaim et al., 2022).

The rank-1 and rank-2 spiking networks studied here can be seen as a particular generalization of spike-coding networks (SCNs; Boerlin et al., 2013; Denève & Machens, 2016), with the spike-threshold boundary being a generalization of the concept of a bounding box put forward in Calaim et al. (2022). As we demonstrate, this generalization frees SCNs from the constraint of linear computation without evoking dendritic or synaptic nonlinearities (Thalmeier et al., 2016; Alemi et al., 2018; Nardin et al., 2021; Mikulasch et al., 2021). We can reinterpret the bounding box as being composed of two components: local regions of the boundary are stable, with lateral inhibition between neurons, but globally the closed nature of the box may lead to errant positive feedback and thus fragility in the form of epileptic (“ping-pong”) behavior (Chalk et al., 2016; Rullán Buxó & Pillow, 2020; Calaim et al., 2022). Our framework eliminates this previous limitation by separating stable and unstable boundaries into two populations in the network. Despite these differences, many of the desirable biological properties and predictions of SCNs are preserved here, including robustness to cell death and inhibition, firing irregularity, E/I balance, and coding efficiency. While it is out of the scope of this current study to explore these properties in detail, in future work it will be interesting to determine how the EI nature of this model leads to unique biological predictions, especially in higher-dimensional networks. More generally, the spike-coding framework has been extended in various ways for other computational aims (Slijkhuis et al., forthcoming; Masset et al., 2022), as well as to incorporate additional biologically plausible details or activity regimes (Schwemmer et al., 2015; Koren & Panzeri, 2022; Safavi et al., 2023), and it will be interesting to see how our perspective here may be integrated with these other studies. Taking into account slower synapses with delays, we are able to relate our work to the dynamics of rate networks, similar to previous approaches (Kadmon et al., 2020; Rullán Buxó & Pillow, 2020). However, our aim here is not to emulate a rate-based system with spikes but rather to demonstrate the relationship between hard spike boundaries and soft rate boundaries.

One major biological implication of our work is the natural distinction between excitatory and inhibitory populations. Many previous models, in contrast, typically start with Dale’s law as a biological constraint, then attempt to achieve performance on par with an unconstrained network, and finally reflect on the potential benefits of such a constraint for, say, stability or robustness (Hennequin et al., 2014; Song et al., 2016; Brendel et al., 2020; Cornford et al., 2021; Haber & Schneidman, 2022). Our framework instead puts forth a particular computational role for Dale’s law: it is the combination of positive feedback from the unstable excitatory boundary with negative feedback from the stable inhibitory boundary that enables flexible input-output mappings. This idea is reminiscent of so-called mixed-feedback systems in control theory, which are beginning to be explored in the neuroscience domain (Sepulchre et al., 2019; Sepulchre, 2022). In addition to the unique computational role of Dale’s law in our model, the separation into excitatory and inhibitory populations also links the dynamical behavior of our model to previous work on EI networks, including inhibition-stabilized networks (Sadeh & Clopath, 2021).

A strict requirement of the rank-2 EI networks that we study is that inhibition should be faster than excitation in order to keep the dynamics tightly around the intersection of the two boundaries. While some evidence suggests that the synaptic time constant of inhibition may be slower than that of excitation (e.g., see the discussion in Fourcaud & Brunel, 2002), other factors, such as the localization of inhibitory synapses close to the peri-somatic zone (Freund & Katona, 2007), or cortical connectivity patterns (short-range for inhibition, long-range for excitation) suggest that inhibition could also be faster than excitation. In addition, self-stabilizing mechanisms within the excitatory population, such as spike-frequency adaptation, could make the speed of inhibition less important.

Overall, our work can be seen as a bridge between mechanistic and functional network models of the brain. Mechanistic models, perhaps best represented by balanced spiking networks (Van Vreeswijk & Sompolinsky, 1996; Amit & Brunel, 1997; Brunel, 2000), capture many important features of cortical dynamics, including asynchronous irregular activity, excitation-inhibition balance, and correlated neural variability (Vogels & Abbott, 2005; Renart et al., 2010; Rosenbaum et al., 2017; Huang et al., 2019). However, such models, and especially their theoretical analysis, have typically been limited to linear computations in large networks (infinitely many neurons in the mean-field limit, or 1000 to 10,000 in practice). In contrast, our work focuses on nonlinear computations in much smaller spiking networks (5 to 50 neurons).

While various frameworks have been proposed to obtain nonlinear computations in the mean-field limit (Renart et al., 2007; Roudi & Latham, 2007; Hansel & Mato, 2013; Lajoie et al., 2016; Ingrosso & Abbott, 2019; Kim & Chow, 2021), one recent study parallels our work in several interesting ways. Baker et al. (2020) utilize networks with “semibalance” or “approximate balance” in order to achieve nonlinear computation. Notably, they observe that active subpopulations retain balanced input, whereas others are hyperpolarized. Despite the size differences, this is precisely what our networks predict as well (see Figure 8d), and future work should explore links between this work and our proposed framework. More generally, this supports the notion that the average balance must be broken or loose enough to enable nonlinear computations (Ahmadian & Miller, 2021).

Functional models of static input-output mappings are perhaps best represented by feedforward rate networks (Hunsberger & Eliasmith, 2015; Yamins & DiCarlo, 2016). However, such networks do not agree with the heavily recurrent connections found in the brain, they do not obey Dale’s law, and they rely on neurons being feature detectors, which conflicts with the idea of balanced inputs (see Figure 1). Here we have shown that function approximation can also be achieved in (low-rank) recurrent networks obeying Dale’s law, and with balanced inputs. Typically, the study of low-rank networks has emphasized the generation of internal dynamics for decision making, memory, and motor control (Eliasmith, 2005; Mastrogiuseppe & Ostojic, 2018; Dubreuil et al., 2022). Though we have deliberately ignored dynamics here, some simple dynamical motifs such as bistability can be achieved even in the rank-2 EI network provided that there are multiple boundary crossings (not explored here). Higher-dimensional systems can have much richer dynamics in the latent space, as is already known from work on rate networks, and will be explored in future work. One aspect of dynamics that we do consider in this work is the existence of noise-amplifying or noise-suppressing connectivities. These properties have also been described in other low-rank networks (Landau & Sompolinsky, 2018, 2021; DePasquale et al., 2023).

Advances in learning algorithms have led to an explosion of research on training functional spiking network models (Abbott et al., 2016; Neftci et al., 2019). Interestingly, many algorithms smooth over the spike-threshold nonlinearity using “surrogate” or pseudo-gradients (Bellec et al., 2018; Neftci et al., 2019; Zenke & Vogels, 2021), or directly train on rate dynamics (DePasquale et al., 2023), which may smooth out or prevent the development of well-formed latent boundaries. It will be interesting to develop scalable training algorithms for the networks discussed here and to see how additional constraints during training may affect the development of spike-threshold boundaries.

In summary, this article puts forth a new perspective on spike-based computation. It offers the potential to better understand the computational regimes and benefits of the balanced state as observed in cortex, and presents a theoretical framework to construct functional spiking neural networks based on these insights. Finally, while the potential to scale up the framework to higher-rank networks with dynamical computations is still to be demonstrated, it may serve as a promising direction for scalable computation with spikes.