## Abstract

Deep feedforward and recurrent neural networks have become successful functional models of the brain, but they neglect obvious biological details such as spikes and Dale’s law. Here we argue that these details are crucial in order to understand how real neural circuits operate. Towards this aim, we put forth a new framework for spike-based computation in low-rank excitatory-inhibitory spiking networks. By considering populations with rank-1 connectivity, we cast each neuron’s spiking threshold as a boundary in a low-dimensional input-output space. We then show how the combined thresholds of a population of inhibitory neurons form a stable boundary in this space, and those of a population of excitatory neurons form an unstable boundary. Combining the two boundaries results in a rank-2 excitatory-inhibitory (EI) network with inhibition-stabilized dynamics at the intersection of the two boundaries. The computation of the resulting networks can be understood as the difference of two convex functions and is thereby capable of approximating arbitrary non-linear input-output mappings. We demonstrate several properties of these networks, including noise suppression and amplification, irregular activity and synaptic balance, as well as how they relate to rate network dynamics in the limit that the boundary becomes soft. Finally, while our work focuses on small networks (5-50 neurons), we discuss potential avenues for scaling up to much larger networks. Overall, our work proposes a new perspective on spiking networks that may serve as a starting point for a mechanistic understanding of biological spike-based computation.

## 1 Introduction

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.

## 2 Spiking Thresholds Form Convex Boundaries in Latent Space

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.

### 2.1 Rank-1 Connectivity Generates a Latent Variable

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 $\tau =1$, so that the unit of time corresponds to $\tau $ (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.

### 2.2 Inhibitory Neurons Form Stable, Attracting Boundaries

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\u2190y+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\u22650$, 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\u22640$, 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).

### 2.3 The Boundary Determines the Input-Output Mapping

### 2.4 The Boundary Determines the Latent Dynamics

### 2.5 Excitatory Neurons Form Unstable, Repellent Boundaries

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).

### 2.6 The Inhibitory Boundary Can Stabilize the Excitatory Boundary

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).

### 2.7 The Rank-2 EI Network Can Approximate Arbitrary, Nonlinear Input-Output Functions

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).

## 3 Synaptic Connections Can Suppress or Amplify Noise

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.

### 3.1 Noise Makes Boundaries Jitter

### 3.2 A Jittery Boundary Causes Irregular Spiking

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$.

### 3.3 Irregular Firing Amplifies Noise in the Decoder Nullspace

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).

### 3.4 The Rank-2 EI Network is Inhibition-Stabilized and Balanced

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).

### 3.5 Mistuned Cross-Connections Amplify Noise Due to Irregular Firing

## 4 Rate Networks Can Approximate Spiking Networks in the Latent Space

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.

### 4.1 Slower Synapses Generate a Finite Boundary

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 $\tau 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, $\tau 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.

### 4.2 Noise Causes a Soft Boundary

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 ($\tau s=0.5\tau $, 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.

### 4.3 The Latent Variables of Spiking and Rate Networks Are Equivalent

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).

### 4.4 The E-I Boundaries Become Nullclines in the Firing-Rate Limit

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\xaf|,y\xaf\u02d9)$-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\xafI|,y\xafE)$, for two neurons (or two homogeneous populations) in Figures 11c and 11d. At a height of $(|dy\xafI/dt|,dy\xafE/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 $\beta $ of the sigmoids, which approach the spiking boundaries (see Figure 11e, dotted lines) for large $\beta $.

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, $\tau E=\tau /4$ and $\tau I=\tau /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).

### 4.5 Synaptic Timescales Impact the Accuracy of the Output

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 $\tau sE=\tau sI=\tau /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.

### 4.6 The Spike-Rate Equivalence Requires a Latent Boundary

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, $\tau 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.

## 5 Discussion

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\u2192\u221e$ 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.

## Appendix A: Reference Table of Model Parameters

Variable(s) . | Description ($X,Y\u2208{E,I}$ indicate population identity) . |
---|---|

$Vi(t),ViX(t)$ | Voltage of neuron $i$ (single-pop. or for pop. $X$) |

$Ti,TiX$ | Threshold of neuron $i$ (single-pop. or for pop. $X$) |

$\tau $ | Membrane time constant (and simulation timescale) |

$\tau s$ | Synaptic time constant |

$\mu $ | Single-neuron cost (“soft” refractory period) |

$c(t)$ | Network input |

$x(t)$ | Filtered network input |

$s(t),sX(t)$ | Spike trains (single pop. or for pop. $X$) |

$r(t),rX(t)$ | Filtered spike trains (single pop. or for pop. $X$ |

$Ei,EiXY$ | Encoder weights (single-pop. or from pop. $Y\u2192X$) |

$Dj,DjX$ | Decoder weights (single-pop. or for pop. $X$) |

$DX(x)$ | Decoder function (single-pop. or for pop. $X$) |

$Wij,WijXY$ | Recurrent weight from $j\u2192i$ (single-pop. or for $Y\u2192X$) |

$y(t),yX(t)$ | Latent variable output (single-pop. or for population $X$) |

$Fi,FiX$ | Feedforward (input) weight (single-pop. or for pop. $X$) |

$B(x)=fcvx(x)+y$ | Population boundary function |

$I()$ | Indicator function ($I(z)=\u221e$ if $z\u22650$, $I(z)=0$ otherwise) |

$\sigma \beta ()$ | Sigmoid function with slope parameter $\beta $ |

Variable(s) . | Description ($X,Y\u2208{E,I}$ indicate population identity) . |
---|---|

$Vi(t),ViX(t)$ | Voltage of neuron $i$ (single-pop. or for pop. $X$) |

$Ti,TiX$ | Threshold of neuron $i$ (single-pop. or for pop. $X$) |

$\tau $ | Membrane time constant (and simulation timescale) |

$\tau s$ | Synaptic time constant |

$\mu $ | Single-neuron cost (“soft” refractory period) |

$c(t)$ | Network input |

$x(t)$ | Filtered network input |

$s(t),sX(t)$ | Spike trains (single pop. or for pop. $X$) |

$r(t),rX(t)$ | Filtered spike trains (single pop. or for pop. $X$ |

$Ei,EiXY$ | Encoder weights (single-pop. or from pop. $Y\u2192X$) |

$Dj,DjX$ | Decoder weights (single-pop. or for pop. $X$) |

$DX(x)$ | Decoder function (single-pop. or for pop. $X$) |

$Wij,WijXY$ | Recurrent weight from $j\u2192i$ (single-pop. or for $Y\u2192X$) |

$y(t),yX(t)$ | Latent variable output (single-pop. or for population $X$) |

$Fi,FiX$ | Feedforward (input) weight (single-pop. or for pop. $X$) |

$B(x)=fcvx(x)+y$ | Population boundary function |

$I()$ | Indicator function ($I(z)=\u221e$ if $z\u22650$, $I(z)=0$ otherwise) |

$\sigma \beta ()$ | Sigmoid function with slope parameter $\beta $ |

## Appendix B: Networks with Randomly-Distributed Parameters

We emphasized in the main text that while all presented example populations were tuned to have a particular population-level boundary, our framework is fully general. This means that any spiking network with low-rank connectivity and sign constraints (reflecting inhibitory or excitatory populations) will have a well-defined boundary. The mathematical reason is that each neuron’s voltage inequality defines a half-space of subthreshold voltages, and the intersection of all of these half-spaces always generates a convex subthreshold set, which is bordered by a convex boundary.

Here, we demonstrate what happens when the neural parameters are not tuned to form a particular boundary. In Figure 13a, we show a network of $N=8$ inhibitory neurons, where each neuron has randomly chosen values for input weights, $Fi$, latent variable encoding weights, $Ei$, and threshold $Ti$ (note that we are only interested in the shape of the boundary here, so the values of $Di$ are not relevant). Two things are apparent. First, there are fully silent neurons, and actually the population boundary is composed of only three of the eight neurons in this particular input domain. Second, there are values of $x$ for which the network itself is silent ($0<x<0.25$).

To avoid these degenerate regimes, we can simply tune each neuron’s threshold $Ti$ (equivalent to tuning a background input current to each neuron) such that it is tangent to a desired boundary (see Figure 13b, gray dotted line). As shown in Figure 13b, this effectively removes the degeneracy of the fully random network and ensures that all neurons participate in forming the boundary. However, we now have a case in which different areas of the boundary may have different amounts of redundancy and some areas with mismatches to the desired boundary. We thus conclude that (smaller) networks with randomly distributed parameters, though still able to form population boundaries, are not ideal for approximating arbitrary boundaries. Though we used the inhibitory population as an example, these results are applicable to the excitatory population and its respective population boundary as well.

## Appendix C: Network Parameters in the Single Rank-1 Population

## Appendix D: Network Parameters in the Rank-2 EI Network

## Appendix E: Difference of Convex Decomposition

The expressibility of difference of convex functions has been explored in various works, perhaps starting with Hartman (1959), and has a rich history of applications in optimization (Horst & Thoai, 1999; Yuille & Rangarajan, 2003; Bačák & Borwein, 2011; Lipp & Boyd, 2016). Due to its importance to the results presented here, we provide a brief overview of this work.

To begin, we review theorem 1 given in Yuille and Rangarajan (2003), which proves that any function with a bounded second derivative can be written as the difference of two convex functions (or equivalently as the sum of a convex and concave function), and also provides a recipe for this decomposition:

(Yuille & Rangarajan, 2003). Let $f(x)$ be an arbitrary continuous function with a bounded second derivative $d2f(x)/dx2$. Then we can always decompose it into the difference of two convex functions.

We begin by selecting an arbitrary convex function $g(x)$, whose second derivative is bounded below by a positive constant $\epsilon >0$. Then there exists a positive constant $\lambda $ such that the second derivative of $q(x)=f(x)+\lambda g(x)$ is nonnegative. Consequently, $q(x)$ is convex. We then define a second convex function, $p(x)=\lambda g(x)$, and confirm that $f(x)=q(x)-p(x)=f(x)+\lambda g(x)-\lambda g(x)$.

A demonstration of this decomposition is shown in Figures 14a to 14d, in which a smooth, nonconvex function $f(x)$ (see Figure 14a) is decomposed into the difference of $q(x)=f(x)+\lambda x2$ (see Figure 14b, red) and $p(x)=\lambda x2$ (see Figure 14b, blue). As illustrated in Figures 14c and 14d, the derivatives of these two convex functions are monotonic, and their second derivatives are strictly nonnegative, thus confirming them to be convex.

We observe that this theorem and proof apply perfectly well when $f(x)$ is continuous and piecewise-linear. However, the theorem leaves open how to choose $g(x)$ and how to work with the piecewise-linear functions required in the spiking framework. Fortunately, piecewise-linear functions can be uniquely decomposed into two piecewise-linear convex functions (Melzer, 1986; Kripfganz & Schulze, 1987; Siahkamari et al., 2020). To illustrate the main idea, we present an additional proof.

Let $f(x)$ be an arbitrary *piecewise-linear* continuous function. Then we can always decompose it into the difference of two piecewise-linear convex functions.

This second approach is demonstrated in Figures 14e to 14h, in which a piecewise-linear, nonconvex function $f(x)$ (see Figure 14e) is decomposed into the difference of two piecewise-linear, convex functions $q(x)$ (see Figure 14f, red) and $p(x)$ (see Figure 14f, blue). As illustrated in Figure 14g, the piecewise-constant gradient of the original function is decomposed into two piecewise-constant, monotonic functions that account for the positive and negative jumps in the gradient. The second-order subdifferential of $f(x)$ consists of a series of delta functions at the edge points. We now see that the two convex functions contain the same edge points, but they are strictly positive.

## Appendix F: Noise and Irregular Firing in the Inhibitory Network

In section 3, we illustrated how small amounts of input current noise lead to irregular firing, but it was left unspecified how much noise is actually needed. To quantify the relationship between noise and irregularity, we simulated the inhibitory population from Figure 6 for varying noise levels and compared it with an equivalent feedforward population in which recurrent weights were set to zero. Focusing on one particular neuron (see Figure 15a, blue line) and input level (gray vertical line), we simulated many trials of the network and recorded spike trains. Then, computing the spike count coefficient of variation (CV) across trials, we show that the recurrent network architecture presented here results in a substantial amount of irregularity even for small amounts of noise (see Figure 15b). The fact that the variability in the network far exceeds the feedforward network demonstrates that the irregularity is not simply due to input noise, but rather reflects the combination of the small amount of noise with the recurrent architecture. Example rasters over trials are shown for two particular noise levels for the two networks (see Figures 15c and 15d).

## Appendix G: Synaptic Dynamics

In sections 2 and 3, we impose two very strict requirements on the synaptic dynamics: we assume synaptic currents to be instantaneous, and assume that inhibitory neurons always fire before excitatory neurons. The notion of instantaneous communication in a spiking network simulation with discrete time steps is tricky. Technically speaking, even in models without synaptic delays, a delay equivalent to the simulation time step, $dt$, is implicit. Event-based simulators circumvent this problem by only integrating voltages up to the next spike. To keep the simplicity of a discrete-time simulation, we instead impose that only one neuron spikes per time step. For the single inhibitory population, the rationale is straightfoward. When the first inhibitory neuron fires a spike, it will inhibit all other neurons in the network and thereby prevent other neurons from spiking (provided $dt$ is sufficiently small). For the EI network, an excitatory spike may instead drive other excitatory and inhibitory neurons above threshold. Here, we make the important assumption that inhibitory neurons will generally be driven above threshold faster than other excitatory neurons, and thereby are able to prevent additional excitatory spikes (see section 5).

However, we generally keep the negative diagonal terms of the connectivity instantaneously fast, as they correspond to voltage resets.

## Code Availability

Code used to generate all figures of this article is available on GitHub at https://github.com/wpodlaski/funcapprox-with-lowrank-EI-spikes. This repository also contains tutorial notebooks for visualizing latent boundaries and function approximation.

## Acknowledgments

We thank Alfonso Renart, Francesca Mastrogiuseppe, and Sander Keemink for detailed comments on the manuscript. We thank all members of the Machens Lab for helpful discussions and feedback, especially Allan Mancoo for discussions on convex optimization, and Bertrand Lacoste for discussions pertaining to decomposition into the difference of two convex functions. This work was supported by the Fundação para a Ciência e a Tecnologia (project FCT-PTDC/BIA-OUT/32077/2017-IC&DT-LISBOA-01-0145-FEDER) and by the Simons Foundation (Simons Collaboration on the Global Brain #543009).

## References

*Nature Neuroscience*

*Neuron*

*Neural Computation*

*Proceedings of the AAAI Conference on Artificial Intelligence*

*Cerebral Cortex*

*Proceedings of the International Conference on Machine Learning*

*Optimization*

*PLOS Computational Biology*

*Current Opinion in Neurobiology*

*Advances in neural information processing systems*

*eLife*

*Advances in neural information processing systems*

*Journal of Physiology*

*PLOS Computational Biology*

*Convex optimization*

*PLOS Computational Biology*

*Journal of Computational Neuroscience*

*eLife*

*eLife*

*Current Opinion in Neurobiology*

*PLOS Computational Biology*

*Proceedings of the 9th International Conference on Learning Representations*

*Nature Neuroscience*

*Theoretical neuroscience: Computational and mathematical modeling of neural systems*

*Nature Neuroscience*

*Neuron*

*Nature Neuroscience*

*Notes and Records of the Royal Society of London*

*Neural Computation*

*Neural engineering: Computation, representation, and dynamics in neurobiological systems*

*Current Opinion in Neurobiology*

*Neural Computation*

*Neuron*

*Current Opinion in Neurobiology*

*Neuron*

*Neuronal dynamics: From single neurons to networks and models of cognition*

*IEEE Transactions on Pattern Analysis and Machine Intelligence*

*36th Conference on Neural Information Processing Systems*

*Journal of Neuroscience*

*Journal of Neuroscience*

*Journal of Neuroscience*

*Pacific Journal of Mathematics*

*Neuron*

*Neuron*

*Journal of Optimization Theory and Applications*

*Neuron*

*Spiking deep networks with LIF neurons.*

*PLOS One*

*Current Opinion in Neurobiology*

*Advances in neural information processing systems*

*Current Opinion in Neurobiology*

*Current Opinion in Neurobiology*

*Neural Computation*

*Advances in neural information processing systems*

*Optimization*

*PLOS Computational Biology*

*PLOS Computational Biology*

*Physical Review Research*

*Nature Reviews Neuroscience*

*37th Conference on Neural Information Processing Systems*

*Advances in neural information processing systems*

*Journal of Cognitive Neuroscience*

*Optimization and Engineering*

*Advances in neural information processing systems*

*Cerebral Cortex*

*Advances in neural information processing systems*

*Neuron*

*Quasidifferential calculus*

*eLife*

*Proceedings of the National Academy of Sciences*

*Journal of Neurophysiology*

*Peer Community Journal*

*IEEE Signal Processing Magazine*

*Current Biology*

*Neural Computation*

*Science*

*Neural Computation*

*Convex analysis*

*Nature Neuroscience*

*PLOS Computational Biology*

*Neuron*

*PLOS Computational Biology*

*Nature Reviews Neuroscience*

*Signatures of criticality in efficient coding networks*

*Current Opinion in Neurobiology*

*Journal of Neuroscience*

*Proceedings of the IEEE*

*Annual Review of Control, Robotics, and Autonomous Systems*

*Proceedings of the National Academy of Sciences*

*Advances in neural information processing systems*

*Journal of Neuroscience*

*PLOS Computational Biology*

*Proceedings of the International Conference on Machine Learning*

*IEEE Transactions on Cognitive and Developmental Systems*

*Journal of Neuroscience*

*PLOS Computational Biology*

*Current Opinion in Neurobiology*

*PLOS Computational Biology*

*Science*

*Journal of Neuroscience*

*Annual Review of Neuroscience*

*Biophysical Journal*

*Nature Neuroscience*

*Neural Computation*

*Neural Computation*