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.

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 1a1c). 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 1d1f). 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.

Figure 1:

The operating regimes and sign constraints of rate (a–c) and balanced spike (d–f) codes are fundamentally different. (a) A rate-coding perspective utilizes the time-averaged frequency-input (f-I) curve (shown in black for a leaky integrate-and-fire (LIF) neuron receiving noisy input) as a basis; such neurons may be activated far above or below threshold (green ellipse). (b) In the absence of additional noise, rate-coding neurons fire regularly in response to an input stimulus (three trials shown) with a rate according to the f-I curve. (c) Function approximation through linear-nonlinear mappings may be achieved in rate-based spiking networks (Eliasmith & Anderson, 2003); here, the black output curve is composed of the weighted sum of the four spiking neuron f-I curves (n1, green; n2, blue; n3, purple; n4, orange), with each neuron having a positive or negative output weight. (d) A balanced-input perspective of a spiking neuron suggests a more localized operating regime with mean input below or at the neuron’s threshold (green ellipse), and thus does not follow the f-I curve (black); this balance is due to the equal strength of excitatory and inhibitory inputs. (e) This regime explains the irregular firing and trial-to-trial variability seen in cortex. (f) It is less well understood how such a balanced, Daleian regime may explain biological computation.

Figure 1:

The operating regimes and sign constraints of rate (a–c) and balanced spike (d–f) codes are fundamentally different. (a) A rate-coding perspective utilizes the time-averaged frequency-input (f-I) curve (shown in black for a leaky integrate-and-fire (LIF) neuron receiving noisy input) as a basis; such neurons may be activated far above or below threshold (green ellipse). (b) In the absence of additional noise, rate-coding neurons fire regularly in response to an input stimulus (three trials shown) with a rate according to the f-I curve. (c) Function approximation through linear-nonlinear mappings may be achieved in rate-based spiking networks (Eliasmith & Anderson, 2003); here, the black output curve is composed of the weighted sum of the four spiking neuron f-I curves (n1, green; n2, blue; n3, purple; n4, orange), with each neuron having a positive or negative output weight. (d) A balanced-input perspective of a spiking neuron suggests a more localized operating regime with mean input below or at the neuron’s threshold (green ellipse), and thus does not follow the f-I curve (black); this balance is due to the equal strength of excitatory and inhibitory inputs. (e) This regime explains the irregular firing and trial-to-trial variability seen in cortex. (f) It is less well understood how such a balanced, Daleian regime may explain biological computation.

Close modal

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.

2.1  Rank-1 Connectivity Generates a Latent Variable

Let us consider a network of N recurrently connected spiking neurons. We denote each neuron’s membrane potential by Vi and assume that the neuron fires a spike when it reaches a threshold, T. Each neuron’s membrane potential will be described as a leaky integrate-and-fire neuron:
τV˙i(t)=-Vi(t)+Fic(t)+j=1NWijsj(t)+bi.
(2.1)

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 model each neuron’s spike train as a series of delta functions, such that si(t)=fδ(t-tif), where tif is the time of the fth spike of the ith neuron. For now, we assume instantaneous communication between neurons without temporal synaptic dynamics or delays (this will be relaxed in section 4). We furthermore assume that the diagonal of the weight matrix, Wii, includes the voltage reset after a spike, thereby setting the reset voltage to
Vi,reset=Ti-Wii.
(2.2)
For our exposition, we will mostly rely on the integrated version of these equations (Gerstner et al., 2014). We first define the filtered input, x(t), and the filtered spike trains, ri(t), as
x˙(t)=-x(t)+c(t),
(2.3)
r˙i(t)=-ri(t)+si(t).
(2.4)
These equations express a filtering or convolution of the input or spike trains, respectively, with a one-sided exponential kernel, h(t)=H(t)exp(-t), where H(t) is the Heaviside function (i.e., H(t)=1 if t0, and H(t)=0 otherwise). Accordingly, we can think of ri(t) as a simple model of a neuron’s postsynaptic potential or as a (single-trial) estimate of a neuron’s instantaneous firing rate.
With these definitions, we can integrate equation 2.1 to obtain
Vi=Fix+j=1NWijrj,
(2.5)
where we have dropped the time index for brevity. By taking the derivative of this equation, inserting the equations for the filtered inputs and spike trains, equations 2.3 and 2.4, and remembering that Ti=T-bi, one retrieves equation 2.1. A key assumption for this section will be that the weight matrix has rank-1, so that we can write
Wij=EiDj,
(2.6)
where we will call the scalars Ei the encoding weights, and Dj the decoding weights. With this definition, the integrated voltage, equation 2.5, becomes
Vi=Fix+j=1NEiDjrj
(2.7)
=Fix+Eij=1NDjrj.
(2.8)
The term inside the parentheses is simply a linear combination of the filtered spike trains, independent of index i. We denote it as
y=j=1NDjrj,
(2.9)
so that the voltage becomes
Vi=Fix+Eiy.
(2.10)
Importantly, the variable y fully controls the dynamics of the network, in that knowledge of y (together with the input x) is sufficient to compute the voltages and, consequently, the spike trains. We will refer to y as the readout, the output, or the latent variable of the network.

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

So far we have not put any sign constraints on the connectivity. Now, considering a population of inhibitory neurons with negative weights, we impose
Wij=EiDj0i,j
(2.11)
by requiring
Ei0i,
(2.12)
Dj0j.
(2.13)
With this sign convention, the negative sign is associated with the decoding weights, and from equation 2.9, the latent variable or readout will also be constrained to be negative, y0. The voltage equation, equation 2.10, remains unchanged, but the positive encoding weights Ei and negative latent variable y ensure that each neuron receives negative feedback from the population.
Our first goal will be to understand the dynamics of the integrate-and-fire network in the joint space of inputs and readouts. Using equations 2.4 and 2.9, we can compute the derivative of the readout,
y˙=-y+j=1NDjsj,
(2.14)
revealing the dynamics to be a weighted, leaky integration of the network’s spike trains. These dynamics therefore fall into two regimes:
y˙=-y,yy+DiifVi=Fix+EiyTi.
(2.15)
Either the readout is leaking toward zero (top equation) or it is changing abruptly due to spikes (bottom equation).

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 yy+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.

Figure 2:

A rank-1 inhibitory population forms a stable, attracting boundary. (a) A single inhibitory neuron’s voltage (see equation 2.10) is visualized in input-output space; the spike threshold (V1=T1; solid blue line, with normal vector (F1,E1)) divides the space into sub- and suprathreshold sets (white and blue shaded area). Gray arrows illustrate dynamics due to the leak; the downward-facing blue arrow illustrates the effect of a spike at the boundary, which sets yy+D1 (dotted blue line). (b) The thresholds of multiple (N=10; colored lines) neurons were tuned to delineate a concave (or negative convex) boundary, y=-fcvx(x)=-x2-12 (dotted gray line). (c) The boundary (blue line) divides the input-output space into a convex subthreshold set (with all voltages below threshold) and a suprathreshold set (with at least one voltage above threshold), forming a stable input-output relationship. Output precision is defined as the distance between the neural thresholds and the resets after a spike, which can be different for each neuron (gray shading; see equation 2.18). (d) Latent variable output (top) and spike raster (bottom) of a simulation of the spiking network from panel c (colors as in panel b) for a time-varying input from x=-1 to 1 (light gray, top; time in units of τ); the spiking simulation output (solid dark gray) follows the true boundary (dashed black line), and each input value is coded by a single neuron. (e) A negative convex surface boundary for a two-dimensional input (see equation 2.16), made up of 36 neurons (colored segments).

Figure 2:

A rank-1 inhibitory population forms a stable, attracting boundary. (a) A single inhibitory neuron’s voltage (see equation 2.10) is visualized in input-output space; the spike threshold (V1=T1; solid blue line, with normal vector (F1,E1)) divides the space into sub- and suprathreshold sets (white and blue shaded area). Gray arrows illustrate dynamics due to the leak; the downward-facing blue arrow illustrates the effect of a spike at the boundary, which sets yy+D1 (dotted blue line). (b) The thresholds of multiple (N=10; colored lines) neurons were tuned to delineate a concave (or negative convex) boundary, y=-fcvx(x)=-x2-12 (dotted gray line). (c) The boundary (blue line) divides the input-output space into a convex subthreshold set (with all voltages below threshold) and a suprathreshold set (with at least one voltage above threshold), forming a stable input-output relationship. Output precision is defined as the distance between the neural thresholds and the resets after a spike, which can be different for each neuron (gray shading; see equation 2.18). (d) Latent variable output (top) and spike raster (bottom) of a simulation of the spiking network from panel c (colors as in panel b) for a time-varying input from x=-1 to 1 (light gray, top; time in units of τ); the spiking simulation output (solid dark gray) follows the true boundary (dashed black line), and each input value is coded by a single neuron. (e) A negative convex surface boundary for a two-dimensional input (see equation 2.16), made up of 36 neurons (colored segments).

Close modal

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 Ei0, 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=EiDj0, 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

The concept of a boundary in the input-output space is central to our developments. This perspective can also be generalized to rank-1 networks with multiple inputs. We write x=(x1,x2,...,xM) for an M-dimensional input, in which case the voltage equation becomes
Vi=Fix+EiyTi.
(2.16)
This inequality still describes a linear boundary, but now in a higher-dimensional input-output space, (x,y)RM+1. Multiple neurons again define different boundaries, each of which divides the input-output space into two half-spaces. The intersection of the individual subthreshold half-spaces yields the population’s subthreshold set, which is guaranteed to be convex. The boundary of this convex set can therefore be written as
y=-fcvx(x),
(2.17)
where fcvx(x) denotes a (piecewise-linear) convex function. The negative sign implies that y is actually a concave function of the inputs x. For xR2, we plot an example boundary in Figure 2e.
The above function serves as an idealized description of the network’s input-output relationship. More precisely, the output of the network jumps back and forth between this threshold boundary and the set of y-values that is reached after a spike (see Figure 2d). Since every input x can be uniquely associated with the neuron whose boundary is exposed, we can define a decoder function, D(x)0, that takes the value of the decoder for the neuron that becomes active for input x. Accordingly, the output takes values in the interval
y=[-fcvx(x)+D(x),-fcvx(x)].
(2.18)
We refer to this deviation from the “true” boundary, -fcvx(x), caused by discrete spiking events, as the precision by which the input is mapped onto the output, as determined by the function D(x) (compare Figure 2c, gray shading, and Figure 2d).

2.4  The Boundary Determines the Latent Dynamics

We can write down the dynamics of y in two equivalent ways. First, we can think of each individual neuron as enforcing its own threshold, and we can think of the boundary as the joint action of all thresholds. Merging the two dynamical regimes, equation 2.15, into a single equation, we obtain
y˙=-y+i=1NDiIFix+Eiy-Ti,
(2.19)
where the indicator function I(·) denotes an infinitely high boundary, I(z)= if z0 and I(z)=0 otherwise. The indicator function is commonly used in convex analysis (Rockafellar, 1970; Boyd & Vandenberghe, 2004) and can be understood as similar to the delta function used in spiking networks, that is, as a limiting case, I(z)=limΔt0H(z)/Δt, where H(z) is the Heaviside function and Δt the integration time step.
Second, we can think of the network as forming a single (globally) stable boundary,
y˙=-y+D(x)IB(x,y),
(2.20)
using the decoder function D(x) and similarly defining a boundary function,
B(x,y)=y+fcvx(x),
(2.21)
where we can recover our definition of the boundary input-output function by setting B(x,y)=0. Equations 2.19 and 2.20 are mathematically equivalent, and both will prove useful later on.

2.5  Excitatory Neurons Form Unstable, Repellent Boundaries

Next, we consider a network of N recurrently connected, excitatory neurons. Much of our exposition follows the same outline as for the inhibitory network, and we will mostly focus on highlighting the differences. Of course, the key difference is positive connectivity, compared to equation 2.11,
Wij=EiDj0,
(2.22)
with constraints on the encoders and decoders,
Ei0i,
(2.23)
Dj0j.
(2.24)
Given these sign conventions for the excitatory network, the latent variable or readout will be constrained to be positive, y0. The voltage equation 2.10, remains the same, but now each neuron receives positive feedback from itself and all other neurons of the excitatory population.

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.

Figure 3:

A rank-1 excitatory population forms an unstable, repellent boundary. (a) Input-output space for a single excitatory neuron (compare with and refer to the caption of Figure 2a). Note that a spike moves the latent variable further into the suprathreshold set, and the leak dynamics decay away from the boundary. (b) The thresholds of multiple (N=10, colored lines) excitatory neurons were tuned to delineate a concave (or negative convex) boundary in input-output space, y=-fcvx(x)=-x2+32 (dotted gray line). (c) The boundary forms a convex subthreshold set (where all voltages are below threshold), analogous to the inhibitory population (see Figure 1). Unlike the inhibitory population, however, the boundary itself is now unstable: it either leaks to zero or explodes.

Figure 3:

A rank-1 excitatory population forms an unstable, repellent boundary. (a) Input-output space for a single excitatory neuron (compare with and refer to the caption of Figure 2a). Note that a spike moves the latent variable further into the suprathreshold set, and the leak dynamics decay away from the boundary. (b) The thresholds of multiple (N=10, colored lines) excitatory neurons were tuned to delineate a concave (or negative convex) boundary in input-output space, y=-fcvx(x)=-x2+32 (dotted gray line). (c) The boundary forms a convex subthreshold set (where all voltages are below threshold), analogous to the inhibitory population (see Figure 1). Unlike the inhibitory population, however, the boundary itself is now unstable: it either leaks to zero or explodes.

Close modal
Just as in the inhibitory network, the set of subthreshold readouts is a convex set, and the boundary is a concave function of the input within a given range (see Figure 3). The same reasoning holds when the input x becomes multidimensional. We can therefore characterize the unstable boundary through a concave function,
y=-fcvx(x),
(2.25)
defined on the set of inputs, x, for which y0. Since this boundary is unstable, it does not describe the input-output mapping of an all-excitatory network. Indeed, the input-output mapping of the network will be determined by the mechanism(s) that stabilize the explosion (see below). Finally, the dynamics of the excitatory network also follow equations 2.19 and 2.20, with the key difference that Di or D(x) are now positive, which inverts the dynamics around the boundary.
Given that the excitatory population is unstable, some other mechanisms need to kick in to stabilize the dynamics. One stabilizing factor is the leak of each neuron, which eventually could counter the explosion of activity. A second approach would be to give each excitatory neuron a self reset and/or refractory period after firing a spike. This can be formalized by adding a diagonal component μ to the connectivity matrix such that
Wii=EiDi-μ,
(2.26)
which serves as a “soft” refractory period. Though technically this breaks the assumed low-rank connectivity structure, it will prove useful, and we return to it in section 3. Finally, a third possibility is to stabilize the network activity with inhibition, which we do next.

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

We now study networks of coupled excitatory (E) and inhibitory (I) neurons. In order to do so, we will first go back to differential equations. Assuming that there are NI inhibitory and NE excitatory neurons, we again use leaky integrate-and-fire neurons to describe the membrane potentials of the two populations:
V˙iI=-ViI+FiIc(t)+j=1NEWijIEsjE(t)+j=1NIWijIIsjI(t),
(2.27)
V˙iE=-ViE+FiEc(t)+j=1NEWijEEsjE(t)+j=1NIWijEIsjI(t).
(2.28)
Besides the self-connections within the subnetworks, WEE and WII, we also introduce cross-connections between the two networks, designated by the matrices WEI and WIE. The thresholds of the two populations are given by TiI and TiE and, as before, may incorporate possible background inputs.
Just as before, we assume that the self-connection matrices are rank-1. We furthermore assume that the cross-connection matrices are likewise rank-1, while sharing the same decoders. Specifically, we set (for all i,j)
WijII=EiIIDjI0,
(2.29)
WijIE=EiIEDjE0,
(2.30)
WijEI=EiEIDjI0,
(2.31)
WijEE=EiEEDjE0,
(2.32)
resulting in a rank-2 excitatory-inhibitory (EI) network. These settings are somewhat restrictive, as the cross-connection matrices do not have their own set of decoders. We defer the analysis of the general case to section 3. All of the encoding weights are assumed to be positive, while the decoders follow previously specified sign constraints, DjI0 and DjE0, ensuring Dale’s law.
Given the above assumptions, we define two latent variable readouts,
yI=j=1NIDjIrjI,
(2.33)
yE=j=1NEDjErjE,
(2.34)
which allows us to integrate the differential equations and obtain
ViI=FiIx+EiIEyE+EiIIyITiI,
(2.35)
ViE=FiEx+EiEEyE+EiEIyITiE.
(2.36)
We see that even in this more complicated network, each neuron can once more be interpreted as a bound in input-output space. The key difference is that we now have a rank-2 network, and the input-output space has become three-dimensional, given by (x,yE,yI)R3. As a consequence, neural thresholds have become planes instead of lines. Furthermore, we now have two distinct sets of spike-threshold boundaries, corresponding to the two populations.
The spiking dynamics of the network can again be understood by focusing on the space of latent variables, which is now two-dimensional. Taking the derivative of the readout equations, we obtain the two-dimensional dynamics of the latent space; compare equation 2.19:
y˙I=-yI+i=1NIDiIIFiIx+EiIEyE+EiIIyI-TiI,
(2.37)
y˙E=-yE+i=1NEDiEIFiEx+EiEEyE+EiEIyI-TiE.
(2.38)
We have already done all the work to understand these equations. In the first equation, for instance, we can simply understand yE as a (time-varying) external input. As a consequence, we can treat the inhibitory population as receiving a two-dimensional input, (x,yE; see Figure 2e), so that the derivations from section 2.2 all remain the same. Similarly, the excitatory population can be treated as receiving a two-dimensional input (x,yI). Accordingly, the inhibitory dynamics generates a stable boundary in (x,yE,yI)-space, and the excitatory dynamics generates an unstable boundary in (x,yE,yI)-space.

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.

Figure 4:

The inhibitory boundary can stabilize the excitatory boundary. (a) The inhibitory population, with latent readout yI, forms a negative convex, attracting boundary in (x,yE) space. (b) The excitatory population, with latent readout yE, forms a negative convex, repellent boundary in (x,yI) space. (c) For fixed x (blue and red lines in panels a and b, respectively, for x=4), the inhibitory and excitatory boundaries can be viewed in (yE,|yI|) activity space, and delineate four distinct regions of the state space: E and I subthreshold (gray), E suprathreshold and I subthreshold (red), E subthreshold and I suprathreshold (blue), and both E and I suprathreshold (purple). (d, e) Single-neuron E and I boundaries with spiking simulation dynamics (black trajectory in panel d; blue and red solid lines in panel e (simulation)) around the boundary crossing (black dot in panel d; blue and red dashed lines in panel e). Dots in panel e show the spikes fired. (f) Heterogeneous E and I boundaries, in which I (respectively, E) neurons (faded blue and red lines) form a piecewise-linear approximation to the idealized convex boundaries (opaque blue and red lines). Dynamics are shown in black, oscillating around the boundary crossing (black dot).

Figure 4:

The inhibitory boundary can stabilize the excitatory boundary. (a) The inhibitory population, with latent readout yI, forms a negative convex, attracting boundary in (x,yE) space. (b) The excitatory population, with latent readout yE, forms a negative convex, repellent boundary in (x,yI) space. (c) For fixed x (blue and red lines in panels a and b, respectively, for x=4), the inhibitory and excitatory boundaries can be viewed in (yE,|yI|) activity space, and delineate four distinct regions of the state space: E and I subthreshold (gray), E suprathreshold and I subthreshold (red), E subthreshold and I suprathreshold (blue), and both E and I suprathreshold (purple). (d, e) Single-neuron E and I boundaries with spiking simulation dynamics (black trajectory in panel d; blue and red solid lines in panel e (simulation)) around the boundary crossing (black dot in panel d; blue and red dashed lines in panel e). Dots in panel e show the spikes fired. (f) Heterogeneous E and I boundaries, in which I (respectively, E) neurons (faded blue and red lines) form a piecewise-linear approximation to the idealized convex boundaries (opaque blue and red lines). Dynamics are shown in black, oscillating around the boundary crossing (black dot).

Close modal
Following equation 2.20, we can also rewrite the dynamics in terms of the boundaries only,
y˙I=-yI+DI(x,yE)IBI(x,yE,yI),
(2.39)
y˙E=-yE+DE(x,yI)IBE(x,yE,yI),
(2.40)
where the boundary functions are given by BI(x,yE,yI)=fcvxI(x,yE)+yI and BE(x,yE,yI)=fcvxE(x,yI)+yE, and the boundaries are retrieved by setting them to zero. To visualize the dynamics of this joint system, we assume a constant input of x=4, which means that the two boundary surfaces reduce to 1D curves in (yE,yI)-space (shown as the blue and red lines in Figures 4a and 4b, respectively). Both 1D curves are plotted together in Figure 4c. Here, we plot the absolute value of the inhibitory readout, |yI|, as it is comparable to the summed population activity (especially if decoders are constant, DI(x,yE)=const), and thereby relates to typical ways of looking at EI networks (Wilson & Cowan, 1972; Dayan & Abbott, 2005; Gerstner et al., 2014).

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.

The latent dynamics of the EI system in Figure 4c will tend to gravitate around the point at which the two boundaries cross. The dynamics resemble a stable limit cycle that, in the limit as the decoders become smaller and smaller, effectively becomes a stable fixed point for certain requirements of the two boundaries. Intuitively, we can understand the stability of the EI system by considering the slopes of the two boundaries in (yE,|yI|)-space (see Figure 4d). Under the rule that inhibitory neurons always fire before excitatory neurons, local stability requires the inhibitory boundary to have a steeper slope, such that it is able to push the dynamics into the subthreshold area when the readouts are above the crossing point. For instance, if the network is composed of one inhibitory neuron and one excitatory neuron (see Figure 4d), this slope condition can be succinctly expressed through the encoding weights as
E1IEE1II>E1EEE1EI,
(2.41)
where both neurons have been labeled with index i=1. This inequality is analogous to the condition on the determinant in a two-dimensional linear stability analysis (Izhikevich, 2007). A more formal treatment of stability can be done after relating the spiking dynamics to those of a rate network (see section 4; Dayan & Abbott, 2005; Izhikevich, 2007). We also note that other stable dynamical regimes may be present in the more general case, such as networks in which the excitatory or both populations are silent. As before, we consider these to be degenerate and do not discuss them here.

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.

Figure 5:

Nonlinear function approximation in the rank-2 EI network. (a) Each of three input levels (illustrated for x=0, 5, and 10) results in a different intersection of the excitatory (red lines) and inhibitory (blue lines) population boundaries, thereby modulating the latent outputs of the network. Note that as x changes, the pair of E and I neurons that form the boundary intersection may change (not shown, but see panel c). The red (or, respectively, blue) shading represents the area for which the excitatory (respectively, inhibitory) boundary is above the other with respect to |yI|, highlighting the changing boundary intersection as a function of x. (b) A desired saw-like output function (top) can be decomposed into the difference of two convex functions q(x) and p(x) (middle), such that the excitatory latent variable approximates the function (red curve, bottom). (c) Visualization of a network of excitatory (NE=3 red planes) and inhibitory (NI=3 blue planes) neurons that approximates the saw-like function following panel b. Note that the spiking dynamics (black faded line; see panel d) follows the boundary intersection as x increases. (d) Simulating the spiking network from panel c for an input x from 0 to 10 (top) confirms that the excitatory (red) and inhibitory (blue) readouts follow the boundaries (bottom). The spike trains demonstrate how different pairs of E and I neurons code for the boundary at different values of x (middle).

Figure 5:

Nonlinear function approximation in the rank-2 EI network. (a) Each of three input levels (illustrated for x=0, 5, and 10) results in a different intersection of the excitatory (red lines) and inhibitory (blue lines) population boundaries, thereby modulating the latent outputs of the network. Note that as x changes, the pair of E and I neurons that form the boundary intersection may change (not shown, but see panel c). The red (or, respectively, blue) shading represents the area for which the excitatory (respectively, inhibitory) boundary is above the other with respect to |yI|, highlighting the changing boundary intersection as a function of x. (b) A desired saw-like output function (top) can be decomposed into the difference of two convex functions q(x) and p(x) (middle), such that the excitatory latent variable approximates the function (red curve, bottom). (c) Visualization of a network of excitatory (NE=3 red planes) and inhibitory (NI=3 blue planes) neurons that approximates the saw-like function following panel b. Note that the spiking dynamics (black faded line; see panel d) follows the boundary intersection as x increases. (d) Simulating the spiking network from panel c for an input x from 0 to 10 (top) confirms that the excitatory (red) and inhibitory (blue) readouts follow the boundaries (bottom). The spike trains demonstrate how different pairs of E and I neurons code for the boundary at different values of x (middle).

Close modal
We assume that the inhibitory population boundary has the form
yI=-fcvxI(x,yE)=-p(x)-ayE,
(2.42)
where p(x) is some convex function and a is a positive constant. Then we assume the excitatory population boundary has the form
yE=-fcvxE(x,yI)=-q(x)-yI,
(2.43)
where q(x) is also a convex function. We note that the stability condition, equation 2.41, requires a>1, which ensures that the inhibitory boundary has a larger slope than the excitatory boundary, and we will use a=2 here. Given that the input-output function will be described by the crossing of the two boundaries, we can rewrite the excitatory latent variable as
yE=-q(x)-yI
(2.44)
=-q(x)+p(x)+2yE.
(2.45)
Finally, rearranging terms, we obtain
yE=q(x)-p(x),
(2.46)
which is the difference of two convex functions. Instead solving for yI, we get
yI=p(x)-2q(x),
(2.47)
which is again the difference of two convex functions. Since any continuous function with a bounded second derivative can be expressed as the difference of two convex functions (Yuille & Rangarajan, 2003), the input-output function of the rank-2 E-I network is fully general (see appendix E for more details). Accordingly, the limitations to computations imposed by the boundaries of the separate populations disappear once excitatory and inhibitory populations are coupled. Furthermore the interpretation of the input-output transformation as computing a difference of two convex functions puts forth an intriguing computational hypothesis for the function of Dale’s law (see section 5).

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.

3.1  Noise Makes Boundaries Jitter

To study the addition of noise, we begin by returning to the single, inhibitory population. Concretely, we assume that all neurons are subject to small, independent white noise currents, ηi(t), which we add to the differential equation, equation 2.1, of each neuron. This white noise will be filtered through the leaky integration of the membrane, so that the integrated voltage, equation 2.8, becomes
Vi=Fix+Eiy+h*η(t)Ti,
(3.1)
where h(·) describes a one-sided exponential kernel and * denotes a convolution. Following Calaim et al. (2022), we move the noise term onto the threshold,
Fix+EiyTi-h*η(t),
(3.2)
and reinterpret this equation geometrically. From this view, the white noise causes random shifts in the precise position (but not the orientation) of a neuron’s boundary in input-output space (see Figure 6a). Additionally, and again following Calaim et al. (2022), we now consider neurons to have an additional self reset, following equation 2.26. This acts as a soft refractory period and can be modeled as an additional threshold term,
Fix+EiyTi-h*η(t)+μri(t),
(3.3)
which increases the neuron’s threshold after spiking (recall that ri(t) is the filtered spike train; see equation 2.4). This can be considered as another mechanism that causes temporary shifts in the threshold of each neuron. In this case, the threshold will jump away from the population boundary by a fixed value after spiking and exponentially decay back to its original location. At the network level, due to both of these effects, we can thus see that the boundaries of all neurons will jitter around their default positions, as illustrated for two noise snapshots in Figure 6a, resulting in a noisy population boundary.
Figure 6:

Noise causes jittery boundaries and irregular firing in the inhibitory population. (a) Current noise makes individual neural threshold boundaries (faded colored lines) jitter up and down independently, resulting in a fluctuating population boundary (opaque multicolored line) around the default boundary function (gray dotted line). Two different realizations of noise are shown (top and bottom) and illustrate that each input value may be coded by a different neuron at different times. (b) Input x (light gray) and latent readout y (dark gray) for a spiking simulation of a network with N=50 neurons, with a small amount of independent current noise added. Default boundary is shown in dashed black (compare with panel a, dotted gray). (c) Spike trains for the simulation in panel b; note that several neurons become active for any given input x. (d, e) Voltages of two example neurons (highlighted in panel c for trial 1) for two trials from the network in panels b and c. Note that the spike trains are variable despite the high precision of the output y.

Figure 6:

Noise causes jittery boundaries and irregular firing in the inhibitory population. (a) Current noise makes individual neural threshold boundaries (faded colored lines) jitter up and down independently, resulting in a fluctuating population boundary (opaque multicolored line) around the default boundary function (gray dotted line). Two different realizations of noise are shown (top and bottom) and illustrate that each input value may be coded by a different neuron at different times. (b) Input x (light gray) and latent readout y (dark gray) for a spiking simulation of a network with N=50 neurons, with a small amount of independent current noise added. Default boundary is shown in dashed black (compare with panel a, dotted gray). (c) Spike trains for the simulation in panel b; note that several neurons become active for any given input x. (d, e) Voltages of two example neurons (highlighted in panel c for trial 1) for two trials from the network in panels b and c. Note that the spike trains are variable despite the high precision of the output y.

Close modal

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 irregularity or noisiness of the individual spike trains of the network has important consequences for decoding. To visualize these consequences, we shift from the low-dimensional latent space back to the N-dimensional space spanned by the filtered spike trains or firing rates of the neurons, r=(r1,r2,...,rN). We first observe that each neuron’s threshold also describes a hyperplane in firing rate space. When the voltage reaches the threshold, the terms of the neuron’s voltage equation (equation 3.3) can be rewritten and solved for y as
y=j=1NDjrj=Ti-h*η(t)+μri(t)-FixEi,
(3.4)
which defines a hyperplane with normal vector D=(D1,...,DN) and offset as defined by the right-hand side. The hyperplanes of different neurons have the same orientation, but different offsets. Two of these hyperplanes are shown in Figure 7a for a two-neuron network.
Figure 7:

Irregular firing amplifies noise in the decoder nullspace in the inhibitory population. (a) In the rank-1 inhibitory network, all neural thresholds are parallel in firing rate space (shown for two example neurons with rates r1, dark blue, and r2, green). For a deterministic network, only one of the neurons will fire (blue, r1, dynamics in black), eventually reaching a fixed point (shown in panel b). (b) The latent readout subspace (light blue) is orthogonal to the thresholds. The gray dotted line shows the subspace orthogonal to the readout (latent nullspace). (c) Simulation of the steady-state spiking dynamics of the deterministic, two-neuron network from panels a and b. Both latent readout and activity projected into the orthogonal subspace are noise-free. (d) For a network contaminated by noise, the neural thresholds fluctuate in firing rate space. As a consequence, both the green and blue neuron participate in the dynamics (black). (e) With both neurons taking random turns in firing, the dynamics diffuse along the orthogonal subspace. (f) For the network with noise, the readout is no longer deterministic, but the noise is relatively small. In comparison, activity projected into the orthogonal subspace fluctuates wildly.

Figure 7:

Irregular firing amplifies noise in the decoder nullspace in the inhibitory population. (a) In the rank-1 inhibitory network, all neural thresholds are parallel in firing rate space (shown for two example neurons with rates r1, dark blue, and r2, green). For a deterministic network, only one of the neurons will fire (blue, r1, dynamics in black), eventually reaching a fixed point (shown in panel b). (b) The latent readout subspace (light blue) is orthogonal to the thresholds. The gray dotted line shows the subspace orthogonal to the readout (latent nullspace). (c) Simulation of the steady-state spiking dynamics of the deterministic, two-neuron network from panels a and b. Both latent readout and activity projected into the orthogonal subspace are noise-free. (d) For a network contaminated by noise, the neural thresholds fluctuate in firing rate space. As a consequence, both the green and blue neuron participate in the dynamics (black). (e) With both neurons taking random turns in firing, the dynamics diffuse along the orthogonal subspace. (f) For the network with noise, the readout is no longer deterministic, but the noise is relatively small. In comparison, activity projected into the orthogonal subspace fluctuates wildly.

Close modal

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

Figure 8:

Irregular firing, balance, and noise control in the rank-2 EI network. (a) Population boundaries for a homogeneous population of NE=50 excitatory (E, red) and NI=50 inhibitory (I, blue) neurons with fully aligned cross-connections and stable spiking dynamics at the boundary crossing (black); one additional inactive excitatory neuron with different parameters is shown (E2, pink). (b) When all inhibitory neurons are stimulated, the original boundary (dashed blue) moves upward (solid blue), and the fixed point is shifted down (spiking simulation in black). (c) Latent readouts (top) and spike rasters (bottom) of the network over time; the inhibitory population is stimulated at t=10, which shifts the boundary intersection (dashed lines) and the spiking dynamics (solid faded lines). (d) Balanced inputs (top) and voltage (bottom, red) for an example excitatory neuron from the homogeneous population (red); the additional excitatory neuron (E2, pink) is hyperpolarized. (e) Illustration of how mistuned decoders can amplify noise, shown in (r1E,r2E) firing rate space (compare with Figure 7); spiking dynamics (faded gray) cause fluctuations in the nullspace of DEE (dotted red), which will be correlated with the latent readout DIE (solid pink). (f) Average mean squared error (MSE) between true boundary and spiking simulation as a function of mistuning (angle between decoders; shading indicates standard deviation). (g) Example simulation for misaligned decoders (indicated by the star in panel f), resulting in noisy fluctuations.

Figure 8:

Irregular firing, balance, and noise control in the rank-2 EI network. (a) Population boundaries for a homogeneous population of NE=50 excitatory (E, red) and NI=50 inhibitory (I, blue) neurons with fully aligned cross-connections and stable spiking dynamics at the boundary crossing (black); one additional inactive excitatory neuron with different parameters is shown (E2, pink). (b) When all inhibitory neurons are stimulated, the original boundary (dashed blue) moves upward (solid blue), and the fixed point is shifted down (spiking simulation in black). (c) Latent readouts (top) and spike rasters (bottom) of the network over time; the inhibitory population is stimulated at t=10, which shifts the boundary intersection (dashed lines) and the spiking dynamics (solid faded lines). (d) Balanced inputs (top) and voltage (bottom, red) for an example excitatory neuron from the homogeneous population (red); the additional excitatory neuron (E2, pink) is hyperpolarized. (e) Illustration of how mistuned decoders can amplify noise, shown in (r1E,r2E) firing rate space (compare with Figure 7); spiking dynamics (faded gray) cause fluctuations in the nullspace of DEE (dotted red), which will be correlated with the latent readout DIE (solid pink). (f) Average mean squared error (MSE) between true boundary and spiking simulation as a function of mistuning (angle between decoders; shading indicates standard deviation). (g) Example simulation for misaligned decoders (indicated by the star in panel f), resulting in noisy fluctuations.

Close modal

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

We now finally consider the general case of the EI network where each population utilizes independent decoders to read out the latent variables. That is, in contrast to the previous shared decoder case in equations 2.30 and 2.31 we now define the cross-connections as
WijIE=EiIEDjIE,
(3.5)
WijEI=EiEIDjEI,
(3.6)
where DjIE and DjEI are potentially distinct from the self-decoders DjE and DjI, respectively. With this redefinition of equations 2.30 and 2.31, we are now looking at a rank-4 EI network, in which all cross- and self-connections can be arbitrary (sign-constrained) rank-1 matrices. The cross-connection matrices can now be “misaligned,” and so they will generally read out a mixture between each population’s own latent readouts, yE and yI, and the orthogonal subspace discussed above (see Figure 8e). To formalize these cross-connection readouts, we introduce two new latent variables,
yIE=j=1NEDjIErjE,
(3.7)
yEI=j=1NIDjEIrjI.
(3.8)
These misaligned readouts can be explicitly split into a part that aligns with the self-decoders and a part that captures noise from the null space. For example, let us focus on the excitatory readout of the inhibitory population. Writing DX=(D1X,...,DNX) for the vector of decoders for all four cases, X{E,EI,IE,I}, and assuming for simplicity that the four decoding vectors are normalized, (DX)DX=1, we obtain
yIE=(DIE)rE
(3.9)
=(DIE)DE(DE)+I-DE(DE)rE
(3.10)
=(DIE)DEyE+DIEI-DE(DE)rE
(3.11)
=αyE+(DIE-αDE)rE.
(3.12)
Here, the first term captures a decreased readout of the correct excitatory latent, yE, since α=(DIE)DE1, and the equality sign applies only when DIE=DE. In turn, the second term captures the random fluctuations from the orthogonal subspace. As the decoders DIE and DE become less and less aligned, α decreases, and the relative power of signal and noise shifts toward noise (see Figures 8f and 8g; see also Landau & Sompolinsky, 2021). The respective networks will thereby amplify noise and contaminate the signal. When the two readouts become orthogonal, α=0, the inhibitory population receives only noise and no longer any signal, yE. In this extreme case, the stability will be compromised (not shown).

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

First, we examine how the boundary is affected when we change our assumptions about postsynaptic currents. For simplicity, we first return to a single, rank-1 inhibitory population. Previously, we demonstrated that the dynamics of y are characterized by the combination of a leak and an infinitely steep boundary (compare equation 2.19 but assuming a one-dimensional input),
y˙=-y+i=1NDiIFix+Eiy-Ti.
(4.1)
For a fixed input x, we illustrate this infinitely steep boundary in Figure 9a, where we plot y versus dy/dt (consider this as a vertical slice through Figures 2a to 2c). The dynamics at the boundary are oscillatory, with leak dynamics, followed by instantaneous spiking (see Figure 2d). The latent readout (see equation 2.14) can also be expressed in terms of the respective spike times, tif,
y˙=-y+i=1NfDiδ(t-tif),
(4.2)
where tif is the fth spike fired by the ith neuron. Here, the delta function models the effect of a spike on the latent variable. Biophysically, the latent variable y is linearly related to the voltage (see equation 2.10) and so the delta function corresponds to an infinitely short, postsynaptic current.
Figure 9:

From hard to soft boundaries in the inhibitory population. (a) Plot of latent y versus its derivative dy/dt for a synaptic delta current (analogous to a vertical slice through Figures 2a and 2c). The threshold of the neuron sitting at y=y0=-2 causes an infinitely steep boundary and divides the subthreshold (white) from suprathreshold (blue) areas. (b) Introduction of finite-time, square pulse synaptic dynamics, with time constant τs (shown for three values of τs: in black, gray, and silver) results in oscillatory dynamics in (y, sk)-space, with jumps in sk of size 1/τs. Arrows indicate equally spaced time intervals. (c) The time-averaged dynamics of y, denoted y¯, can be qualitatively approximated by assuming fixed synaptic input of amplitude 1/τs when the boundary is crossed (at y¯=y¯0), resulting in sk(y¯) as a scaled Heaviside function. (d) The resulting plot of the average latent y¯ versus its derivative for the approximation to the square pulse synapse, sk(y¯), from panel c, where the jump in the derivative is a negative Heaviside function scaled by Dk/τs. (e) Averaging the spiking dynamics for a noisy threshold boundary from panel d (gray, simulation) results in a soft boundary well approximated by a sigmoid plus leak (blue). Trials were simulated with random initial conditions for y. (f) The soft boundary equivalent to the hard inhibitory network boundary from Figures 2b to 2d and Figure 6, in which neurons are arranged along a quadratic function.

Figure 9:

From hard to soft boundaries in the inhibitory population. (a) Plot of latent y versus its derivative dy/dt for a synaptic delta current (analogous to a vertical slice through Figures 2a and 2c). The threshold of the neuron sitting at y=y0=-2 causes an infinitely steep boundary and divides the subthreshold (white) from suprathreshold (blue) areas. (b) Introduction of finite-time, square pulse synaptic dynamics, with time constant τs (shown for three values of τs: in black, gray, and silver) results in oscillatory dynamics in (y, sk)-space, with jumps in sk of size 1/τs. Arrows indicate equally spaced time intervals. (c) The time-averaged dynamics of y, denoted y¯, can be qualitatively approximated by assuming fixed synaptic input of amplitude 1/τs when the boundary is crossed (at y¯=y¯0), resulting in sk(y¯) as a scaled Heaviside function. (d) The resulting plot of the average latent y¯ versus its derivative for the approximation to the square pulse synapse, sk(y¯), from panel c, where the jump in the derivative is a negative Heaviside function scaled by Dk/τs. (e) Averaging the spiking dynamics for a noisy threshold boundary from panel d (gray, simulation) results in a soft boundary well approximated by a sigmoid plus leak (blue). Trials were simulated with random initial conditions for y. (f) The soft boundary equivalent to the hard inhibitory network boundary from Figures 2b to 2d and Figure 6, in which neurons are arranged along a quadratic function.

Close modal
In reality, of course, the postsynaptic current generated by a spike will last for a finite amount of time. We can explicitly include these synaptic dynamics by replacing the delta function in equation 4.2 with a model of finite synaptic dynamics, such as exponential decay or a square pulse (Gerstner et al., 2014). Taking the square pulse as an example, we write its synaptic dynamics as
α(t)=1τsH(t)H(τs-t),
(4.3)
where the synaptic time constant τs is given in units of the membrane time constant, since τ=1 in our case. Equation 4.3 scales α(t) with 1/τs, such that it always integrates to one and recovers the Dirac delta in the limit as τs0. Repeating the derivations in sections 2.1 and 2.2, the latent dynamics becomes (compare equation 2.14)
y˙=-y+i=1NfDiα(t-tif)
(4.4)
=-y+i=1NDisi(t),
(4.5)
where in the last step we redefined the “spike train,” si(t), which should now primarily be interpreted as a synaptic current input, consisting of a sum of postsynaptic currents.

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.

To simplify the dynamics, we will aim to qualitatively describe an approximation of y, which we will denote as y¯, by replacing the temporal dynamics in sk(t) with a simple static function s¯k(y¯). We will take advantage of the fact that the square pulse synapse model only takes on one of two discrete values. Qualitatively, we see that sk sits at zero until the boundary is reached, when it jumps up to a nonzero value momentarily, sk1/τs, at which point y moves negatively back into the subthreshold area, until sk jumps back to zero. The timescale τs affects the size of the jump in sk, as well as the size of the change in y. This change in y is bounded by Dk (if τs0), but gets smaller and slower as τs increases. In order to replace this oscillatory dynamics by an effective time-averaged dynamics, we assume a functional form of s¯k as a scaled Heaviside function: it jumps up by 1/τs whenever the boundary is crossed and jumps down to zero when y moves back into the subthreshold set (see Figure 9c; three values of τs shown, corresponding to Figure 9b). As a consequence, we can approximate each neuron’s si as
s¯i(y¯)=1/τsifFix+Eiy¯-Ti>00otherwise.
(4.6)
This approximation allows us to effectively remove the history dependence of the synapses. Back at the level of the latent variable dynamics, we can now write
y¯˙=-y¯+i=1NDis¯i(y¯),=-y¯+i=1NDiτsH(Fix+Eiy¯-Ti).
(4.7)
We have thereby replaced the previously oscillatory dynamics at the boundary with a single fixed point, given when the boundary is reached (strictly speaking, the fixed point of the average dynamics does depend slightly on the timescale τs and will be different from the original y0, but we have here ignored this effect for simplicity).

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.

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

In Figure 9e, we simulated a neuron with normally distributed threshold noise added, and plot the trial-averaged output against the trial-averaged derivative. With a slight abuse of notation, we refer to this trial-averaged latent output as y¯ as well. We see that the main effect of the noise is to soften the boundary. Mathematically, the reason for this softening is that a Heaviside function with input subject to gaussian noise takes the shape of the sigmoidal error function (erf). Due to its common use in rate networks and qualitative similarity to the error function, we choose to use the logistic function,
σβ(u)=11+e-βu,
(4.8)
where the parameter β determines the steepness and is inversely proportional to the noise. For β, we recover the deterministic Heaviside function. Apart from the parameters mentioned before—strength of the decoding weights, the inverse of the synaptic timescale, number of neurons around the boundary—the steepness of the boundary therefore also depends inversely on the level of the noise. When we now put this picture back into the full input-output, (x-y)-space as before, we can visualize the same network boundary as in Figures 2b to 2d and Figure 6, but now for the soft case, which is shown in Figure 9f.
We thus have the following dynamics for the trial-averaged activity of a network with finite synaptic dynamics,
y¯˙=-y¯+i=1NDiτsσβ(Fix+Eiy¯-Ti).
(4.9)
This equation therefore describes the boundary by a single stable fixed point with locally asymmetric dynamics: fast attraction due to a steep boundary on one side and shallow attraction due to the leak on the other side (see Figure 9e). The slower the synapses and the larger the noise, the less strong this asymmetry becomes, until the notion of a boundary ceases to be useful in describing the system dynamics.

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.

Figure 10:

The inhibitory network with slow synapses. (a–c) The network of N=50 neurons from Figure 6 with slow synapses (square pulse; τs=0.5τ) yields stable coding (a) and irregular firing (b,c). (d–f) A rate version of the network also has stable coding (d) and has single-unit rate activities (e) that reproduce the trial-averaged activity of the spiking network (f; gaussian-smoothed spikes for 100 trials of the system in panels a to c).

Figure 10:

The inhibitory network with slow synapses. (a–c) The network of N=50 neurons from Figure 6 with slow synapses (square pulse; τs=0.5τ) yields stable coding (a) and irregular firing (b,c). (d–f) A rate version of the network also has stable coding (d) and has single-unit rate activities (e) that reproduce the trial-averaged activity of the spiking network (f; gaussian-smoothed spikes for 100 trials of the system in panels a to c).

Close modal

4.3  The Latent Variables of Spiking and Rate Networks Are Equivalent

The dynamics of the trial-averaged latents has a simple relation to classical firing rate networks. Let us use the notation rF,i to denote the activation of a unit i of a firing rate network. Assuming a one-dimensional input x, the equation for an all-inhibitory rate network can be written as
r˙F,i=-rF,i+σβFix+jWijrF,j-Ti,
(4.10)
where σβ(·) is the sigmoidal activation function of each neuron, equation 4.8, and the threshold term Ti simply serves as a constant negative offset. Setting the inhibitory connection matrix to be rank-1, Wij=EiDj, with Dj0, and defining the readout as yF=1τsiDirF,i, we can compute the derivative of yF to obtain
y˙F=-yF+i=1NDiτsσβFix+EiyF-Ti.
(4.11)
We note that this equation is identical to equations 4.12 and 4.13 if we identify the readout yF with the trial-averaged readout y¯ from the spiking network. Indeed, when we simulate this firing-rate version of the network, we see that it tracks the boundary well, and the firing rates of the individual neurons closely match the trial-averaged rates of the spiking network (see Figures 10d to 10f)).

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¯|,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.

Figure 11:

Dynamics of the rank-2 EI network with soft boundaries. (a) Stable inhibitory dynamics (solid line), composed of a sigmoidal boundary and a leak (dashed lines). Replotted from Figure 9e, using the absolute value of the inhibitory output, |y¯I|. (b) Unstable excitatory dynamics, again composed of a sigmoidal boundary and a leak. (c, d) The soft inhibitory and excitatory boundaries for the rank-2 EI network in 2D latent space (homogeneous populations). Nullclines are indicated by the black curves. (e) I and E nullclines of the homogeneous rate-based rank-2 EI system for three different values of β, controlling the slope of the sigmoids. The boundaries of the spike-based system are shown as dotted lines for comparison. (f) I and E nullclines for a heterogeneous rank-2 EI system, in which neurons are arranged along quadratic functions as in Figures 4c and 4f.

Figure 11:

Dynamics of the rank-2 EI network with soft boundaries. (a) Stable inhibitory dynamics (solid line), composed of a sigmoidal boundary and a leak (dashed lines). Replotted from Figure 9e, using the absolute value of the inhibitory output, |y¯I|. (b) Unstable excitatory dynamics, again composed of a sigmoidal boundary and a leak. (c, d) The soft inhibitory and excitatory boundaries for the rank-2 EI network in 2D latent space (homogeneous populations). Nullclines are indicated by the black curves. (e) I and E nullclines of the homogeneous rate-based rank-2 EI system for three different values of β, controlling the slope of the sigmoids. The boundaries of the spike-based system are shown as dotted lines for comparison. (f) I and E nullclines for a heterogeneous rank-2 EI system, in which neurons are arranged along quadratic functions as in Figures 4c and 4f.

Close modal
When combining the two populations, we once more assume that they share the same decoders, following equations 2.29 to 2.32. Starting from equation 4.11, we then obtain the following equations for the trial-averaged dynamics:
y¯˙I=-y¯I+i=1NIDiIτsIσβFiIx+EiIEy¯E+EiIIy¯I-TiI,
(4.12)
y¯˙E=-y¯E+i=1NEDiEτsEσβFiEx+EiEEy¯E+EiEIy¯I-TiE.
(4.13)
Compared to the original equation for the spiking network, equations 2.37 and 2.38, we again see that we have simply replaced the indicator function I(·) by a sigmoidal function, σβ(·), divided by the timescales of the synaptic current pulses, τsI and τsE, respectively.

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

Figure 12:

Function approximation in the rank-2 EI network with slow (spike-based) synapses and rate dynamics (compare with Figure 5). (a) The example from Figure 5 with slow synapses (τsE=τ/4, τsI=τ/10) yields a reliable yet noisier approximation to the boundary crossing. (b) A rate-based version of the network with stable, reliable dynamics, and with rates that match the trial-averaged activity of the spiking network (not shown). (c) Slowing down the inhibitory synapses sufficiently (τsE=τ/4, τsI=τ/4) leads to more oscillatory dynamics with less reliable coding.

Figure 12:

Function approximation in the rank-2 EI network with slow (spike-based) synapses and rate dynamics (compare with Figure 5). (a) The example from Figure 5 with slow synapses (τsE=τ/4, τsI=τ/10) yields a reliable yet noisier approximation to the boundary crossing. (b) A rate-based version of the network with stable, reliable dynamics, and with rates that match the trial-averaged activity of the spiking network (not shown). (c) Slowing down the inhibitory synapses sufficiently (τsE=τ/4, τsI=τ/4) leads to more oscillatory dynamics with less reliable coding.

Close modal

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

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, τ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.

Variable(s)Description (X,Y{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
τ Membrane time constant (and simulation timescale) 
τs Synaptic time constant 
μ 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. YX
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 ji (single-pop. or for YX
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)= if z0, I(z)=0 otherwise) 
σβ() Sigmoid function with slope parameter β 
Variable(s)Description (X,Y{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
τ Membrane time constant (and simulation timescale) 
τs Synaptic time constant 
μ 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. YX
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 ji (single-pop. or for YX
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)= if z0, I(z)=0 otherwise) 
σβ() Sigmoid function with slope parameter β 

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

Figure 13:

A rank-1 inhibitory population with randomly distributed parameters. (a, b) The thresholds of multiple (N=8, colored lines) neurons are visualized in input-output space, delineating stable, negative-convex boundaries. (a) For each neuron i, parameter values were drawn as independent random variables from gaussian distributions of the form FiN(0,1), EiN(1,14), and TiN(14,116). Two properties are apparent from this “random” network: (1) silent neurons: the population boundary is composed of only three neurons, with the other five remaining silent, and (2) silent activity regimes: for x=0 to x=0.25, the population boundary sits at y=0, meaning that the network will be fully silent. (b) Same network from panel a, but now with each threshold parameter Ti tuned such that each neuron’s boundary is tangent to the concave boundary (gray dotted line) as in Figure 2. Tuning the thresholds largely solves the problems from panel a with silent neurons and activity regimes. However, neurons are not equally distributed along the boundary in this case, so each input level will have varying amounts of coding redundancy.

Figure 13:

A rank-1 inhibitory population with randomly distributed parameters. (a, b) The thresholds of multiple (N=8, colored lines) neurons are visualized in input-output space, delineating stable, negative-convex boundaries. (a) For each neuron i, parameter values were drawn as independent random variables from gaussian distributions of the form FiN(0,1), EiN(1,14), and TiN(14,116). Two properties are apparent from this “random” network: (1) silent neurons: the population boundary is composed of only three neurons, with the other five remaining silent, and (2) silent activity regimes: for x=0 to x=0.25, the population boundary sits at y=0, meaning that the network will be fully silent. (b) Same network from panel a, but now with each threshold parameter Ti tuned such that each neuron’s boundary is tangent to the concave boundary (gray dotted line) as in Figure 2. Tuning the thresholds largely solves the problems from panel a with silent neurons and activity regimes. However, neurons are not equally distributed along the boundary in this case, so each input level will have varying amounts of coding redundancy.

Close modal

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.

For the single rank-1 population with 1D input and 1D output, neural parameters were chosen such that each neuron was tangent to a particular boundary curve, denoted y=-fcvx(x). Given a network of N neurons, boundaries were distributed uniformly along the curve in a particular input interval [xA,xB], resulting in a tangent point (xi,yi) associated with each neuron. Then, each neuron’s voltage equation (see equation 2.10) at Vi=Ti was used to compute the parameters Fi, Ei, and Ti. Given the redundancy of the three parameters, we arbitrarily set the encoding weight, Ei, to be 1 for each neuron. We thus have
Ei=1,
(C.1)
Fi=-ddxfcvx(xi),
(C.2)
Ti=Fixi+Eiyi.
(C.3)
Then, the decoding weight Di was chosen to achieve a particular jump size in the boundary at each value, and with sign constraint according to the population identity (Di<0 for inhibition and Di>0 for excitation).
To give a concrete example, let’s consider the boundary in Figures 2b and 2c, with boundary function y=-x2-12. One neuron, say, neuron i=3, is tangent to the curve at the point (xi=-0.5,yi=-0.75). This neuron’s parameters were set to
Ei=1,
(C.4)
Fi=-ddxfcvx(xi)=2xi=-1,
(C.5)
Ti=Fixi+Eiyi=-0.25,
(C.6)
Di=-0.35.
(C.7)
Note the negative threshold, with the interpretation that the neuron is spontaneously active due to an additional background current (see equation 2.1). To constrain voltages to all be on the same scale, in practice we use a default threshold Ti0=1 for all neurons, and include a bias current bi=Ti0-Ti. For the case of higher, M-dimensional input, as in Figure 2e and equation 2.16, the procedure is analogous: a tangent point is designated for each neuron in M+1-dimensional space, and the feedforward weights simply become an M-dimensional vector with the gradient of the function -fcvx(xi).
Following a similar methodology as above for the single rank-1 population, here we consider fitting neural parameters of each population to a boundary in (x,yE,yI)-space. Specifically, we denote the inhibitory boundary as yI=-fcvxI(x,yE) and the excitatory boundary as yE=-fcvxE(x,yI). Taking the inhibitory boundary as an example, neurons’ tangent points were distributed in a 2D grid in intervals [xA,xB] and [yAE,yBE], such that each neuron had an associated point (xi,yiE,yiI). Then, as before, encoders for self-connections, EII, were set arbitrarily to 1 for all neurons. All parameters thus were set as
EiII=1,
(D.1)
FiI=-xfcvx(xi,yiE),
(D.2)
EiIE=-yEfcvx(xi,yiE),
(D.3)
TiI=FiIxi+EIEyiE+EiIIyiI.
(D.4)
An analogous procedure was followed for the excitatory population. Once again, in principle decoders could be set arbitrarily following sign conventions (DiI<0, DiE>0). However, we found that in practice, decoders could be optimized to ensure that the two populations take turns in spiking. We generally followed a heuristic that the effective direction that the decoders point toward should aim to counteract the leak dynamics at the boundary crossing point—in other words, DI/DEy*I/y*E, where (y*E,y*I) is the crossing point of the two boundaries.

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:

Theorem 1

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

Proof.

We begin by selecting an arbitrary convex function g(x), whose second derivative is bounded below by a positive constant ε>0. Then there exists a positive constant λ such that the second derivative of q(x)=f(x)+λg(x) is nonnegative. Consequently, q(x) is convex. We then define a second convex function, p(x)=λg(x), and confirm that f(x)=q(x)-p(x)=f(x)+λg(x)-λ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)+λx2 (see Figure 14b, red) and p(x)=λ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.

Figure 14:

Difference of convex function decomposition. (a–d) A smooth, C2 function f(x) (a) is decomposed into the difference of two convex functions q(x) and p(x) (b); first- (c) and second-order (d) gradients show that q'(x) and p'(x) are monotone and q''(x) and p''(x) are nonnegative, confirming them as convex. (e–h) A continuous, piecewise-linear function fPL(x) (e, black curve) is given as an approximation to the smooth curve (panel a; gray curve in e), and is decomposed into the difference of two piecewise-linear convex functions qPL(x) and pPL(x) (f); first- (g) and second-order (h) gradients show that qPL'(x) and pPL'(x) are piecewise-constant and monotone, and qPL''(x) and pPL''(x) are nonnegative, and consist of a series of scaled delta functions.

Figure 14:

Difference of convex function decomposition. (a–d) A smooth, C2 function f(x) (a) is decomposed into the difference of two convex functions q(x) and p(x) (b); first- (c) and second-order (d) gradients show that q'(x) and p'(x) are monotone and q''(x) and p''(x) are nonnegative, confirming them as convex. (e–h) A continuous, piecewise-linear function fPL(x) (e, black curve) is given as an approximation to the smooth curve (panel a; gray curve in e), and is decomposed into the difference of two piecewise-linear convex functions qPL(x) and pPL(x) (f); first- (g) and second-order (h) gradients show that qPL'(x) and pPL'(x) are piecewise-constant and monotone, and qPL''(x) and pPL''(x) are nonnegative, and consist of a series of scaled delta functions.

Close modal

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.

Theorem 2.

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.

Proof.
We consider f(x) in a particular interval of x[a,b]. Since f(x) is piecewise linear, its subdifferential is piecewise constant, with a set of n “edge points” x1 to xn (including the boundary points), which denote the discontinuities and a set of “jumps” v1 to vn-1, where vi indicates how much f'(x) jumps at xi, that is, vi=(f(xi+1)-f(xi))/(xi+1-xi). We note that these sets of points and jumps will fully describe the function f(x) such that we can write it as
f(x)=i=1nH(x-xi)vix,
(E.1)
where H() denotes the Heaviside function. We now wish to decompose f(x) into two piecewise-linear convex functions q(x) and p(x) using the same set of edge points. This is straightforward to do if we group the edge points based on the sign of each jump, vi. We can then write the two functions as
q(x)=i=1nH(x-xi)H(vi)vix,
(E.2)
p(x)=i=1nH(x-xi)H(-vi)|vi|x.
(E.3)
It is then straightforward to show that f(x)=q(x)-p(x).

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.

The above approach is in fact exactly what was used for the results in Figure 5. Once the boundary functions are determined as a function of x, the last step is to define the full boundary surfaces in (x,yE,yI)-space. For simplicity, and as already mentioned in the main text, we assume that the I boundary is linear in yE and that the E boundary is linear in yI, such that the two boundaries can be decomposed into the functional forms:
yI=-p(x)-2yE,
(E.4)
yE=-q(x)-yI.
(E.5)
Then, solving for yE, we confirm that
yE=q(x)-p(x).
(E.6)

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

Figure 15:

Noise and irregular firing in the inhibitory population. (a) The thresholds of multiple neurons (colored lines, shown for N=20) are visualized in input-output space, delineating stable, negative-convex boundaries. The spike trains of one neuron (blue line) were analyzed for a fixed input x=-0.5 (gray vertical line). (b) Spike count coefficient of variation (CV) plotted as a function of the voltage noise standard deviation for the default inhibitory network with recurrent inhibition (black) and comparison with a feedforward variant of the network (identical except that recurrent connections were set to zero). Spike count CV was computed over 500 trials of a constant stimulus x=-0.5 for 5τ time steps.

Figure 15:

Noise and irregular firing in the inhibitory population. (a) The thresholds of multiple neurons (colored lines, shown for N=20) are visualized in input-output space, delineating stable, negative-convex boundaries. The spike trains of one neuron (blue line) were analyzed for a fixed input x=-0.5 (gray vertical line). (b) Spike count coefficient of variation (CV) plotted as a function of the voltage noise standard deviation for the default inhibitory network with recurrent inhibition (black) and comparison with a feedforward variant of the network (identical except that recurrent connections were set to zero). Spike count CV was computed over 500 trials of a constant stimulus x=-0.5 for 5τ time steps.

Close modal

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

In section 4, we eliminate these ad hoc rules and instead introduce finite synaptic dynamics, denoted with the variable si(t), in the form of exponential decay or square pulses. With this modification, the voltage equation becomes
V˙i(t)=-Vi(t)+Fic(t)+j=1NWijsj(t),
(G.1)
where sj(t) is the series of postsynaptic currents due to neuron j.

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

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.

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

Abbott
,
L. F.
,
DePasquale
,
B.
, &
Memmesheimer
,
R.-M
. (
2016
).
Building functional networks of spiking model neurons
.
Nature Neuroscience
,
19
(
3
),
350
355
.
Ahmadian
,
Y.
, &
Miller
,
K. D.
(
2021
).
What is the dynamical regime of cerebral cortex?
Neuron
,
109
(
21
),
3373
3391
.
Ahmadian
,
Y.
,
Rubin
,
D. B.
, &
Miller
,
K. D
. (
2013
).
Analysis of the stabilized supralinear network
.
Neural Computation
,
25
(
8
),
1994
2037
.
Alemi
,
A.
,
Machens
,
C.
,
Deneve
,
S.
, &
Slotine
,
J.-J.
(
2018
).
Learning nonlinear dynamics in efficient, balanced spiking networks using local plasticity rules
. In
Proceedings of the AAAI Conference on Artificial Intelligence
,
32
.
Amit
,
D. J.
, &
Brunel
,
N
. (
1997
).
Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex
.
Cerebral Cortex
,
7
(
3
),
237
252
.
Amos
,
B.
, &
Kolter
,
J. Z
. (
2017
).
OptNet: Differentiable optimization as a layer in neural networks
. In
Proceedings of the International Conference on Machine Learning
(pp.
136
145
).
Bačák
,
M.
, &
Borwein
,
J. M
. (
2011
).
On difference convexity of locally Lipschitz functions
.
Optimization
,
60
(
8–9
),
961
978
.
Baker
,
C.
,
Zhu
,
V.
, &
Rosenbaum
,
R
. (
2020
).
Nonlinear stimulus representations in neural circuits with approximate excitatory-inhibitory balance
.
PLOS Computational Biology
,
16
(
9
), e1008192.
Barak
,
O
. (
2017
).
Recurrent neural networks as versatile tools of neuroscience research
.
Current Opinion in Neurobiology
,
46
,
1
6
.
Barrett
,
D. G.
,
Denève
,
S.
, &
Machens
,
C. K.
(
2013
).
Firing rate predictions in optimal balanced networks
. In
C. J. C.
Burges
,
L.
Bottou
,
M.
Welling
,
Z.
Ghahramani
, &
K. Q.
Weinberger
(Eds.)
,
Advances in neural information processing systems
,
26
.
Curran
.
Barrett
,
D. G.
,
Deneve
,
S.
, &
Machens
,
C. K
. (
2016
).
Optimal compensation for neuron loss
.
eLife
,
5
, e12454.
Bellec
,
G.
,
Salaj
,
D.
,
Subramoney
,
A.
,
Legenstein
,
R.
, &
Maass
,
W.
(
2018
).
Long short-term memory and learning-to-learn in networks of spiking neurons
. In
S.
Bengio
,
H.
Wallach
,
H.
Larochelle
,
K.
Grauman
,
N.
Cesa-Bianchi
, &
R.
Garnett
(Eds.)
,
Advances in neural information processing systems
,
31
.
Curran
.
Bernáez Timón
,
L.
,
Ekelmans
,
P.
,
Kraynyukova
,
N.
,
Rose
,
T.
,
Busse
,
L.
, &
Tchumatchenko
,
T.
(
2022
).
How to incorporate biological insights into network models and why it matters
.
Journal of Physiology
,
601
(
15
),
3037
3053
.
Boerlin
,
M.
,
Machens
,
C. K.
, &
Denève
,
S
. (
2013
).
Predictive coding of dynamical variables in balanced spiking networks
.
PLOS Computational Biology
,
9
(
11
), e1003258.
Boyd
,
S. P.
, &
Vandenberghe
,
L.
(
2004
).
Convex optimization
.
Cambridge University Press
.
Brendel
,
W.
,
Bourdoukan
,
R.
,
Vertechi
,
P.
,
Machens
,
C. K.
, &
Denève
,
S
. (
2020
).
Learning to represent signals spike by spike
.
PLOS Computational Biology
,
16
(
3
), e1007692.
Brunel
,
N
. (
2000
).
Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons
.
Journal of Computational Neuroscience
,
8
,
183
208
.
Calaim
,
N.
,
Dehmelt
,
F. A.
,
Gonçalves
,
P. J.
, &
Machens
,
C. K
. (
2022
).
The geometry of robustness in spiking neural networks
.
eLife
,
11
, e73276.
Chalk
,
M.
,
Gutkin
,
B.
, &
Deneve
,
S
. (
2016
).
Neural oscillations as a signature of efficient coding in the presence of synaptic delays
.
eLife
,
5
, e13824.
Chung
,
S.
, &
Abbott
,
L
. (
2021
).
Neural population geometry: An approach for understanding biological and artificial neural networks
.
Current Opinion in Neurobiology
,
70
,
137
144
.
Cimeša
,
L.
,
Ciric
,
L.
, &
Ostojic
,
S
. (
2023
).
Geometry of population activity in spiking networks with low-rank structure
.
PLOS Computational Biology
,
19
(
8
), e1011315.
Cornford
,
J.
,
Kalajdzievski
,
D.
,
Leite
,
M.
,
Lamarquette
,
A.
,
Kullmann
,
D. M.
, &
Richards
,
B.
(
2021
). Learning to live with Dale’s principle: ANNs with separate excitatory and inhibitory units. In
Proceedings of the 9th International Conference on Learning Representations
(ICLR 2021).
Cunningham
,
J. P.
, &
Yu
,
B. M
. (
2014
).
Dimensionality reduction for large-scale neural recordings
.
Nature Neuroscience
,
17
(
11
),
1500
1509
.
Dayan
,
P.
, &
Abbott
,
L. F.
(
2005
).
Theoretical neuroscience: Computational and mathematical modeling of neural systems
.
MIT Press
.
Denève
,
S.
, &
Machens
,
C. K
. (
2016
).
Efficient codes and balanced networks
.
Nature Neuroscience
,
19
(
3
),
375
382
.
DePasquale
,
B.
,
Sussillo
,
D.
,
Abbott
,
L.
, &
Churchland
,
M. M.
(
2023
).
The centrality of population-level factors to network computation is demonstrated by a versatile approach for training spiking networks
.
Neuron
,
3
(
5
), 631.
Dubreuil
,
A.
,
Valente
,
A.
,
Beiran
,
M.
,
Mastrogiuseppe
,
F.
, &
Ostojic
,
S
. (
2022
).
The role of population structure in computations through neural dynamics
.
Nature Neuroscience
,
25
(
6
),
783
794
.
Eccles
,
J. C
. (
1976
).
From electrical to chemical transmission in the central nervous system: The closing address of the Sir Henry Dale Centennial Symposium Cambridge, 19 September 1975
.
Notes and Records of the Royal Society of London
,
30
(
2
),
219
230
.
Eliasmith
,
C
. (
2005
).
A unified approach to building and controlling spiking attractor networks
.
Neural Computation
,
17
(
6
),
1276
1314
.
Eliasmith
,
C.
, &
Anderson
,
C. H.
(
2003
).
Neural engineering: Computation, representation, and dynamics in neurobiological systems
.
MIT Press
.
Eliasmith
,
C.
, &
Trujillo
,
O
. (
2014
).
The use and abuse of large-scale brain models
.
Current Opinion in Neurobiology
,
25
,
1
6
.
Fourcaud
,
N.
, &
Brunel
,
N
. (
2002
).
Dynamics of the firing probability of noisy integrate-and-fire neurons
.
Neural Computation
,
14
(
9
),
2057
2110
.
Freund
,
T. F.
, &
Katona
,
I
. (
2007
).
Perisomatic inhibition
.
Neuron
,
56
(
1
),
33
42
.
Fusi
,
S.
,
Miller
,
E. K.
, &
Rigotti
,
M
. (
2016
).
Why neurons mix: High dimensionality for higher cognition
.
Current Opinion in Neurobiology
,
37
,
66
74
.
Gallego
,
J. A.
,
Perich
,
M. G.
,
Miller
,
L. E.
, &
Solla
,
S. A
. (
2017
).
Neural manifolds for the control of movement
.
Neuron
,
94
(
5
),
978
984
.
Gerstner
,
W.
,
Kistler
,
W. M.
,
Naud
,
R.
, &
Paninski
,
L.
(
2014
).
Neuronal dynamics: From single neurons to networks and models of cognition
.
Cambridge University Press
.
Gould
,
S.
,
Hartley
,
R.
, &
Campbell
,
D
. (
2021
).
Deep declarative networks
.
IEEE Transactions on Pattern Analysis and Machine Intelligence
,
44
(
8
),
3988
4004
.
Haber
,
A.
, &
Schneidman
,
E.
(
2022
). The computational and learning benefits of Daleian neural networks. In
36th Conference on Neural Information Processing Systems
(NeurIPS 2022).
Haider
,
B.
,
Duque
,
A.
,
Hasenstaub
,
A. R.
, &
McCormick
,
D. A
. (
2006
).
Neocortical network activity in vivo is generated through a dynamic balance of excitation and inhibition
.
Journal of Neuroscience
,
26
(
17
),
4535
4545
.
Hansel
,
D.
, &
Mato
,
G
. (
2013
).
Short-term plasticity explains irregular persistent activity in working memory tasks
.
Journal of Neuroscience
,
33
(
1
),
133
149
.
Hansel
,
D.
, &
van Vreeswijk
,
C
. (
2002
).
How noise contributes to contrast invariance of orientation tuning in cat visual cortex
.
Journal of Neuroscience
,
22
(
12
),
5118
5128
.
Hartman
,
P
. (
1959
).
On functions representable as a difference of convex functions
.
Pacific Journal of Mathematics
,
9
(
3
),
707
713
.
Hennequin
,
G.
,
Ahmadian
,
Y.
,
Rubin
,
D. B.
,
Lengyel
,
M.
, &
Miller
,
K. D
. (
2018
).
The dynamical regime of sensory cortex: Stable dynamics around a single stimulus-tuned attractor account for patterns of noise variability
.
Neuron
,
98
(
4
),
846
860
.
Hennequin
,
G.
,
Vogels
,
T. P.
, &
Gerstner
,
W
. (
2014
).
Optimal control of transient dynamics in balanced networks supports generation of complex movements
.
Neuron
,
82
(
6
),
1394
1406
.
Horst
,
R.
, &
Thoai
,
N. V
. (
1999
).
DC programming: Overview
.
Journal of Optimization Theory and Applications
,
103
,
1
43
.
Huang
,
C.
,
Ruff
,
D. A.
,
Pyle
,
R.
,
Rosenbaum
,
R.
,
Cohen
,
M. R.
, &
Doiron
,
B
. (
2019
).
Circuit models of low-dimensional shared variability in cortical networks
.
Neuron
,
101
(
2
),
337
348
.
Hunsberger
,
E.
, &
Eliasmith
,
C.
(
2015
).
Spiking deep networks with LIF neurons.
.
Ingrosso
,
A.
, &
Abbott
,
L
. (
2019
).
Training dynamically balanced excitatory-inhibitory networks
.
PLOS One
,
14
(
8
), e0220547.
Izhikevich
,
E. M.
(
2007
).
Dynamical systems in neuroscience
.
MIT Press
.
Jazayeri
,
M.
, &
Ostojic
,
S
. (
2021
).
Interpreting neural computations by examining intrinsic and embedding dimensionality of neural activity
.
Current Opinion in Neurobiology
,
70
,
113
120
.
Kadmon
,
J.
,
Timcheck
,
J.
, &
Ganguli
,
S.
(
2020
).
Predictive coding in balanced neural networks with noise, chaos and delays
. In
H.
Larochelle
,
M.
Ranzato
,
R.
Hadsell
,
M. F.
Balcan
, &
H.
Lin
(Eds.)
,
Advances in neural information processing systems
,
33
(pp.
16677
16688
),
Curran
.
Kao
,
T.-C.
, &
Hennequin
,
G
. (
2019
).
Neuroscience out of control: Control-theoretic perspectives on neural circuit dynamics
.
Current Opinion in Neurobiology
,
58
,
122
129
.
Keemink
,
S. W.
, &
Machens
,
C. K
. (
2019
).
Decoding and encoding (de)mixed population responses
.
Current Opinion in Neurobiology
,
58
,
112
121
.
Kim
,
C. M.
, &
Chow
,
C. C
. (
2021
).
Training spiking neural networks in the strong coupling regime
.
Neural Computation
,
33
(
5
),
1199
1233
.
Koren
,
V.
, &
Panzeri
,
S
. (
2022
).
Biologically plausible solutions for spiking networks with efficient coding
. In
S.
Koyejo
,
S.
Mohamed
,
A.
Agarwal
,
D.
Belgrave
,
K.
Cho
, &
A.
Oh
(Eds.)
,
Advances in neural information processing systems
,
35
(pp.
20607
20620
).
Curran
.
Kripfganz
,
A.
, &
Schulze
,
R
. (
1987
).
Piecewise affine functions as a difference of two convex functions
.
Optimization
,
18
(
1
),
23
29
.
Lajoie
,
G.
,
Lin
,
K. K.
,
Thivierge
,
J.-P.
, &
Shea-Brown
,
E
. (
2016
).
Encoding in balanced networks: Revisiting spike patterns and chaos in stimulus-driven systems
.
PLOS Computational Biology
,
12
(
12
), e1005258.
Landau
,
I. D.
, &
Sompolinsky
,
H
. (
2018
).
Coherent chaos in a recurrent neural network with structured connectivity
.
PLOS Computational Biology
,
14
(
12
), e1006309.
Landau
,
I. D.
, &
Sompolinsky
,
H
. (
2021
).
Macroscopic fluctuations emerge in balanced networks with incomplete recurrent alignment
.
Physical Review Research
,
3
(
2
), 023171.
Langdon
,
C.
,
Genkin
,
M.
, &
Engel
,
T. A
. (
2023
).
A unifying perspective on neural manifolds and circuits for cognition
.
Nature Reviews Neuroscience
,
24
,
363
377
.
Li
,
P.
,
Cornford
,
J.
,
Ghosh
,
A.
, &
Richards
,
B.
(
2023
). Learning better with Dale’s Law: A spectral perspective. In
37th Conference on Neural Information Processing Systems
(NeurIPS 2023).
Li
,
Q.
, &
Pehlevan
,
C
. (
2020
).
Minimax dynamics of optimally balanced spiking networks of excitatory and inhibitory neurons
. In
H.
Larochelle
,
M.
Ranzato
,
R.
Hadsell
,
M. F.
Balcan
, &
H.
Lin
(Eds.)
,
Advances in neural information processing systems
,
33
(pp.
4894
4904
).
Curran
.
Lindsay
,
G. W
. (
2021
).
Convolutional neural networks as a model of the visual system: Past, present, and future
.
Journal of Cognitive Neuroscience
,
33
(
10
),
2017
2031
.
Lipp
,
T.
, &
Boyd
,
S
. (
2016
).
Variations and extension of the convex–concave procedure
.
Optimization and Engineering
,
17
,
263
287
.
Mancoo
,
A.
,
Keemink
,
S.
, &
Machens
,
C. K
. (
2020
).
Understanding spiking networks through convex optimization
. In
H.
Larochelle
,
M.
Ranzato
,
R.
Hadsell
,
M. F.
Balcan
, &
H.
Lin
(Eds.)
,
Advances in neural information processing systems
,
33
(pp.
8824
8835
).
Curran
.
Martin
,
K. A.
(
1994
).
A brief history of the “feature detector.”
Cerebral Cortex
,
4
(
1
),
1
7
.
Masset
,
P.
,
Zavatone-Veth
,
J.
,
Connor
,
J. P.
,
Murthy
,
V.
, &
Pehlevan
,
C
. (
2022
).
Natural gradient enables fast sampling in spiking neural networks
. In
S.
Koyejo
,
S.
Mohamed
,
A.
Agarwal
,
D.
Belgrave
,
K.
Cho
, &
A.
Oh
(Eds.)
,
Advances in neural information processing systems
,
35
(pp.
22018
22034
).
Curran
.
Mastrogiuseppe
,
F.
, &
Ostojic
,
S
. (
2018
).
Linking connectivity, dynamics, and computations in low-rank recurrent neural networks
.
Neuron
,
99
(
3
),
609
623
.
Melzer
,
D
. (
1986
).
On the expressibility of piecewise-linear continuous functions as the difference of two piecewise-linear convex functions
. In
V. F.
Demyanov
&
L. C. D.
Dixon
(Eds.)
,
Quasidifferential calculus
(pp.
118
134
).
Springer
.
Miconi
,
T
. (
2017
).
Biologically plausible learning in recurrent neural networks reproduces neural dynamics observed during cognitive tasks
.
eLife
,
6
, e20899.
Mikulasch
,
F. A.
,
Rudelt
,
L.
, &
Priesemann
,
V
. (
2021
).
Local dendritic balance enables learning of efficient representations in networks of spiking neurons
.
Proceedings of the National Academy of Sciences
,
118
(
50
), e2021925118.
Miller
,
K. D.
, &
Troyer
,
T. W
. (
2002
).
Neural noise can explain expansive, power-law nonlinearities in neural response functions
.
Journal of Neurophysiology
,
87
(
2
),
653
659
.
Nardin
,
M.
,
Phillips
,
J. W.
,
Podlaski
,
W. F.
, &
Keemink
,
S. W
. (
2021
).
Nonlinear computations in spiking neural networks through multiplicative synapses
.
Peer Community Journal
,
1
.
Neftci
,
E. O.
,
Mostafa
,
H.
, &
Zenke
,
F
. (
2019
).
Surrogate gradient learning in spiking neural networks: Bringing the power of gradient-based optimization to spiking neural networks
.
IEEE Signal Processing Magazine
,
36
(
6
),
51
63
.
Pang
,
R.
,
Lansdell
,
B. J.
, &
Fairhall
,
A. L
. (
2016
).
Dimensionality reduction in neuroscience
.
Current Biology
,
26
(
14
),
R656
R660
.
Parisien
,
C.
,
Anderson
,
C. H.
, &
Eliasmith
,
C
. (
2008
).
Solving the problem of negative synaptic weights in cortical models
.
Neural Computation
,
20
(
6
),
1473
1494
.
Renart
,
A.
,
De La Rocha
,
J.
,
Bartho
,
P.
,
Hollender
,
L.
,
Parga
,
N.
,
Reyes
,
A.
, &
Harris
,
K. D
. (
2010
).
The asynchronous state in cortical circuits
.
Science
,
327
(
5965
),
587
590
.
Renart
,
A.
,
Moreno-Bote
,
R.
,
Wang
,
X.-J.
, &
Parga
,
N
. (
2007
).
Mean-driven and fluctuation-driven persistent activity in recurrent networks
.
Neural Computation
,
19
(
1
),
1
46
.
Rockafellar
,
R. T.
(
1970
).
Convex analysis
.
Princeton University Press
.
Rosenbaum
,
R.
,
Smith
,
M. A.
,
Kohn
,
A.
,
Rubin
,
J. E.
, &
Doiron
,
B
. (
2017
).
The spatial structure of correlated neuronal variability
.
Nature Neuroscience
,
20
(
1
),
107
114
.
Roudi
,
Y.
, &
Latham
,
P. E
. (
2007
).
A balanced memory network
.
PLOS Computational Biology
,
3
(
9
), e141.
Rubin
,
D. B.
,
Van Hooser
,
S. D.
, &
Miller
,
K. D
. (
2015
).
The stabilized supralinear network: A unifying circuit motif underlying multi-input integration in sensory cortex
.
Neuron
,
85
(
2
),
402
417
.
Rullán Buxó
,
C. E.
, &
Pillow
,
J. W
. (
2020
).
Poisson balanced spiking networks
.
PLOS Computational Biology
,
16
(
11
), e1008261.
Sadeh
,
S.
, &
Clopath
,
C
. (
2021
).
Inhibitory stabilization and cortical computation
.
Nature Reviews Neuroscience
,
22
(
1
),
21
37
.
Safavi
,
S.
,
Chalk
,
M.
,
Logothetis
,
N.
, &
Levina
,
A.
(
2023
).
Signatures of criticality in efficient coding networks
.
bioRxiv:2023–02
.
Saxena
,
S.
, &
Cunningham
,
J. P
. (
2019
).
Towards the neural population doctrine
.
Current Opinion in Neurobiology
,
55
,
103
111
.
Schwemmer
,
M. A.
,
Fairhall
,
A. L.
,
Denéve
,
S.
, &
Shea-Brown
,
E. T
. (
2015
).
Constructing precisely computing networks with biophysical spiking neurons
.
Journal of Neuroscience
,
35
(
28
),
10112
10134
.
Sepulchre
,
R.
(
2022
).
Spiking control systems
.
Proceedings of the IEEE
,
110
(
5
).
Sepulchre
,
R.
,
Drion
,
G.
, &
Franci
,
A
. (
2019
).
Control across scales by positive and negative feedback
.
Annual Review of Control, Robotics, and Autonomous Systems
,
2
,
89
113
.
Seung
,
H. S
. (
1996
).
How the brain keeps the eyes still
.
Proceedings of the National Academy of Sciences
,
93
(
23
),
13339
13344
.
Seung
,
H. S.
,
Richardson
,
T.
,
Lagarias
,
J.
, &
Hopfield
,
J. J.
(
1997
).
Minimax and Hamiltonian dynamics of excitatory-inhibitory networks
. In
M. I.
Jordan
,
M. J.
Kearns
, &
S. A.
Solla
(Eds.)
,
Advances in neural information processing systems
,
10
.
MIT Press
.
Shadlen
,
M. N.
, &
Newsome
,
W. T
. (
1998
).
The variable discharge of cortical neurons: Implications for connectivity, computation, and information coding
.
Journal of Neuroscience
,
18
(
10
),
3870
3896
.
Shao
,
Y.
, &
Ostojic
,
S
. (
2023
).
Relating local connectivity and global dynamics in recurrent excitatory-inhibitory networks
.
PLOS Computational Biology
,
19
(
1
), e1010855.
Siahkamari
,
A.
,
Gangrade
,
A.
,
Kulis
,
B.
, &
Saligrama
,
V
. (
2020
).
Piecewise linear re- gression via a difference of convex functions
. In
Proceedings of the International Conference on Machine Learning
(pp.
8895
8904
).
Slijkhuis
,
F. S.
,
Keemink
,
S. W.
, &
Lanillos
,
P.
(
forthcoming
). Closed-form control with spike coding networks.
IEEE Transactions on Cognitive and Developmental Systems
.
Softky
,
W. R.
, &
Koch
,
C
. (
1993
).
The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs
.
Journal of Neuroscience
,
13
(
1
),
334
350
.
Song
,
H. F.
,
Yang
,
G. R.
, &
Wang
,
X.-J
. (
2016
).
Training excitatory-inhibitory recurrent neural networks for cognitive tasks: A simple and flexible framework
.
PLOS Computational Biology
,
12
(
2
), e1004792.
Sussillo
,
D
. (
2014
).
Neural circuits as computational dynamical systems
.
Current Opinion in Neurobiology
,
25
,
156
163
.
Thalmeier
,
D.
,
Uhlmann
,
M.
,
Kappen
,
H. J.
, &
Memmesheimer
,
R.-M
. (
2016
).
Learning universal computations with spikes
.
PLOS Computational Biology
,
12
(
6
), e1004895.
Van Vreeswijk
,
C.
, &
Sompolinsky
,
H
. (
1996
).
Chaos in neuronal networks with balanced excitatory and inhibitory activity
.
Science
,
274
(
5293
),
1724
1726
.
Vogels
,
T. P.
, &
Abbott
,
L. F
. (
2005
).
Signal propagation and logic gating in networks of integrate-and-fire neurons
.
Journal of Neuroscience
,
25
(
46
),
10786
10795
.
Vyas
,
S.
,
Golub
,
M. D.
,
Sussillo
,
D.
, &
Shenoy
,
K. V
. (
2020
).
Computation through neural population dynamics
.
Annual Review of Neuroscience
,
43
,
249
275
.
Wilson
,
H. R.
, &
Cowan
,
J. D
. (
1972
).
Excitatory and inhibitory interactions in localized populations of model neurons
.
Biophysical Journal
,
12
(
1
),
1
24
.
Yamins
,
D. L.
, &
DiCarlo
,
J. J
. (
2016
).
Using goal-driven deep learning models to understand sensory cortex
.
Nature Neuroscience
,
19
(
3
),
356
365
.
Yuille
,
A. L.
, &
Rangarajan
,
A
. (
2003
).
The concave-convex procedure
.
Neural Computation
,
15
(
4
),
915
936
.
Zenke
,
F.
, &
Vogels
,
T. P
. (
2021
).
The remarkable robustness of surrogate gradient learning for instilling complex function in spiking neural networks
.
Neural Computation
,
33
(
4
),
899
925
.
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. For a full description of the license, please visit https://creativecommons.org/licenses/by/4.0/legalcode