The geometry of robustness in spiking neural networks
Abstract
Neural systems are remarkably robust against various perturbations, a phenomenon that still requires a clear explanation. Here, we graphically illustrate how neural networks can become robust. We study spiking networks that generate lowdimensional representations, and we show that the neurons’ subthreshold voltages are confined to a convex region in a lowerdimensional voltage subspace, which we call a 'bounding box'. Any changes in network parameters (such as number of neurons, dimensionality of inputs, firing thresholds, synaptic weights, or transmission delays) can all be understood as deformations of this bounding box. Using these insights, we show that functionality is preserved as long as perturbations do not destroy the integrity of the bounding box. We suggest that the principles underlying robustness in these networks — lowdimensional representations, heterogeneity of tuning, and precise negative feedback — may be key to understanding the robustness of neural systems at the circuit level.
Editor's evaluation
The article introduces a geometrical interpretation for the dynamics and function of certain spiking networks, based on the earlier work of Machens and Deneve. Given that spiking networks are notoriously hard to understand, the approach could prove useful for many computational neuroscientists. Here, that visualization tool serves to assess how fragile the network is to perturbation of its parameters, such as neuronal death, or spurious noise in excitation and inhibition.
https://doi.org/10.7554/eLife.73276.sa0Introduction
The ability to maintain functionality despite perturbations is one of the defining properties of biological systems, from molecular signaling pathways to whole ecosystems (Csete and Doyle, 2002; Kitano, 2004; Whitacre, 2012; Félix and Barkoulas, 2015). Neural systems likewise withstand a certain amount of damage or external disturbances, which is evident from lesion studies or neurodegenerative diseases (Morrison and Hof, 1997; Bredesen et al., 2006; Palop et al., 2006), as well as from perturbation experiments (Wolff and Ölveczky, 2018; Li et al., 2016; Trouche et al., 2016; Fetsch et al., 2018). However, the mechanisms that underlie this robustness are not entirely clear. Indeed, most models of neural networks, when faced with partial damage, lose their functionality quite rapidly (Figure 1A–C; Seung et al., 2000; Koulakov et al., 2002; Li et al., 2016). Beyond its biological interest, understanding the robustness of neural systems is also crucial for the correct interpretation of experiments that seek to manipulate neural circuits (Wolff and Ölveczky, 2018).
Robustness has sometimes been attributed to various singlecell mechanisms, such as those that stabilize the dynamics of the stomatogastric ganglion of crustaceans against temperature fluctuations (O’Leary and Marder, 2016; Haddad and Marder, 2018), or the oculomotor integrator against instabilities in positive feedback loops (Koulakov et al., 2002; Goldman et al., 2003). On the circuitlevel, robustness has been tied to the excitatoryinhibitory (EI) balance of synaptic inputs, either by linking such balance with the efficiency of neural population codes (Boerlin et al., 2013; Bourdoukan et al., 2012; Barrett et al., 2013), or by using it as a corrective feedback for integrator models (Lim and Goldman, 2013). The corresponding spiking networks can maintain functional representations of their inputs by rebalancing when faced with perturbations such as neuron loss (Barrett et al., 2016; Lim and Goldman, 2013). Figure 1D illustrates such a robust, spiking network model.
Here, we illustrate how circuits can be made robust through simple geometric insights that tie the lowdimensional representations found in many population recordings (Saxena and Cunningham, 2019; Keemink and Machens, 2019; Vyas et al., 2020) to the biophysics of neurons, such as their voltages, thresholds, and synaptic inputs. We therefore provide a principled theory on how networks may have become robust to the many perturbations encountered in nature. We use this theory to illustrate two effects. First, we show that the resulting robustness mechanisms include the balanced regime, but are not limited to it. Indeed, networks can be robust without exhibiting any EI balance. Second, we predict a surprising asymmetry to perturbations: we find that robust networks are insensitive to broad inhibitory perturbations, yet quite sensitive to small excitatory perturbations, even if the latter are restricted to single neurons in large networks. This heightened sensitivity may explain the ability of animals to recognize exceedingly small, excitatory perturbations (Houweling and Brecht, 2008; Huber et al., 2008; Dalgleish et al., 2020).
Beyond questions of robustness, our work also provides a new way of thinking about spiking networks, which complements and extends classical approaches such as meanfield or attractor dynamics. To simplify our exposition, we focus on generic networks of integrateandfire neurons, rather than modeling a specific system. Consequently, we ignore part of the biological complexity (e.g. Dale’s law, more complex computations or dynamics), and defer explanations on how the framework may generalize to more realistic network models to the discussion and the methods.
Results
Our first assumption is that neural networks generate lowdimensional representations of sensory or motor signals, which can be extracted from the spike trains of a neural population through filtering and summation. Here, ‘lowdimensional’ simply means that the number of signals (or dimensions) represented is far less than the number of neurons in a circuit, so that there exists a certain amount of redundancy. Such redundant representations have been observed in many brain circuits (Saxena and Cunningham, 2019; Keemink and Machens, 2019), and are an integral part of most mid to largescale network models (Vogels et al., 2005; Eliasmith and Anderson, 2004; Barak, 2017). However, they do not per se guarantee robustness as shown in an example network in Figure 1C.
Passive redundancy
A classical example of a redundant, but nonrobust representation is a sensory layer of $N$ independent neurons acting as feature detectors. Here, each neuron receives the same $M$ timevarying input signals, $({x}_{1}(t),{x}_{2}(t),\dots ,{x}_{M}(t))=\mathbf{x}(t)$. Each signal is weighted at a neuron’s synapse, and the resulting synaptic currents are then summed in the neuron’s soma. If the synaptic weights are chosen differently for different neurons, the population generates a distributed code, whose redundancy we define as the number of neurons per signal dimension, or $\rho =N/M$. We will call this redundancy 'passive' as each neuron fires completely independent of what other neurons are doing.
While actual sensory systems are obviously more complex, this layer of independent feature detectors still serves as a useful baseline. For instance, in such a layer, perturbing a set of neurons by exciting or inhibiting them will have an effect on the representation that is directly proportional to the number of neurons perturbed. Passive redundancy therefore leads to a gradual decline of functionality (Figure 1B) or a gradual response of a system to any perturbation.
Autoencoder with lowdimensional readouts
To create robustness to perturbations, neurons cannot act independently, but rather need to coordinate their firing. We will now consider an example network that is generating a sensory representation at the population level, such that this representation is optimal with respect to a given linear readout (Figure 2A). We will focus on this simple scenario in order to highlight the mechanisms that endow networks with robustness. In the Discussion and Materials and methods, we will point out how to transfer these insights to more general networks.
Just as above, we consider a network of $N$ neurons which receive an $M$dimensional vector of input signals, $\mathbf{x}(t)$. The task of the network is to be an autoencoder, that is, to generate spike trains such that the input signals can be read out by a downstream area. We assume a linear readout, which filters each spike train with an exponential filter, similar to the postsynaptic potentials generated in a single synapse. Then, the filtered spike trains are weighted and summed, similar to the passive summation in a dendritic tree. Formally, we write
where ${r}_{k}(t)$ is the filtered spike train of the $k$th neuron, $N$ is the number of neurons, $\widehat{\mathbf{x}}(t)=({\widehat{x}}_{1}(t),{\widehat{x}}_{2}(t),\mathrm{\dots},{\widehat{x}}_{M}(t))$ is the vector of readouts, distinguished from the input signals by a hat, and ${\mathbf{D}}_{k}=({D}_{1k},{D}_{2k},\mathrm{\dots},{D}_{Mk})$ is the decoding vector of the $k$th neuron, whose individual elements contain the respective decoding weights.
We can depict the geometrical consequences of this decoding mechanism by imagining a network of five neurons that is encoding two signals. At a given point in time, we can illustrate both the input signals and the readout produced by the network as two points in signal space (Figure 2B, black cross and gray dot). Now let us imagine that one of the neurons, say neuron i, spikes. When that happens, the spike causes a jump in its filtered output spike train. In turn, and according to Equation 1, the vector of readouts, $\hat{\mathbf{x}}$, jumps in the direction of the decoding vector, $\mathbf{D}}_{i$, as illustrated in Figure 2B. Since the direction and magnitude of this jump are determined by the fixed readout weights, they are independent of the past spike history or the current values of the readouts. After this jump, and until another neuron fires, all components of the readout will decay. Geometrically, this decay corresponds to a movement of the readout towards the origin of the coordinate system.
Coordinated redundancy and the error bounding box
We furthermore assume that a neuron spikes only when its spike moves the readout closer to the desired signal, $\mathbf{x}$. For each neuron, this spike rule divides the whole signal space into two regions: a 'spike' halfspace where the readout error decreases if the neuron spikes, and a 'nospike' halfspace where the readout error increases if the neuron spikes (Figure 2B). The boundary between these two half spaces is the neuron’s spiking threshold, as seen in signal space. Consequently, the neuron’s voltage, ${V}_{i}$, must be at threshold, ${T}_{i}$, whenever the readout reaches this boundary, and the voltage must be below or above threshold on either side of it. We therefore identify the neuron’s voltage with the geometric projection of the readout error onto the decoding vector of the neuron,
where, without loss of generality, we have assumed that ${\mathbf{D}}_{i}$ has unit length (see Materials and methods, 'Coordinated spiking and the bounding box'). The effect of this definition is illustrated in Figure 2E, where the voltage increases or decreases with distance to the boundary. Accordingly, the voltage measures part of the error, given here by the distance of the readout to the neuron’s boundary.
In addition to its functional interpretation, the voltage equation has a simple biophysical interpretation, as illustrated in Figure 2C. Here, the two input signals, x_{1} and x_{2}, get weighted by two synaptic weights, ${D}_{1i}$ and ${D}_{2i}$, leading to two postsynaptic voltages that are then summed in the dendritic tree of neuron i. At the same time, the two readouts, ${\widehat{x}}_{1}$ and ${\widehat{x}}_{2}$, are fed back into the neuron via two exactly opposite synaptic weights, ${D}_{1i}$ and ${D}_{2i}$, thereby giving rise to the required subtraction. As a consequence, the neuron’s voltage becomes the projection of the readout error, as prescribed above. When the neuron’s voltage reaches the voltage threshold, ${T}_{i}$, the neuron fires a spike, which changes the readout, $\widehat{\mathbf{x}}$. In turn, this change is fed back into the neuron’s dendritic tree and leads to an effective reset of the voltage after a spike, as shown in Figure 2D.
One neuron alone can only improve the readout along one specific direction in signal space and thus cannot correct the readout for all possible input signals (Figure 2D, arrow). To properly bound the error, we therefore need several neurons, and each neuron needs to mix the input signals in different ways. As a consequence, the error will be corrected along different directions in signal space. A second neuron, say neuron j, is added in Figure 2F–H. Following the logic above, its voltage is given by ${V}_{j}={\mathbf{D}}_{j}^{\top}(\mathbf{x}\widehat{\mathbf{x}})$, and the respective voltage isoclines are shown in Figure 2H. We see that the voltage of neuron j jumps when neuron i spikes. Mathematically, the size of this jump is simply given by the dot product of the two decoding vectors, ${\mathbf{D}}_{j}^{\top}{\mathbf{D}}_{i}$. Biophysically, such a jump could be caused by negative feedback through the readout units, but it could also arise through a direct synaptic connection between the two neurons, in which case ${\mathrm{\Omega}}_{ji}={\mathbf{D}}_{j}^{\top}{\mathbf{D}}_{i}$ corresponds to the synaptic weight from neuron i to neuron j.
Finally, if we add three more neurons, and give them different sets of decoding weights, the network as a whole can restrict the readout to a bounded region in signal space (a polygon in two dimensions), as shown in Figure 2I–K. We will call this bounded region the 'error bounding box' or simply the 'bounding box.' Its overall size determines the error tolerance of the network. To highlight the structure of this network, we can change Equation 2 by inserting the definition of the readout, Equation 1, to obtain
Here, the term ${\mathrm{\Omega}}_{ik}={\mathbf{D}}_{i}^{\top}{\mathbf{D}}_{k}$ can be interpreted as a lateral connection between neurons i and k in the network (Figure 2I). The diagonal elements of the respective connectivity matrix, ${\mathrm{\Omega}}_{ii}$, can be interpreted as the hyperpolarization of the membrane voltage following a spike. While the connectivity of the network is symmetric, this assumption can be relaxed (see Materials and methods, 'Generalization of the bounding box II'). The connectivity term shows that information about the readout can be relayed through lateral connections and selfresets (Figure 2I), rather than through explicit negative feedback from a downstream layer. In either case, the feedback causes neurons to coordinate their firing. We will refer to this mechanism as 'coordinated spike coding' (Boahen, 2017) or 'coordinated redundancy'.
As shown previously (Bourdoukan et al., 2012; Boerlin et al., 2013), the temporal derivative of the above equation yields a network of currentbased, leaky integrateandfire neurons (see Materials and methods, 'Coordinated spiking and the bounding box'). We emphasize that there are two distinct situations that cause neurons to emit spikes. First, the readout always leaks towards the origin, and when it hits one of the boundaries, the appropriate neuron fires and resets the readout into the bounding box. Second, any change in the input signal, $\mathbf{x}$, causes a shift of the entire bounding box, since the signal is always at the centre of the box. A sudden shift may therefore cause the readout to fall outside of the box, in which case neurons whose boundaries have been crossed will fire to get the readout back into the box.
When the signal dimensionality is increased to $M=3$, the thresholds of the neurons become planes, and the bounding box is determined by the intersection of all planes, thereby becoming a threedimensional object such as a soccer ball. We strongly recommend to view Video 1 for an animation of the operation of the network in two and three dimensions, which highlights the relation between the bounding box and the resulting spike trains produced by the network.
The bounding box limits the coding error
The maximum coding error is limited by the size of the error bounding box, simply because the readout cannot deviate from the signal beyond the borders of the box. The size of the box is determined by the neurons’ thresholds. For simplicity, we will assume that all thresholds are identical, which could be regulated through homeostatic mechanisms (Turrigiano, 2012). The more general scenario is explained in Materials and methods, 'Generalization of the bounding box I'.
Beyond changing the maximum allowable coding error, the size of the error bounding box also influences the resulting code in more subtle ways. First, the coding error can be split into a systematic bias and random fluctuations (see Appendix 1—figure 1A). As the box becomes wider, the systematic bias increases. This bias can be largely eliminated by rescaling the readouts with a constant factor. We will sometimes use this corrected readout (see Materials and methods, 'Readout biases and corrections'), but note that the corrected readout is not confined to stay within the bounding box. Second, if the box becomes very narrow, the readout can eventually jump beyond the boundary of the opposite side, thereby crossing the threshold(s) of oppositely tuned neurons (see also Appendix 1—figure 1B–D). By default, we will assume that the bounding box is sufficiently wide to avoid this effect.
Finally, while the bounding box may seem like a fairly abstract construction, we note that it also has a simple, physical manifestation. Since the neurons’ voltages are constrained to live in an $M$dimensional subspace, a constraint given by the input dimensionality, the error bounding box delineates the borders of this voltage subspace, which is illustrated in Appendix 1—figure 2.
Robustness to inhibitory, sensitivity to excitatory perturbations
We will now study how the network reacts to perturbations by contrasting two example networks at opposite ends of a spectrum. In the first network, neurons are independent. For a twodimensional signal, we obtain this scenario when the bounding box consists of four neurons forming a square (Figure 3A, left), in which case neighbouring decoding vectors are orthogonal to each other, and their recurrent connections disappear (${\mathbf{D}}_{i}^{\top}{\mathbf{D}}_{k}=0$, see also Materials and methods, 'Generalization of the bounding box III'). The second network consists of $N=21$ randomly tuned neurons with equidistant thresholds, in which case the bounding box approximates the shape of a circle (Figure 3B, left).
The first perturbation we consider is the death of a single neuron. Throughout an organism’s life, cells, including neurons, can undergo the process of cell death or apoptosis if they are damaged or unfit, as may happen in diseased states (Moreno et al., 2015; Morrison and Hof, 1997; Bredesen et al., 2006; Coelho et al., 2018). Geometrically, the death of a neuron is equivalent to the removal of its corresponding face from the bounding box (Figure 3A and B, and Video 2). When the bounding box is breached on one side, the readout can no longer contain changes of the input signal along the open direction. This is precisely what happens in the case without redundancy (Figure 3A, right). In contrast, with coordinated redundancy, the removal of a single neuron has an almost imperceptible impact on the shape of the bounding box (Figure 3B, right). Consequently, the coding error remains bounded with essentially unchanged precision. The bounding box provides therefore a straightforward and intuitive explanation for the robustness against neuron loss observed in these spiking networks (Barrett et al., 2016; Boerlin et al., 2013).
The second perturbation we consider is a change in the excitability of one neuron. Such a change could come about through an experimentally injected current, for example, via patch clamp or optogenetics, or because of intrinsic plasticity. In either case, a change in excitability is equivalent to a change in the neuron’s threshold (see Materials and methods, 'Perturbations'). Within the bounding box picture, an inhibitory perturbation or decrease in excitability leads to an outward shift of the neuron’s threshold, and an excitatory perturbation or increase in excitability leads to an inward shift (Figure 3C and D). Without redundancy, the bounding box expands or shrinks, respectively (Figure 3C). At first sight, changing the box size increases or decreases the maximum error of the readout. More subtly, however, it also introduces a bias in the corrected average readout (Figure 3C, arrows). With coordinated redundancy, inhibitory and excitatory perturbations do not have opposite effects. Whereas an excitatory perturbation has an effect equivalent in size to a nonredundant system, as the threshold that is pushed inwards shrinks the box (Figure 3D, left), an inhibitory perturbation does not affect the system at all, because the outward shift of the threshold is compensated by the presence of other neurons (Figure 3D, right, see also Video 2). We note that a strong excitatory perturbation could cause the readout to move beyond the boundary of the opposite side, thereby leading to the undesirable 'pingpong' effect (Appendix 1—figure 1). This effect can be avoided if the bounding box is sufficiently wide.
In most experimental settings, one will perturb many neurons at once. We illustrate the effects of two such perturbations in Figure 3E and F. In Figure 3E, left, we excite six randomly selected neurons so that their boundaries are pushed inwards. As a consequence, the bounding box shrinks equally from all sides (Figure 3E, right). Just as in Figure 3D, left panel, the respective inward shifts cause (inputsignaldependent) biases in the corrected readout. In contrast, when we randomly inhibit neurons as in Figure 3E, right, the respective boundaries are pushed outwards. The bounding box barely changes, as the remaining neurons contain the readout error with essentially unchanged precision, and the system remains functional. In Figure 3F, we excite and inhibit neurons with similar tuning. In this case, the inhibitory perturbation is so large that there are no longer any neurons that can compensate. As a consequence, the bounding box expands in the perturbed direction and the corrected readout becomes biased, even for inhibitory perturbations.
In summary, we observe that coordinated redundancy endows the system with robustness against perturbations that act inhibitorily on the neurons, such as neuron death or increases in spiking thresholds. The function of the system remains intact until the bounding box breaks open. However, coordinated redundancy also makes the system highly sensitive towards any excitatory perturbations or decreases in spiking thresholds. Indeed, perturbing only a single neuron is sufficient to generate a change in the readout. These results contrast with passively redundant systems, whose representation changes gradually with either excitatory or inhibitory perturbations.
We note that there is circumstantial evidence that cortical systems are indeed highly sensitive to excitatory perturbations, and potentially less sensitive to inhibitory perturbations. For instance, animals can detect excitatory currents injected into a few pyramidal cells in somatosensory cortex (Houweling and Brecht, 2008; Huber et al., 2008; Dalgleish et al., 2020). To our knowledge, no study has shown that cortical systems are similarly sensitive to inhibitory perturbations in one or a few neurons. Rather, several studies have shown that neural systems can compensate for inhibitory perturbations in a sizeable fraction of pyramidal cells (Li et al., 2016; Fetsch et al., 2018; Trouche et al., 2016).
The neurophysiological signatures of perturbations
Besides insights into a system’s functionality, the bounding box also allows us to see immediately how perturbations affect the firing of the unperturbed neurons. For instance, when we excite a single neuron, its threshold moves inwards (Figure 3D) and occludes the thresholds of the neighboring neurons, that is, neurons with similar selectivity. Since the readout can no longer reach these neurons, they stop firing, as shown in a simulation of the network in Figure 4A. Conversely, if we inhibit one or more neurons, their thresholds become hidden and they no longer participate in containing the readout (Figure 3D–F). As a consequence, the surrounding neurons have to pick up the bill and fire more, so that their firing rates will increase, as shown in Figure 4A. Biophysically, these effects are of course mediated by an increase or decrease of lateral or recurrent inhibition. However, the bounding box provides a simple visualization of the effect and its purpose.
The bounding box also allows us to visualize the relation between EI balance and robustness. We can again study two extreme examples, as illustrated in Figure 4B and C. Here, both bounding boxes are intact and contain the readout. The box in Figure 4B has low redundancy ($N=5$ neurons for $M=2$ signals), whereas the box in Figure 4C has high redundancy ($N=100$ neurons for $M=2$ signals). In turn, we can visualize the synaptic current balance in these boxes by highlighting how the movements of the readout give rise to excitatory or inhibitory currents.
In the low redundancy case, the readout initially decays, and thereby moves towards the threshold of one of the neurons (shown in green). A movement towards threshold corresponds to a depolarization of the respective membrane potential and therefore an excitatory (or disinhibitory) current, here illustrated in red. The respective neuron spikes, and the readout jumps away from the threshold, which is mediated by a hyperpolarizing selfreset current, here illustrated in green. Since the excitatory drive is cancelled by selfresets following spiking, the synaptic inputs are not balanced, and the neuron fires spikes in a relatively regular rhythm (Figure 4B, right panel, see Materials and methods, 'Excitationinhibition balance').
The situation is quite different for the highredundancy network (Figure 4C). Here, the readout decays towards the thresholds of multiple neurons (some of which are highlighted with gray lines), but only one of these neurons will fire. When that happens, all neurons in the vicinity immediately receive inhibitory currents that signal the concomitant change in the readout. These inhibitory currents thereby cancel the excitatory feedforward drive, and the respective neurons experience a tight EI balance, leading to sparse and irregular spike trains (Figure 4C, left panel). We note that this tight balance only holds in the neurons that are sufficiently close to the neuron whose threshold is crossed. Neurons further away will experience only small current fluctuations, or will, on average, be hyperpolarized.
As a result, we see that EI balance occurs in networks that are sufficiently redundant, but not in networks with no or low redundancy. Nonetheless, even networks with low redundancy have a measure of robustness: for instance, the network in Figure 4B is robust against the loss of one neuron. While previous work has suggested that networks recover functionality when perturbed by dynamically reestablishing EI balance (Lim and Goldman, 2013; Boerlin et al., 2013; Barrett et al., 2016), our considerations here show that robustness extends beyond the regime of EI balance. Figure 4D illustrates this result by contrasting the performance and balance of networks as a function of their redundancy.
Scaling up
While the simple toy networks we have studied so far are useful for illustration and intuition, biological neural networks, and especially cortical networks, consist of thousands of neurons that are thought to represent hundreds of signals simultaneously. To get closer to the biological reality, we therefore also need to study larger and more powerful networks. As we have shown above, the features of networks are tightly linked to the shape of the bounding box. For threedimensional input signals, the threshold of each neuron becomes a plane, and the bounding box becomes a polyhedron (see also Appendix 1—figure 3A and Video 1). For higher dimensional signals, the precise shape of the bounding box is hard to visualize. However, if we assume that the number of neurons scales linearly with the number of signal dimensions, $N=\rho M$, one can show that the resulting, higher dimensional bounding boxes are somewhat closer to a hypercube than a hypersphere. (Some insights on the geometry of such higher dimensional bounding boxes can be found in Appendix 1—figure 3.)
In Figure 3, we saw that all bounding boxes are sensitive to excitatory perturbations, but that only noredundancy (or very lowredundancy) bounding boxes are sensitive to inhibitory perturbations. A key question when scaling up is therefore whether larger networks with finite redundancy become sensitive to inhibitory perturbations, or whether they remain insensitive. An extreme form of inhibitory perturbation is the loss of neurons. Figure 5 shows that, even in high dimensions, networks are robust against such perturbations. As before, we assume that the decoding vectors of the neurons ${\mathbf{D}}_{i}$ are of similar length, but otherwise random, and that the thresholds of all neurons are the same. As the death or birth of random neurons simply corresponds to a change in the overall redundancy of the network, we can understand how network performance and statistics change by simply moving along the redundancy axis in Figure 5A–C. We observe that changing the redundancy over a broad range ($\rho \approx 550$) has negligible effects on the performance (Figure 5A). This contrasts with networks of independent neurons in which performance scales linearly with any change in redundancy for a fixed readout. We furthermore observe that decreasing redundancy, leads to higher firing rates (Figure 5B) and more regular firing, or lower CVs (Figure 5C). The decrease of CVs here simply reflects a decrease in the number of spike patterns that can represent the constant input signal, given the smaller pool of neurons in the network. In other words, when we kill neurons, the neural code becomes less redundant, and the spike patterns of individual neurons lose some of their apparent randomness. Conversely, as the network size increases, so does the number of possible spike patterns, with the consequent increase of CV. As the number of neurons keeps increasing, it becomes more and more likely that the network has neurons optimally tuned to a given input signal, contributing to a decrease of the CV. Therefore, the increase and subsequent decrease in CV with increasing redundancy is the result of these two counteracting effects (Figure 5C).
In summary, when we scale up, networks with some redundancy remain robust to partial, inhibitory perturbations, even though the firing statistics of the neurons change.
Natural perturbations
Biological systems should also be robust against the mistuning of any of their components. We will now show that many types of parameter mistuning can be understood as deformations of the bounding box. As shown in Figure 3, the simplest type of perturbation is a change in a neuron’s spiking threshold: an increase of a neuron’s spiking threshold will push the corresponding face of the bounding box outwards, and a decrease will push the face inwards.
While permanent changes in the threshold can come about through changes in conductances or reversal potentials, a neuron can also suffer from temporary changes in its effective spiking threshold through, for example, noise. Biological systems are constantly subject to noise at multiple levels such as sensory transduction noise, ion channel noise (Faisal et al., 2008), or 'background' synaptic activity (Destexhe et al., 2001; Fellous et al., 2003). We can study the impact of such noise by injecting small, random currents into each neuron. These currents change how close the voltage of a neuron is to its spiking threshold. With regard to spike generation, the resulting voltage fluctuations are thus equivalent to fluctuations of the threshold, or random movements of all of the faces of the bounding box around their unperturbed positions (Figure 6A, see also Video 2).
For networks with low redundancy, $\rho $, small voltage fluctuations cause only minor deformations of the bounding box. In turn, the error tolerance remains roughly the same, and network performance is not affected (Figure 6B, left; Figure 6E and F). However, for networks with high redundancy, $\rho $, small voltage fluctuations can cause a fatal collapse of the system (Figure 6B, middle). The key reason is that the effective size of the bounding box is determined by the position of the thresholds that have moved furthest into the box. As more and more neurons are added, the likelihood that some of them have very decreased thresholds increases, and the effective size of the bounding box shrinks. In turn, the probability that the network moves into an 'epileptic seizure' (due to the 'pingpong' effect, see Appendix 1—figure 1) increases as well. While the readouts may still be contained in this scenario (Figure 6E), the excessive number of spikes fired (Figure 6F) comes at a high metabolic cost and would be detrimental to biological systems. To avoid this failure mode, neurons need to lower their excitability, which in turn increases the size of the bounding box for a fixed redundancy (Figure 6B, right panel). Such a 'wide box' will be more resilient towards noise (Figure 6B, right panel, Figure 6E and F). More generally, our results suggest that more redundant networks may require a better control or suppression of intrinsic sources of noise than less redundant networks.
Next, we will study perturbations of a neuron’s reset potential, that is, the voltage reached directly after a spike. This voltage should ideally be ${V}_{i,\text{reset}}={T}_{i}{\mathbf{D}}_{i}^{\top}{\mathbf{D}}_{i}$. Biophysically, when the neuron resets to a voltage above (below) this ideal reset potential, then its postspike voltage is temporarily closer (further) from threshold. In terms of the neuron’s spiking output, a change in its reset voltage is therefore equivalent to a (temporary) change in its threshold. Within the bounding box, a reset voltage above (below) the optimal reset will lead to a push of the neuron’s threshold inwards (outwards) (Figure 6C). However, because of the voltage leak, the threshold will then decay back to its normal position. Video 2 illustrates this effect in a system with a twodimensional input. We note that positive and negative changes to the default reset potential will lead to asymmetric effects on robustness like those observed for excitatory and inhibitory perturbations. Specifically, if the resets become too small, and if the leak is insufficiently fast, then successive spiking of a single neuron will draw its threshold inwards, thereby leading to a collapse of the bounding box.
Finally, we study perturbations of the synaptic connectivity in the network. Synapses could be permanently mistuned or they could be temporarily mistuned, for instance through transmission failures or through stochastic fluctuations in the release of neurotransmitters (Faisal et al., 2008). From a geometric perspective, a mistuned synapse causes a temporary change in the threshold of the postsynaptic neuron whenever a presynaptic spike arrives (Figure 6D). We again note an asymmetry: an excitatory synapse with decreased strength (or an inhibitory synapse with increased strength) leads to an outward move of the postsynaptic neuron’s threshold, which is generally harmless. In turn, an excitatory synapse with increased strength (or an inhibitory synapse with decreased strength) leads to an inward move, which could be a temporarily harmful perturbation. Accordingly, random synaptic failures in excitatory synapses (but not inhibitory synapses) leave the bounding box functionally intact.
When all synapses in the network are randomly mistuned, then each spike fired will cause a random, but transient deformation of the bounding box (see Video 2). Overall, we find that more redundant networks (with consequently more synapses) are typically more vulnerable to these perturbations. Just as for voltage noise, the amount of deformation of the bounding box therefore increases with the number of neurons. For large perturbations, the synaptic noise eventually leads to inefficient networks with high spike rate (Figure 6G and H). As shown in Appendix 1—figure 4, the effects of (voltage or synaptic) noise on the networks hold independent of the signal dimensionality.
Synaptic delays
So far, we have assumed that the propagation of action potentials is instantaneous. However, lateral excitation and inhibition in biological networks incur delays on the order of milliseconds. Previous work has shown that networks which coordinate their spiking as suggested here are extremely sensitive to delays when neurons are similarly tuned (Chalk et al., 2016; Rullán Buxó and Pillow, 2020). Indeed, when spikes are delayed, voltages no longer reflect an accurate estimate of the coding error. For neurons with identical decoders, delays can lead to uninformed spikes that actually increase the coding error (Figure 7A and B). With networks that represent Mdimensional signals at once, the effects of delays are more complex. However, the bounding box allows us to visualize them and explain how they can, in principle, be avoided. Below, we study the impact of these delays, which apply directly to recurrent excitation and inhibition. We also apply the same delays to the network readout for mathematical convenience, but those do not affect the network dynamics (see Materials and methods).
To visualize the effect of a synaptic delay, we show the readout dynamics around a single spike in Figure 7C (see also Video 2). After hitting the threshold, the spiking neuron resets its own voltage immediately. However, due to the delay, neither a hypothetical readout unit nor other neurons in the network are aware of the spike. From the network perspective, the voltage of the spiking neuron appears temporarily too low (or its threshold too high), which we can visualize as an outward jump of its boundary (Figure 7C, second and third panels). When the spike finally arrives, the readout and all other voltages are updated, and the voltage of the firing neuron once again agrees with the network state. In our visualization, its boundary thus returns to its default position (Figure 7C, fourth panel).
The visualization illustrates that the effect of a delayed spike depends on the bounding box shape. In Figure 7C, the nearly orthogonal tuning of neighbouring neurons makes the delayed spike harmless. The situation is different if neurons are more similarly tuned as in Figure 7D. Here, a second neuron might cross its threshold before the delayed spike arrives. As a consequence, it also fires a spike, and its boundary also retracts (Figure 7D, third panel). Eventually, both spikes arrive, the readout is updated, and the bounding box regains its original shape (Figure 7D, fourth panel). At this point, the readout may overshoot, cross an opposite boundary, and trigger further 'pingpong' spikes. The resulting 'epileptic seizures' are essentially unavoidable in highly redundant networks with synaptic delays. Consequently, we identify two problems with synaptic delays. The first problem is that multiple thresholds are crossed simultaneously. The second problem is that the resulting strong change in the readout can cause a 'pingpong' of uninformed spikes.
For the first problem, we note that the impact of synaptic delays depends on the angles of neighbouring neurons, as shown in Figure 7C and D. For higher signal dimensions and fixed redundancy, these angles become more orthogonal (Appendix 1—figure 3D), which alleviates the detrimental effect of delays. Numerically, however, we find that this effect is not sufficient to avoid all uninformed spikes (for signal spaces with up to $M=50$ dimensions), and the networks still degenerate into ‘pingpong’.
To avoid the second problem, we need to eliminate the crossing of opposite thresholds by widening the box, which can prevent 'pingpong' (Figure 7E, Appendix 1—figure 1). However, permanently widening the bounding box in all directions can reduce coding accuracy, even when the readout is properly rescaled (Figure 7G, Appendix 1—figure 5, see Materials and methods, 'Iterative adaptation of parameters to avoid pingpong'). A better solution is therefore to 'widen' the box only temporarily. For instance, if we eliminate excitatory connections between pairs of neurons that are practically antipodes, we are setting the respective synapses to zero. This change in the network connectivity can be understood as a specific and targeted synaptic perturbation (Figure 6D), whose effect is to expand the thresholds of a neuron’s antipodes whenever it fires a spike, thereby temporarily widening the box exactly in the direction in which the readout overshoots (Figure 7F). As a consequence, the networks become less likely to initiate pingpong. Moreover, as their widening is local and only temporary, performance is less affected. Indeed, for higher dimensional systems and biologically plausible delays (1–2 ms), performance of networks with delays reaches the performance of networks without delays (Appendix 1—figure 5). The rapid increase in firing due to pingpong is avoided as well (see also Appendix 1—figure 5).
Adding computations
The bounding box provides a useful tool even if we endow the networks with a set of slower connections to perform linear or nonlinear computations (Boerlin et al., 2013; Savin and Deneve, 2014; Thalmeier et al., 2016). Indeed, the simulation in Figure 1D used these slower connections to generate oscillatory dynamics (see Materials and methods, 'Generalization of the bounding box IV'). This extension to networks that generate persistent activity or dynamical patterns works because the mechanisms underlying the encoding of the signals into spike trains are decoupled from the mechanisms that generate the dynamics of the signals (or readouts). In fact, the extra currents generated by the slow recurrent connections can be seen as a perturbation of the bounding box thresholds. This perturbation shifts the bounding box in the space of readouts as illustrated in Appendix 1—figure 6.
Discussion
In this study, we characterized the functioning of networks with coordinated redundancy under normal conditions and under a diversity of perturbations, using a simple, geometric visualization, the bounding box. The bounding box delimits the error that such a network tolerates in approximating a set of input signals, and its geometry is found to be largely determined by the properties of its decoders. It allows us to visualize and thus understand the dynamics of coordinated spiking networks, including the firing of every single spike. We showed how various perturbations of the network can be mapped onto shape deformations of this bounding box. As long as the box stays intact, the network’s performance is essentially unaffected, in that downstream readouts of the network’s outputs will not notice the perturbation.
In many respects, the bounding box is a 'toy model' (a deliberately abstract model), which we see mainly as a tool to conceptualize and highlight generic circuit mechanisms, rather than an attempt to model any specific system. Nonetheless, it is worthwhile to point out that the bounding box is also a spiking version of classical sparse coding models of V1 (Olshausen and Field, 1996). Indeed, previous work has demonstrated that these networks can explain various perturbation experiments in V1 (Barrett et al., 2016). So, besides shedding light on the robustness of coordinated spike codes, the bounding box can also be seen as a simple model of a sensory system.
Robustness of networks with coordinated spike coding
Several overarching principles have been identified that allow systems to be robust (Csete and Doyle, 2002; Kitano, 2004; Whitacre, 2012; Félix and Barkoulas, 2015). These include (1) negative feedback, to correct perturbations and recover functionality; (2) heterogeneity of components, to avoid common modes of failure; and (3) modularity or 'bowtie' architectures, to create alternative pathways or solutions in the case of a perturbation. Furthermore, (4) making a system robust against certain perturbations almost always involves a tradeoff, in that the system becomes fragile against other perturbations.
These core themes can also be found in the networks we studied here. (1) Negative feedback exists through extensive lateral connectivity (or, alternatively, through actual feedback of the readout, as in Figure 2F), and is precisely tuned such that it automatically corrects any perturbations. (2) Individual neurons are heterogeneous and thereby allow the system (as visualized by the bounding box) to remain functional for all types of input signals. (3) Since neuron space is always larger than signal space, there are many alternative neural codes ('alternative pathways') that give rise to the same linear readout, thus embodying a bowtie architecture whose core is the signaling space. Shrinking the network’s redundancy, for example, by killing neurons, in turn eliminates these alternative codes and leads to more regular and reliable spike trains. (4) Furthermore, the networks are fragile against any perturbation that leads to a shrinking of the box. Paradoxically, this fragility may become more relevant if a system becomes more redundant. These four themes may relate the robustness of the networks studied here to the more general topic of tissue robustness (Kitano, 2004).
Coordinated redundancy allows the construction of robust subcircuits, that can selfcorrect problems instead of passing them on, so that downstream networks remain unaffected. These observations remain correct even if we move beyond the simple autoencoder networks that we have studied here. Indeed, we could generalize the connectivities we consider or abandon the idea that the readout must match the input without changing the robustness of the networks (see Materials and methods, 'Coordinated spiking and the bounding box'). We could also add slower recurrent synapses which allows to generate dynamics within the networks (Boerlin et al., 2013; Savin and Deneve, 2014; Thalmeier et al., 2016), as explained above.
Fragility of networks with coordinated spike coding
Despite their strong robustness, networks with coordinated redundancy are also surprisingly fragile against any perturbations that cause an effective shrinking of the box, and thereby lead to the pingpong effect. These problems can be ameliorated by widening the box, which brings networks back into workable regimes if they represent highdimensional signals with limited redundancy. However, a true 'fix' of this problem can only be achieved if neurons with opposite decoder weights (which are connected through excitatory connections) are prohibited. Such a change would break the symmetric treatment of excitatory and inhibitory connections, which causes neurons to both excite and inhibit different downstream partners, thereby violating Dale’s law. Future work will need to reconsider these issues which seem to be tightly connected. (We note that Boerlin et al., 2013 developed networks that obey Dale’s law, but did so without fixing the issue of the pingpong effect.)
Structural robustness of neural networks
Historically, the study of networklevel mechanisms of robustness has received relatively little attention. A key focus has been the robustness of network attractors, defined as the ability of a system to remain in the same attractor landscape despite perturbations. For instance, systems such as the oculomotor integrator or the head direction system can be described as continuous attractors (Seung, 1996; Zhang, 1996). Such continuous attractors are structurally unstable, in that even small perturbations in single neurons can lead to rapid dynamic drifts (Seung, 1996; Zhang, 1996). However, this fragility to perturbations is not observed in biological neural networks.
In order to achieve the required robustness, several biophysical mechanisms have been proposed to enhance continuous attractors models, e.g. bistability at the somatic level (Koulakov et al., 2002) or dendritic level (Goldman et al., 2003). More recent work proposed networklevel mechanisms based on derivative feedback, in order to solve the problem of robustness for continuous attractor networks (Lim and Goldman, 2013). In our work, the problem is solved because perturbations such as neuron loss, noise, or tuning of synapses are compensated through the fast, lateral connections. As a consequence, perturbations of the singleneuron level (spiking) are uncoupled from perturbations of the populationlevel (readout). Consequently, only perturbations that manage to disturb the linear readout can impact the network attractor dynamics.
Models of neural networks implementing point attractors, such as the Hopfield model (Hopfield, 1982), are typically considered structurally robust, meaning that perturbations up to certain magnitudes of their parameters and the introduction of dynamics noise do not disrupt the attractor. We note, however, that perturbations in these networks lead to changes in neurons’ firing rates, which may still cause changes in putative downstream linear readouts. From the point of view of a downstream observer, perturbations are therefore not compensated within classical attractor networks. The induced perturbations may be inconsequential, however, when the downstream readout is taken to be a classifier; only the combined system of attractor network and classifier readout can then be seen as a 'robust module', that is, a module that keeps problems to itself, rather than spreading them to all who listen.
Similar observations apply to studies of the robustness of deep networks against various perturbations such as the loss of neurons (Morcos, 2018; Barrett et al., 2019). In these cases, the network’s robustness is evaluated with respect to the output of a final classification step, such as the identification of an object. Indeed, a lot of work has been dedicated to making this final output robust to small perturbations, especially perturbations applied to the inputs (Szegedy, 2013; Biggio, 2013; Carlini, 2019; Brendel et al., 2020). Based on the arguments above, we similarly expect that the problem of making a graded output robust will be harder and fundamentally different.
Insights on spiking networks
Spiking networks have traditionally been quite hard to understand, except for special cases (Maass and Bishop, 1999; Vogels et al., 2005; Gerstner et al., 2014). Here, we have shown how the dynamics of (coordinated) spike coding networks can be understood within a lowerdimensional signal space, which is tightly linked to linear readouts. Since (lowdimensional) linear readouts are a ubiquitous finding in recordings from neural populations, we may speculate that our signal space is roughly equivalent to the latent subspaces discovered by linear projections of neural activities, as, for example, obtained through dimensionality reduction methods (Cunningham and Yu, 2014; Keemink and Machens, 2019). This link between a space of neural activities and a space of (latent) signals is common to all network models based on lowrank connectivities (Eliasmith, 2005; Seung, 1996; Mastrogiuseppe and Ostojic, 2018). In contrast to these studies, however, and in line with (Boerlin et al., 2013), our work focuses on spiking networks and introduces a third space, the voltage space, which represents the system’s coding errors. As we have shown here, the coding errors are confined to an error bounding box. Accordingly, the bounding box finds its physical—and in principle measurable—manifestation in a lowdimensional subspace of a network’s voltage space (Appendix 1—figure 2).
We believe that the links we have made here—which allow us to jointly visualize a lowdimensional signal space, the spiking activity, and the subthreshold voltages—may provide useful insights into the functioning of spiking networks in the brain, and may well be expanded beyond the confines of the current study.
Materials and methods
Coordinated spiking and the bounding box
Request a detailed protocolMathematically, our networks can be derived from a single objective function that quantifies coding accuracy. Stepbystep derivation for the autoencoder networks can be found in Barrett et al., 2016; networks that additionally involve a set of slow connections are derived in Boerlin et al., 2013. Here, we focus on the autoencoder networks which contain all the crucial elements needed to understand the spiking dynamics of the networks. Instead of starting with an objective function, we take a slightly different perspective in our derivation here, which ties more directly into our geometric interpretations.
In short, we assume that a network of $N$ neurons encodes an $M$dimensional input signal $\mathbf{x}(t)$, in its spike trains $\mathbf{s}(t)$, such that the signal can be read out from the filtered spike trains,
Here, $\hat{\mathbf{x}}(t)$ is the $M$dimensional linear readout or signal estimate, the $M\times N$ matrix $\mathbf{D}$ contains the decoding weights (and each column corresponds to a decoding vector $\mathbf{D}}_{i$), the filtered spike trains are represented by an $N$dimensional vector $\mathbf{r}(t)$, and $\lambda $ determines the filtering time constant.
The key idea of coordinated spike coding is to derive a spiking rule that bounds the difference between the input signal $\mathbf{x}$, and the linear readout $\hat{\mathbf{x}}$,
where $\parallel \cdot \parallel $ denotes the Euclidean distance or L2 norm and $T$ determines the maximally allowed difference. In the network implementation, we approximate this bound (which defines a hypersphere) by a set of linear bounds or inequalities, one for each neuron i,
For simplicity, we assume that the decoding vectors $\mathbf{D}}_{i$ have unit norm. Each inequality defines a halfspace of solutions for the readout $\hat{\mathbf{x}}$. For properly chosen $\mathbf{D}}_{i$, the intersection of all of these halfspaces is nonempty and bounded, and thus forms the interior of the bounding box. Geometrically, the equations define a polytope $B=\{\hat{\mathbf{x}}\in {\mathbb{R}}^{M}\phantom{\rule{thickmathspace}{0ex}}{\mathbf{D}}^{\phantom{\rule{0.083em}{0ex}}\mathsf{T}}\phantom{\rule{negativethinmathspace}{0ex}}\left(\mathbf{x}\hat{\mathbf{x}}\right)<\mathbf{T}\}$. If the thresholds are chosen sufficiently large, then crossing a bound and firing a spike keeps the readout inside the bounding box.
The dynamics of the network are obtained by identifying the lefthand side of the above equation with the neuron’s voltage, ${V}_{i}$, and then taking the temporal derivative (Boerlin et al., 2013; Barrett et al., 2016). If we also add some noise to the resulting equations, we obtain,
which describes a network of leaky integrateandfire neurons. The first term on the righthand side is the leak, the second term corresponds to the feedforward input signals to the network, the third term captures the fast recurrent connectivity, with synaptic weights $\mathrm{\Omega}}_{ij}={\mathbf{D}}_{i}^{\phantom{\rule{0.083em}{0ex}}\mathsf{T}}{\mathbf{D}}_{j$, and the fourth term is added white current noise with standard deviation ${\sigma}_{V}$. When the voltage ${V}_{i}$ reaches the threshold $T$, the selfconnection $\mathrm{\Omega}}_{ii}={\mathbf{D}}_{i}^{\phantom{\rule{0.083em}{0ex}}\mathsf{T}}{\mathbf{D}}_{i$ causes a reset of the voltage to ${V}_{\text{reset}}=T+{\mathrm{\Omega}}_{ii}$. For biological plausibility, we also consider a small refractory period of ${\tau}_{\text{ref}}=2\text{ms}$ for each neuron. We implemented this refractory period by simply omitting any spikes coming from the saturated neuron during this period.
Mathematically, the voltages are thereby confined to a subspace given by the image of the transposed decoder matrix, ${\mathbf{D}}^{\top}$. The dynamics within this voltage subspace are then bounded according to Equation 7, which can be seen as a physical manifestation of the bounding box (see also Appendix 1—figure 2).
Generalization of the bounding box I: Heterogeneous thresholds
Request a detailed protocolIn our exposition, we generally assume that all decoding vectors are of the same length, and all thresholds are identical. For isotropically distributed input signals and isotropically distributed decoding vectors, this scenario will cause all neurons to fire the same average number of spikes over time. Indeed, to the extent that homeostatic plasticity sets synaptic weights and firing thresholds to guarantee this outcome (Turrigiano, 2012), a network will automatically revert to a spherically symmetric bounding box for such input signals (see also Brendel et al., 2020).
If input signals are not isotropically distributed then homeostatic plasticity would essentially lower the thresholds of neurons that receive overall less inputs, and it would increase the thresholds of neurons that receive overall more inputs. In turn, the bounding box would take more elliptical shapes. We have not considered this scenario here for simplicity, but the key findings on robustness will hold in this case, as well.
Generalization of the bounding box II: Asymmetric connectivities
Request a detailed protocolIn the main text, we have assumed that the readout always jumps orthogonal to the threshold boundary (or face) of a neuron. This assumption leads to symmetric connectivities in the network, given by $\mathrm{\Omega}}_{ij}={\mathbf{D}}_{i}^{\phantom{\rule{0.083em}{0ex}}\mathsf{T}}{\mathbf{D}}_{j$. However, our results on robustness also hold if we decouple the orientation of a neuron’s face from the direction of the readout jump. This can be achieved if we define the voltage as ${V}_{i}={\mathbf{F}}_{i}(\mathbf{x}\hat{\mathbf{x}})$, where $\mathbf{F}}_{i$ denotes the norm vector of a bounding box face, but then let the readout jump in the direction $\mathbf{D}}_{i$. A nonorthogonal jump with respect to the face then simply requires $\mathbf{D}}_{i}\ne {\mathbf{F}}_{i$. Indeed, for elliptically shaped bounding boxes, nonorthogonal jumps of the readout can be advantageous. The more general dynamical equation for the networks is then given by
and was first described in Brendel et al., 2020. In principle, these generalized networks include all spiking networks with lowrank connectivities. However, the bounding box interpretation is most useful when each spike is reset back into the bounding box, which will only happen if the net effect of a spike on neighboring neurons is inhibitory. Spikes that cause (temporary) jumps out of the box, and therefore have a net excitatory and erroramplifying effect, will be considered in future work.
Generalization of the bounding box III: Opening the box
Request a detailed protocolThe equation for the synaptic connectivity, $\mathrm{\Omega}}_{ij}={\mathbf{D}}_{i}^{\phantom{\rule{0.083em}{0ex}}\mathsf{T}}{\mathbf{D}}_{j$, implies that neurons with similar decoding vectors inhibit each other, neurons with orthogonal decoding vectors are unconnected, and neurons with opposite decoding vectors excite each other. Consequently, if the bounding box is a (hyper)cube, then almost all neurons are unconnected, except for neurons whose faces are opposite to each other. The excitatory connections between these neurons ensure that their voltages remain in sync. However, in practice, those voltages do not need to be tied, and the excitatory connections can therefore also be eliminated (as in Figure 7), which can help against the pingpong effect.
Alternatively, we can choose decoding vectors such that all synapses are inhibitory, ${\mathrm{\Omega}}_{ij}\le 0$. In this case, the bounding box remains open on one side. The network no longer represents the input signal, but rather computes a piecewise linear function of the input (Mancoo, 2020). In turn, the network’s new function (piecewise linear output) will now remain robust against perturbations for exactly the same reasons explained before. Indeed, the reader may notice that most of the results on robustness do not require the bounding box to be closed.
Generalization of the bounding box IV: Slow connections
Request a detailed protocolThroughout the manuscript, we focused on autoencoder networks. However, as illustrated in Figure 1 and derived in Boerlin et al., 2013, by introducing a second set of slower connections, we can endue these networks with computations,
in which case the network approximates the dynamical system:
We note that all results on robustness hold for these more complicated networks as well. Indeed, the robustness of the autoencoder networks relies on the fast recurrent connections, which are present in these architectures as well. Due to the time scale separation, these mechanisms do not interfere with the slower recurrent connections, which create the slow dynamics of the readouts (see also Appendix 1—figure 6).
Readout biases and corrections
Request a detailed protocolWhen one of the neurons fires, its spike changes the readout, which jumps into the bounding box. In previous work (Boerlin et al., 2013, Barrett et al., 2016), the neurons’ thresholds were linked with the length of the jumps through the equation ${T}_{i}=\Vert {\mathbf{D}}_{i}{\Vert}^{2}/2$. Accordingly, the jumps were generally taken to reach the opposing face of the bounding box, creating a tight error bounding box around $\mathbf{x}$. This setting guarantees that the timeaveraged readout matches the input signal.
When the jumps are significantly shorter than the average bounding box width, however, the timeaveraged readout will be biased away from the input signal (see Appendix 1—figure 1). However, in many cases, this bias can be corrected by rescaling the readout. For instance, if the bounding box is shaped like a hypersphere (i.e. in the limit of an infinite number of neurons $N$), and assuming a constant (or slowly varying) stimulus, we can correct the readout as
where the angular brackets denote the timeaveraged readout. Accordingly, in this case the bias only affects the length of the readout vectors, but not their direction.
If the bounding box is shaped like a hypercube, we alternatively correct the readout bias by assuming that a downstream decoder area has access to the identity of spiking neurons in the recent past. In this case, the downstream area can simply correct the readout according to the following equation:
where $S$ is the set of active neurons for a given fixed time window in the past.
In all other cases, we empirically found that we can apply a correction to the readout using a similar scaling as in Equation 11 where $\u27e8\Vert \mathbf{D}\mathbf{r}\Vert \u27e9\approx \Vert \mathbf{D}\mathbf{r}(t)\Vert$. In other words, in most cases, the bias mainly affects the length of the readout vectors, whereas their direction is less affected.
In Figure 1, we used networks that involve an extra set of slow recurrent connections (Boerlin et al., 2013). In this case, we additionally scaled the slow recurrent connectivity matrix ${\mathrm{\Omega}}_{\text{slow}}$ with the same scaling factor as the corrected readout in Equation 11:
Geometry of highdimensional bounding boxes
Request a detailed protocolThe dimensionality of the bounding box is determined by the dimensionality $M$ of the input signal. Throughout the illustrations in the Results section, we mostly used twodimensional bounding boxes for graphical convenience. In order to illustrate some properties of higherdimensional error bounding boxes (Appendix 1—figure 3), we compared their behavior against that of hyperspheres and hypercubes. We defined the equivalent hypersphere as $\{\mathit{p}\in {\mathbb{R}}^{M}:{\parallel \mathit{p}\parallel}_{2}\le T\}$ and the equivalent hypercube as $\{\mathit{p}\in {\mathbb{R}}^{M}:{\parallel \mathit{p}\parallel}_{\mathrm{\infty}}\le T\}$, where ${\parallel \mathit{p}\parallel}_{2}=\sqrt{{p}_{1}^{2}+\mathrm{\dots}+{p}_{n}^{2}}$ and ${\parallel \mathit{p}\parallel}_{\mathrm{\infty}}={\mathrm{max}}_{i}{p}_{i}$. In practice, we chose the smallest box size, $T=0.5$ (Appendix 1—figure 3).
For a first comparison, we took the intersection between the border of the $M$dimensional polytope $B$ and a random twodimensional plane containing the centre of the polytope. We computed such intersections numerically by first choosing two random and orthogonal directions $u$ and $v$ in the full space defining the twodimensional plane. Then for each $\theta \in [0,2\pi ]$ , we defined a ray in the twodimensional plane, $w(\rho )=\rho \mathrm{cos}(\theta )u+\rho \mathrm{sin}(\theta )v$, and then plotted
For a second comparison, we found the distribution of angles between neighbouring neurons by first randomly choosing one neuron, and then moving along the surface of the $M$polytope in a random direction, until we found a point that belongs to the face of a different neuron. We then computed the angle between the decoding weights of those two neurons.
We tested whether the results obtained for random decoding vectors hold for more structured decoding vectors as well. For instance, if we want to represent natural visual scenes, we may consider that the receptive fields of simple cells in V1 roughly correspond to the decoding vectors of our neurons (Olshausen and Field, 1996; Barrett et al., 2016). We illustrated a highdimensional bounding box with a set of Gabor patches defined as
where $\stackrel{~}{x}=x\mathrm{cos}\theta +y\mathrm{sin}\theta $ and $\stackrel{~}{y}=x\mathrm{sin}\theta +y\mathrm{cos}\theta $. For our purposes, we randomly chose the Gabor parameters: $\lambda $, the wavelength of the sinusoidal stripe pattern, was sampled uniformly from $\{3,5,10\}$ Hz; $\theta $, the orientation of the stripes, was sampled uniformly in $[0,2\pi ]$, the standard deviation of the Gaussian envelope, was sampled uniformly from $\{1,1.5\}$, the spatial aspect ratio, was sampled uniformly from $\{1,1.5\}$.
Finally we randomly centred the resulting Gabor patch in one of 9 different locations on the 13 × 13 grid. We computed the angle (in the 169dimensional space) between the Gabor patches and found that roughly 80% of the neurons are quasiorthogonal (their angle falls between 85 and 95 degrees) to a given example patch (Appendix 1—figure 3E).
Perturbations
Perturbations to the excitability of a neuron, be it due to changes of the spiking threshold, changes of the reset potential, synaptic weights, etc., can all be formulated as extra currents, ${p}_{i}(t)$, which capture the temporal evolution of the perturbation. Adding a current to the voltage dynamics is equivalent to a transient change in the neuronal thresholds,
Here, $\mathbf{p}(t)$ denotes the vector of current perturbations, and $h\ast \mathbf{p}$ denotes a convolution of the perturbation currents with an exponential kernel, $h(t)$. Note that moving the perturbation onto the threshold does not change the spiking behavior of the neuron. Appendix 1—table 1 includes the range of perturbations used throughout this manuscript.
Voltage noise
Request a detailed protocolWe implement voltage noise as an extra random current on the voltage dynamics. This extra current follows a Wiener process scaled by ${\sigma}_{V}$, which denotes the standard deviation of the noise process with Gaussian increments (see Equation 8). In the absence of recurrence,
so that the leaky integration with time constant $\lambda $ biases the random walk of the thresholds back towards their default values. For stationary inputs, the thresholds therefore follow an OrnsteinUhlenbeck process.
Synaptic perturbations
Request a detailed protocolWe perturb synapses between different neurons ($i\ne j$) by a multiplicative noise term
where ${u}_{i,j}\sim \mathcal{U}(1,1)$. Here, the parameter ${\delta}_{\mathrm{\Omega}}$ is the maximum weight change in percentage of each synapse, which in Figure 6G is referred to as maximal synaptic scaling.
Synaptic delays
Request a detailed protocolWe implement delayed recurrent connections with the same constant delay length $\theta \ge 0$ for all pairs of neurons. Regardless of whether or not lateral excitation and inhibition are delayed in this way, the selfreset of a neuron onto itself remains instantaneous. Equation 3 thus becomes
where ${\delta}_{ik}$ is Kronecker’s delta. We assume that the decoder readout is equally delayed.
Parameter choices
Request a detailed protocolThe spiking networks presented here depend on several parameters:
The number of neurons in the network, $N$.
The number of signals fed into the network, $M$, also called the dimensionality of the signal.
The $M\times N$ matrix of decoding weights, ${D}_{ik}$, where each column ${\mathbf{D}}_{k}$, corresponds to the decoding weights of one neuron.
The inverse time constant of the exponential decay of the readout, $\lambda $.
The threshold (or error tolerances) of the neurons, $T$.
The refractory period, ${\tau}_{\text{ref}}$.
The current noise, ${\sigma}_{V}$.
These parameters fully define both the dynamics and architecture – in terms of feedforward and recurrent connectivity – of the networks, as well as the geometry of the bounding box. We studied networks with various number of neurons $N$ and input dimensionality $M$. The decoding weights of each neuron were drawn from an $M$dimensional standard normal distribution,
and then normalized,
such that each neuronal decoding vector is of length 1. We then did a sweep on the remaining parameters ($\lambda $, $T$, ${\tau}_{\text{ref}}$, ${\sigma}_{V}$), to narrow down the range of parameters that roughly matches key observational constraints, such as low median firing rates ($\sim 5$ Hz), as found in cortex (Hromádka et al., 2008; Wohrer et al., 2013; Figure 5B), and coefficients of variation of interspike intervals close to one for each neuron, corresponding to Poissonlike spike statistics (Figure 5C). Appendix 1—table 1 displays the range of parameters used to simulate baseline and perturbed networks.
Input signal
Request a detailed protocolWe used two different types of inputs throughout our simulations. The results shown in Figure 4 and Appendix 1—figure 5 are for a circular, 2dimensional signal,
with constant amplitude $a$ and constant frequency $\omega $.
For all other simulations shown in figure panels, the input signal was a constant signal with additive noise. More precisely, for each trial, we sampled a single point in input space from an $M$dimensional Gaussian distribution,
The input signal ramps linearly from zero to this point ${\mathit{x}}_{0}$ during the first 400ms. For the rest of the trial, the input to the neurons is set to slowly vary around this chosen input vector. To generate the slow variability, we sampled from an $M$dimensional Gaussian distribution as many times as there were time steps in the rest of the trial; we then twicefiltered the samples with a moving average window of 1s for each dimension of $\mathbf{x}$, and for each dimension of $\mathbf{x}$ and across time, we normalized the individual slow variabilities to not exceed ${\eta}_{x}=0.5$ in magnitude. This procedure was chosen to mimic experimental trialtotrial noise.
Metrics and network benchmarking
To compare the behavior of our networks under baseline conditions to those under the different perturbations, we need reliable measures of both coding accuracy and firing statistics. Below, we describe the measures used in this study.
Distributions of firing rates and coefficients of variation
Request a detailed protocolWe measured the timeaveraged firing rate for a given neuron by dividing the total number of spikes by the total duration of a trial. The coefficient of variation (CV) of a single spike train is computed as the ratio of the standard deviation of the interspike intervals (ISI) to their mean
We recorded the full distributions of both the firing rates and CVs for a given network, pooling across neurons and different trials.
Network performance
Request a detailed protocolWhen our aim was to compare the relative network performance with and without the different perturbations, we opted to use a simple Euclidean distance or L2 norm to measure the average error of each network:
To compute the relative performance, we divided the error of the perturbed network by the error of the equivalent, unperturbed network using the formula
where ${E}_{\text{dead}}$ is the error of a nonfunctional or dead network ($\hat{\mathbf{x}}(t)=0$ or $E}_{\text{dead}}=\u27e8\Vert \mathbf{x}(t)\Vert {\u27e9}_{t$). We included this case to provide a baseline, worstcase scenario.
A key limitation with most error measures is that they scale in various ways with dimensionality. This becomes an issue in Figure 5 as this hinders the comparison of errors across different signal dimensionalities. For this particular case, we chose to measure the coding errors in a dimensionalityindependent way by pooling together the errors in each individual signal component, ${x}_{i}{\widehat{x}}_{i}$. We can then compute the median of this aggregated distribution in order to consistently compare the performance of these networks across different signal dimensionalities.
Excitationinhibition balance
Request a detailed protocolIn order to compute the EI balance of a given neuron $j$, we divided the total synaptic input throughout a given trial into its positive (${C}_{j}^{+}$) and negative (${C}_{j}^{}$) components
The normalized EI difference b_{j} of a neuron $j$ was then computed as
In other words, ${b}_{j}=0$ if a neuron is perfectly balanced, $0<{b}_{j}\le 1$ if a neuron receives more excitation than inhibition and $1\le {b}_{j}<0$ if a neuron receives more inhibition than excitation.
Benchmarking
Request a detailed protocolTo fully compare the behavior of the networks under baseline conditions to those under the different perturbations, we adopted the following benchmarking procedure: each simulated trial with a perturbation is compared to an otherwise identical trial without perturbation. For each trial, we generated a new network with a different random distribution of decoding weights, random input signal, and random voltage noise. These parameters were used for both the perturbed and unperturbed trial. We then applied our $M$dimensional input signal $\mathbf{x}$ as described above, and recorded coding error and spiking statistics for both perturbed and unperturbed trial. This procedure was repeated multiple times (${N}_{\text{trials}}\ge 20$), each repetition resulting in different network connectivity, inputs, and injected current noise, and each pair of trials returning one performance value as defined above.
We choose this benchmarking procedure to sample the space of input signals in an unbiased way. This ensures that network performance is not accidentally dominated by a perfect match, or mismatch, between the fixed decoding weights and a given random input. Particularly bad mismatches may still lead to high decoding errors, but because our error measure considers the median response, these extremes do not bias our benchmarking procedure.
Number of simulations
Request a detailed protocolFigure 1D shows a single trial. Figure 5 shows a total of 29,400 trials. Figure 6E and F show 16,830 pairs of trials, and Appendix 1—figure 4 shows 4996 pairs. Panels Figure 6G and H consist of 840 trials each. Figure 7G show 18,000 pairs of trials, or 200 pairs per data point, and Appendix 1—figure 5 shows 1 perturbed trial per row.
Numerical implementation
Request a detailed protocolWe numerically solve the differential equations (Equation 8) describing the temporal evolution of membrane voltage by the forward EulerMaruyama method. Because of finite simulation time steps, more than one neural threshold may be crossed during the same step, and more than one neuron may thus be eligible to spike. This problem can be circumvented by decreasing the time step, which, however, increases simulation time. To avoid this tradeoff, we essentially slow down time whenever multiple neurons crossed threshold (Appendix 1—algorithm 1). Note that when considering finite delays $\theta $, delayed lateral recurrence arrives only at the end of each time step (Appendix 1—algorithm 2).
We implemented these methods in both MATLAB and Python, and both sets of code can be used interchangeably. Our code for simulation, analysis and figure generation, as well as sample data files can be found at https://github.com/machenslab/boundingbox (copy archived at swh:1:rev:d9ce2cf52e833ecf67dccc796bd8c9dc505f2e00, Calaim, 2022), under a Creative Commons CC BYNCSA 4.0 license.
Iterative adaptation of parameters to avoid pingpong
Request a detailed protocolIn networks with delays, we can avoid pingpong either by increasing box size or by removing a number of strongest excitatory connections. In both cases, we compute the minimum required value offline using an iterative procedure (Appendix 1—algorithm 3). Note that trials must be sufficiently long to avoid falsenegative reports of pingpong.
Movie visualization
Request a detailed protocolAll movies (Videos 1 and 2) were produced in Python, with the exception of the threedimensional visualization of a polytope, for which we used the bensolve toolbox for MATLAB (Löhne and Weißing, 2017).
FORCElearning rate network and perturbations
Request a detailed protocolWe trained a recurrent network of 1000 rate units using FORCE learning (Sussillo and Abbott, 2009) in the absence of any perturbation. The network dynamics are described by the following system of differential equations:
where $\displaystyle \mathbf{r}=\mathrm{tanh}(\mathbf{x})$ corresponds to the firing rates, and where $\tau =10$ms, $g=1.5$, and $\mathbf{J}$ is a sparse random matrix whose elements are zero with probability $1p$. Each nonzero element is drawn independently from a Gaussian distribution with zero mean and variance equal to $(1000p{)}^{1}$. The entries of the matrix $\mathbf{J}}_{z$ are uniformly distributed from –1 to 1.
We then applied one of three perturbations to the fully trained network: neuron death, rate noise, or synaptic perturbation. We emulated neural loss by setting the respective neural activities to zero, i.e. ${x}_{i}=0$. The rate noise perturbation was simulated by injecting white noise within the inputoutput nonlinearity and its magnitude was chosen so that fluctuations on the network activities were of the same order of magnitude as the ones simulated for coordinated spiking networks. Finally, we simulated synaptic perturbations following the same procedure and magnitude as for the coordinated spiking networks, i.e., each element of the recurrent connectivity matrix was changed randomly up to 2.5% of its value.
Appendix 1
Appendix 1—algorithm 1. Numerical implementation of a general network with voltage noise ${\sigma}_{V}$ and refractory period ${\tau}_{\text{ref}}$. 
$K\leftarrow \left\{k\phantom{\rule{thickmathspace}{0ex}}\phantom{\rule{thickmathspace}{0ex}}k\in \mathbb{N}:1\le k\le N\right\}$ //all neurons initialise for $t\mathrm{=}\mathrm{0}$ to $t}_{\mathrm{m}\mathrm{a}\mathrm{x}$ in steps $\mathrm{\Delta}\mathit{}t$ do $R\leftarrow \left\{k\phantom{\rule{thickmathspace}{0ex}}\phantom{\rule{thickmathspace}{0ex}}k\in K:t\underset{{t}^{\mathrm{\prime}}<t}{\text{arg,max}}({s}_{k}({t}^{\mathrm{\prime}})=1)<{\tau}_{\text{ref}}\right\}$ //in refraction $C\leftarrow \left\{k\phantom{\rule{thickmathspace}{0ex}}\phantom{\rule{thickmathspace}{0ex}}k\in K\mathrm{\setminus}R:{V}_{k}(t)>{T}_{k}(t)\right\}$ //spike candidates while $C\mathrm{\ne}\mathrm{\varnothing}$ do $w\leftarrow \underset{k\in C}{\text{arg\hspace{0.17em}max}}\left({V}_{k}(t){T}_{k}(t)\right)$ // furthest above threshold ${s}_{w}(t)\leftarrow 1$ //spike $\mathit{V}(t)\leftarrow \mathit{V}(t){\mathit{D}}^{T}{\mathit{D}}_{w}$ //instant recurrence $R\leftarrow R\cup \{w\}$ // refraction $C\leftarrow \left\{k\phantom{\rule{thickmathspace}{0ex}}\phantom{\rule{thickmathspace}{0ex}}k\in K\mathrm{\setminus}R:{V}_{k}(t)>{T}_{k}(t)\right\}$ // spike candidates end sample $\mathit{\eta}(t)\sim \mathcal{N}(\mathbf{0},{\sigma}_{V}\mathit{I})$ $\mathit{V}(t+\mathrm{\Delta}t)\leftarrow \mathit{V}(t)+\mathrm{\Delta}t\left(\lambda \mathit{V}(t)+\lambda \mathit{D}\mathit{x}(t)\right)+\sqrt{\mathrm{\Delta}t}\mathit{\eta}(t)$ end 
Appendix 1—algorithm 2. Numerical implementation of a general network with finite delays $\theta $, refractory period ${\tau}_{\text{ref}}$, current noise ${\sigma}_{V}$, timevarying synaptic noise $\mathrm{\Delta}\mathbf{\Omega}(t)$ and timevarying optogenetic currents $\mathit{p}(t)$. 
$K\leftarrow \{kk\in \mathbb{N}:1\le k\le N\}$ // all neurons initialise ${V}_{k}(0)\forall k\in K$ $\mathbf{\Omega}={\mathit{D}}^{T}\mathit{D}$ // standard recurrent matrix if $\theta >0$ then ${\mathbf{\Omega}}^{f}=\text{diag}(\mathbf{\Omega})$ // instant selfreset vector ${\mathbf{\Omega}}^{\theta}=\mathbf{\Omega}\text{diag}({\mathbf{\Omega}}^{f})$ // delayed recurrence matrix end for $t\mathrm{=}\mathrm{0}$ to $t}_{\mathrm{m}\mathrm{a}\mathrm{x}$ in steps $\mathrm{\Delta}\mathit{}t$ do sample ${\mathbf{\Omega}}^{\ast}(t)\leftarrow \mathbf{\Omega}+\mathrm{\Delta}\mathbf{\Omega}(t)$ // synaptic noise if $\theta >0$ then ${\mathbf{\Omega}}^{f}=\text{diag}({\mathbf{\Omega}}^{\ast}(t))$ // instant selfreset vector ${\mathbf{\Omega}}^{\theta}={\mathbf{\Omega}}^{\ast}(t)\text{diag}({\mathbf{\Omega}}^{f})$ // delayed recurrence matrix end $R\leftarrow \left\{k\phantom{\rule{thickmathspace}{0ex}}\phantom{\rule{thickmathspace}{0ex}}k\in K:t\underset{{t}^{\mathrm{\prime}}<t}{\text{arg,max}}({s}_{k}({t}^{\mathrm{\prime}})=1)<{\tau}_{\text{ref}}\right\}$ // in refraction $C\leftarrow \left\{k\phantom{\rule{thickmathspace}{0ex}}\phantom{\rule{thickmathspace}{0ex}}k\in K\mathrm{\setminus}R:{V}_{k}(t)>{T}_{k}(t)\right\}$ // spike candidates while $C\mathrm{\ne}\mathrm{\varnothing}$ do $w\leftarrow \underset{k\in C}{\text{arg\hspace{0.17em}max}}\left({V}_{k}(t){T}_{k}(t)\right)$ // furthest above threshold ${s}_{w}(t)\leftarrow 1$ // spike if $\theta >0$ then ${V}_{w}(t)\leftarrow {V}_{w}(t){\mathbf{\Omega}}_{w}^{f}$ // instant selfreset else $\mathit{V}(t)\leftarrow \mathit{V}(t){\mathbf{\Omega}}_{w}^{*}$ // instant recurrence end $R\leftarrow R\cup \{w\}$ // refraction $C\leftarrow \left\{k\phantom{\rule{thickmathspace}{0ex}}\phantom{\rule{thickmathspace}{0ex}}k\in K\mathrm{\setminus}R:{V}_{k}(t)>{T}_{k}(t)\right\}$ // spike candidates end $\mathrm{\Delta}\mathit{V}=\mathrm{\Delta}t\left(\lambda \mathit{V}(t)+\lambda \mathit{D}\mathit{x}(t)\right)$ // dynamics unperturbed network sample $\mathit{\eta}(t)\sim \mathcal{N}(\mathbf{0},{\sigma}_{V}\mathit{I})$ $\mathrm{\Delta}\mathit{V}\leftarrow \mathrm{\Delta}\mathit{V}+\sqrt{\mathrm{\Delta}t}\mathit{\eta}(t)$ // current noise $\mathrm{\Delta}\mathit{V}\leftarrow \mathrm{\Delta}\mathit{V}+\mathrm{\Delta}t\mathit{p}(t)$ // optogenetic currents if $\theta >0$ then $\mathrm{\Delta}\mathit{V}\leftarrow \mathrm{\Delta}\mathit{V}{\mathbf{\Omega}}^{\theta}\mathit{s}(t+\mathrm{\Delta}t\theta )$ // delayed recurrence end $\mathit{V}(t+\mathrm{\Delta}t)\leftarrow \mathit{V}(t)+\mathrm{\Delta}\mathit{V}$ end 
Appendix 1—algorithm 3. Numerical search for the "safe width" of a bounding box, avoiding pingpong. Typical parameters are ${T}_{\text{min}}=0.55$, $\alpha =1.5$, $\beta =0.95$, $\gamma =0.1$, $\u03f5=0.05\cdot 2\theta $, $N=100$. In each trial, all neurons $j$ have the same threshold ${T}_{j}$, and the box is thus widened or narrowed symmetrically. 
initialise $T\leftarrow {T}_{\text{min}}>0$ // current box width initialise ${T}^{\ast}\leftarrow 0$ // best box width so far initialise $k\leftarrow 0$ // trial counter while $k<K$ do $k\leftarrow k+1$ simulate network with $N$ neurons and box width $T$ for $1<j\le N$ do ${\mathrm{\Theta}}_{j}\leftarrow \{t{s}_{j}(t)=1\}$ // spike times $S}_{j}\leftarrow \{t{t}^{\mathrm{\prime}}\phantom{\rule{thickmathspace}{0ex}}\phantom{\rule{thickmathspace}{0ex}}t,{t}^{\mathrm{\prime}}\in {\mathrm{\Theta}}_{j}\phantom{\rule{thickmathspace}{0ex}}\wedge \phantom{\rule{thickmathspace}{0ex}}t=\underset{x}{\mathrm{a}\mathrm{r}\mathrm{g}\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{i}\mathrm{n}}\left(x>{t}^{\mathrm{\prime}}\right)\$ // intervals end $S\leftarrow {\bigcup}_{j=1}^{N}{S}_{j}$ // pool interspike intervals $A\leftarrow \left\{a\in S\phantom{\rule{thickmathspace}{0ex}}\phantom{\rule{thickmathspace}{0ex}}2\theta \u03f5<a<2\theta +\u03f5\right\}$ // SISIs near doubledelay $P\leftarrow \frac{\leftA\right}{\leftS\right}>\gamma$ // Boolean: pingpong present? if $P$ then if ${w}^{\ast}>0$ then $w\leftarrow {T}^{\ast}$ // use previous estimate... $k\leftarrow K$ //...and quit else $T\leftarrow \alpha T$ // increase box size $k\leftarrow 0$ // restart trial counter end else if $k\mathrm{=}N$ then ${T}^{\ast}\leftarrow w$ // update best estimate $T\leftarrow \beta T$ // slightly decrease box size $k\leftarrow 0$ // restart trial counter end end 
Data availability
The current manuscript is a computational study, so no data have been generated for this manuscript. Modelling code is uploaded on https://github.com/machenslab/boundingbox, (copy archived at swh:1:rev:d9ce2cf52e833ecf67dccc796bd8c9dc505f2e00).
References

Recurrent neural networks as versatile tools of neuroscience researchCurrent Opinion in Neurobiology 46:1–6.https://doi.org/10.1016/j.conb.2017.06.003

ConferenceFiring rate predictions in optimal balanced networksAdvances in Neural Information Processing Systems 26.

Analyzing biological and artificial neural networks: challenges with opportunities for synergy?Current Opinion in Neurobiology 55:55–64.https://doi.org/10.1016/j.conb.2019.01.007

ConferenceEvasion attacks against machine learning at test timeJoint European conference on machine learning and knowledge discovery in databases. pp. 387–402.https://doi.org/10.1007/9783642387098

A neuromorph’s prospectusComputing in Science & Engineering 19:14–28.https://doi.org/10.1109/MCSE.2017.33

Predictive coding of dynamical variables in balanced spiking networksPLOS Computational Biology 9:e1003258.https://doi.org/10.1371/journal.pcbi.1003258

ConferenceLearning optimal spikebased representationsAdvances in Neural Information Processing Systems 25.

Learning to represent signals spike by spikePLOS Computational Biology 16:e1007692.https://doi.org/10.1371/journal.pcbi.1007692

Reverse engineering of biological complexityScience 295:1664–1669.https://doi.org/10.1126/science.1069981

Dimensionality reduction for largescale neural recordingsNature Neuroscience 17:1500–1509.https://doi.org/10.1038/nn.3776

BookNeural Engineering: Computation, Representation, and Dynamics in Neurobiological SystemsMIT Press.

A unified approach to building and controlling spiking attractor networksNeural Computation 17:1276–1314.https://doi.org/10.1162/0899766053630332

Pervasive robustness in biological systemsNature Reviews. Genetics 16:483–496.https://doi.org/10.1038/nrg3949

BookNeuronal Dynamics: From Single Neurons to Networks and Models of CognitionCambridge University Press.https://doi.org/10.1017/CBO9781107447615

Decoding and encoding (de)mixed population responsesCurrent Opinion in Neurobiology 58:112–121.https://doi.org/10.1016/j.conb.2019.09.004

Balanced cortical microcircuitry for maintaining information in working memoryNature Neuroscience 16:1306–1314.https://doi.org/10.1038/nn.3492

The vector linear program solver Bensolve – notes on theoretical backgroundEuropean Journal of Operational Research 260:807–813.https://doi.org/10.1016/j.ejor.2016.02.039

ConferenceUnderstanding spiking networks through convex optimizationAdvances in Neural Information Processing Systems 33.

Poisson balanced spiking networksPLOS Computational Biology 16:e1008261.https://doi.org/10.1371/journal.pcbi.1008261

ConferenceSpatiotemporal representations of uncertainty in spiking neural networksAdvances in Neural Information Processing Systems 27.

Towards the neural population doctrineCurrent Opinion in Neurobiology 55:103–111.https://doi.org/10.1016/j.conb.2019.02.002

Learning Universal Computations with SpikesPLOS Computational Biology 12:e1004895.https://doi.org/10.1371/journal.pcbi.1004895

Recoding a cocaineplace memory engram to a neutral engram in the hippocampusNature Neuroscience 19:564–567.https://doi.org/10.1038/nn.4250

Homeostatic synaptic plasticity: local and global mechanisms for stabilizing neuronal functionCold Spring Harbor Perspectives in Biology 4:a005736.https://doi.org/10.1101/cshperspect.a005736

Neural network dynamicsAnnual Review of Neuroscience 28:357–376.https://doi.org/10.1146/annurev.neuro.28.061604.135637

Computation Through Neural Population DynamicsAnnual Review of Neuroscience 43:249–275.https://doi.org/10.1146/annurevneuro092619094115

Populationwide distributions of neural activity during perceptual decisionmakingProgress in Neurobiology 103:156–193.https://doi.org/10.1016/j.pneurobio.2012.09.004

The promise and perils of causal circuit manipulationsCurrent Opinion in Neurobiology 49:84–94.https://doi.org/10.1016/j.conb.2018.01.004

Representation of spatial orientation by the intrinsic dynamics of the headdirection cell ensemble: a theoryThe Journal of Neuroscience 16:2112–2126.
Decision letter

Markus MeisterReviewing Editor; California Institute of Technology, United States

Michael J FrankSenior Editor; Brown University, United States

Brent DoironReviewer; University of Pittsburgh, United States

Markus MeisterReviewer; California Institute of Technology, United States
Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Decision letter after peer review:
Thank you for submitting your article "Robustness in spiking networks: a geometric perspective" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Markus Meister as Reviewing Editor and Reviewer #3, and the evaluation has been overseen by Michael Frank as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Brent Doiron (Reviewer #2).
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential revisions:
Domain of applicability:
1. The authors analyze a very specific network, with carefully tuned feedback designed for a specific cost function (Equation 3). Do such tuned networks exist in Nature? The reader would benefit from knowing what the domain of applicability really is.
2. How would the framework generalize to settings other than an autoencoder, e.g., for the dynamical systems approach in Boerlin et al.?
Consistency of the bounding box picture:
3. A central assumption of the framework is a coherence between the readout from the network and the dynamics and connectivity of the network. This allows the bounding box to follow the input vector around, so one can easily assess decoding errors. However, it seems that some of the perturbations break the assumptions that lead to the bounding box. For example, exciting a single neuron decreases its threshold and thereby shrinks the bounding box (Figure 3D). As the bounding box limits the error in the readout, this should mean that the error decreases. However, the opposite is the case (Figure 4A). Presumably the discrepancy arises from changing the threshold without changing the readout at the same time, so the bounding box loses its meaning. By contrast the loss of a neuron, or its inhibition, leaves the remaining network adjusted properly, so the bounding box retains its meaning. The authors should address this question, and whether the bounding box loses its meaning under some perturbations and not others.
4. A related question applies to delays. Where exactly are the delays? In the readout and then the network is optimised for that? Or is the network the optimal one assuming no delays and then delays are added in the recurrent connections and/or the readout? Clarifying that will make it easier for the reader to follow the argument.
5. In the section on perturbations of the reset potential: Is there again an asymmetry, e.g., do changes up or down in the optimal reset potential have very different consequences?
Robustness claims:
6. There are frequent comparisons between networks with few neurons and networks with many neurons. The idea that networks with more neurons are more robust seems kind of obvious. But the authors analyze a very specific network tuned for coordinated redundancy (Equation 3). How does that network compare to one with the same number of neurons but no coordinated redundancy? This would seem a more interesting comparison than varying the number of neurons.
7. Intuitively one might expect that the performance is fragile to changes in the recurrent coupling, since the decoding vector is tied to the feedforward and recurrent weights (Equation 10). Because of this the claimed robustness to perturbations in synaptic weight seems surprising. Is the 'synapse scaling factor' in Figure 6G,H the parameter \δ_{\Omega} in Equation 17? The way it is defined it is the maximal change, however most synapses are perturbed to a much smaller extent. Also, for \δ approaching 0.2 the performance only drops by 20% (Figure 6G). However, for such large \δ's the firing rates in the network seem unreasonably high (Figure 6H), suggesting it no longer performs efficiently. Overall it is hard to judge from the current presentation how robust the network is to synaptic perturbation. Perhaps the firing rates could be bounded, or the performance penalized for using very high rates?
Relation to plasticity:
8. Prior work (e.g. Bourdoukan 2012) proposed a Hebbian learning process that could tune up such a network automatically. If this plasticity is "always on" how would it respond to the various perturbations being considered? How might that contribute to robustness on a longer time scale? How does the size/shape of the bounding box depend on the parameters of the learning rule? Can one get looser and tighter boxes, as needed here for certain types of robustness?
https://doi.org/10.7554/eLife.73276.sa1Author response
Essential revisions:
Domain of applicability:
1. The authors analyze a very specific network, with carefully tuned feedback designed for a specific cost function (Equation 3). Do such tuned networks exist in Nature? The reader would benefit from knowing what the domain of applicability really is.
The reviewers are correct in that our networks require carefully tuned inhibitory feedback. We first notice that this is not a new proposal. For instance, in previous work we have shown that the loss functions we minimize are identical to classical sparse coding networks (Olshausen and Field, 1996). Indeed, the type of connectivity we find in our networks is exactly the same as that in sparse coding rate networks (compare e.g., Figure 10.6 in Dayan Abbott with Equation (8) in our paper). So we could say a potential applicability are sensory systems, such as primary visual cortex. Do such tuned networks really exist in nature? That is a more difficult question to answer. There is currently an ongoing debate on whether inhibitory tuning is precise or broad, and so we would say that the jury is still out (see e.g., Najafi et al., (2020) Neuron 105:16579, for an argument that inhibitory tuning is precise or selective). We have sought to clarify these points in the discussion where we now write:
"In many respects, the bounding box is a "toy model" (borrowing the term from physics, in the sense of a deliberately abstract model), which we see mainly as a tool to conceptualize and highlight generic circuit mechanism, rather than an attempt to model any specific system. Nonetheless, it is worthwhile to point out that the bounding box is also a spiking version of classical sparse coding models of V1 [44]. Indeed, previous work has demonstrated that these spiking networks can explain various perturbation experiments in V1 [21]. So, besides shedding light on the robustness of coordinated spike codes, the bounding box can also be seen as a simple model of a sensory system."
2. How would the framework generalize to settings other than an autoencoder, e.g., for the dynamical systems approach in Boerlin et al.?
The framework fully generalizes. In the case of an autoencoder, the positions of bounding box faces are fully defined by the input. For the more general dynamical systems approach of Boerlin et al., there still is a bounding box that arises from the fast recurrent connectivity, and the position of its faces evolves in signal space according to the chosen dynamical system. See Supplementary Figure 6 for a graphical explanation of this comparison. We also note that the simulation in Figure 1D of the main paper is actually a dynamical system à la Boerlin et al., and not an autoencoder.
To address these concerns, we added Supplementary Figure 6, and we introduced a new section at the end of the Results to clarify that our framework generalizes to Boerlin et al., Specifically, we wrote:
"The bounding box provides a useful tool even if we endow the networks with a set of slower connections to perform linear or nonlinear computations [17, 42, 43]. Indeed, the simulation in Figure 1D used these slower connections to generate oscillatory dynamics (see Methods, section ’Generalisation of the bounding box IV’, for mathematical details). This extension to networks that generate persistent activity or dynamical patterns works because the mechanisms underlying the encoding of the signals into spike trains are decoupled from the mechanisms that generate the dynamics of the signals (or readouts). Accordingly, the extra currents generated by the slow recurrent connections can be seen as a perturbation of the bounding box thresholds. This perturbation shifts the bounding box in the space of readouts as illustrated in Supplementary Figure 6."
Consistency of the bounding box picture:
3. A central assumption of the framework is a coherence between the readout from the network and the dynamics and connectivity of the network. This allows the bounding box to follow the input vector around, so one can easily assess decoding errors. However, it seems that some of the perturbations break the assumptions that lead to the bounding box. For example, exciting a single neuron decreases its threshold and thereby shrinks the bounding box (Figure 3D). As the bounding box limits the error in the readout, this should mean that the error decreases. However, the opposite is the case (Figure 4A). Presumably the discrepancy arises from changing the threshold without changing the readout at the same time, so the bounding box loses its meaning. By contrast the loss of a neuron, or its inhibition, leaves the remaining network adjusted properly, so the bounding box retains its meaning. The authors should address this question, and whether the bounding box loses its meaning under some perturbations and not others.
The reviewers are correct that the readout is a central pillar of our analysis. There is one important subtlety, though. For wide bounding boxes, i.e., for bounding boxes in which the readout jump caused by a spike does not transverse the whole box, we need to correct the readout for a bias (see Author response image 1, replotted from Supplementary Figure 1A). This biascorrection is done by rescaling the readout with a scalar. The bias correction achieves that the average rescaled readout will match the signal.
The bounding box only limits the coding error of the actual readout, but not the rescaled, bias corrected readout. In other words, the bounding box only has meaning for the original readout, not the corrected readout. When we shrink the bounding box from one side, we decrease the maximum error of the original readout. The biascorrected, rescaled readout, however, will now, on average, no longer match the signal. That is the effect we plot in Figure 3D.
To make these issues clearer, we now clearly specify throughout the text whether we use the readout or the corrected readout. When explaining readout and corrected readout, we explain the limits of the bounding box picture: "This bias can be largely eliminated by rescaling the readouts with a constant factor. We will sometimes use this corrected readout (see Methods, ’Readout biases’), but note that the corrected readout is not confined to stay within the bounding box."
We also added the following sentence when explaining changes in the excitability of neurons:
"At first sight, changing the box size increases or decreases the maximum error of the readout. More subtly, however, it also introduces a bias in the corrected average readout (Figure 3C, arrows)."
4. A related question applies to delays. Where exactly are the delays? In the readout and then the network is optimised for that? Or is the network the optimal one assuming no delays and then delays are added in the recurrent connections and/or the readout? Clarifying that will make it easier for the reader to follow the argument.
Crucially, the delays are in the recurrent connections. The reviewers’ second assertion is correct: Network connectivity is optimized for the delayfree case, and the delays are then added to the recurrent connections. In addition to these recurrent delays, we have also assumed a delay of identical length in the readout. This latter delay is not crucial, but facilitates the geometric visualisation, as the arrival of a delayed recurrent spike and the update of the readout thus happen at the exact same time and can thus be shown in the same figure panel.
We now clarify our choices in the Results section, where we have inserted the following statement:
"Below, we study the impact of these delays, which apply directly to recurrent excitation and inhibition. We also apply the same delays to the network readout for mathematical convenience, but those do not affect network dynamics (see Methods)."
5. In the section on perturbations of the reset potential: Is there again an asymmetry, e.g., do changes up or down in the optimal reset potential have very different consequences?
The reviewers are correct. A change in reset potential leads to a transient change in the corresponding face of the bounding box, and inward changes (a transient shrinkage of the box) can have very different effects than an outward change. Indeed, a reset that is too small can cause a successive shrinking of the box if the respective neuron fires repeatedly.
We now included the following statement in the Results section:
“We note that positive and negative changes to the default reset potential will lead to asymmetric effects on robustness like those observed for excitatory and inhibitory perturbations. Specifically, if the resets become too small, and if the leak is insufficiently fast, then successive spiking of a single neuron will draw its threshold inwards, thereby leading to a collapse of the bounding box."
Robustness claims:
6. There are frequent comparisons between networks with few neurons and networks with many neurons. The idea that networks with more neurons are more robust seems kind of obvious. But the authors analyze a very specific network tuned for coordinated redundancy (Equation 3). How does that network compare to one with the same number of neurons but no coordinated redundancy? This would seem a more interesting comparison than varying the number of neurons.
We agree with the reviewer that a more systematic comparison with other types of networks would be preferable, but we have found that there are two key problems. First, there is no agreed, default network model that we could compare our framework against. Second, almost all standard network models are not built to be robust, so that systematic comparisons are not very illuminating.
We have addressed these issues in the following way:
First, we note that we include a comparison with networks without coordinated redundancy. We simply compare our networks against a set of unconnected neurons (Section ’passive redundancy’). A set of unconnected neurons is a simple example, as any perturbation effect scales linearly with the number of neurons perturbed. Here, increasing the number of neurons will not guard against the elimination of a certain fraction of the network (say 25%). In the bounding box picture, however, increasing the number of neurons (or the redundancy) will indeed increase robustness against killing one fourth of the network, as smaller networks are not necessarily robust. (We note that in both cases, we assume the same downstream readout before and after the perturbation.) This illustrates that more redundancy is not trivially better, at least not if the perturbation size scales with the network size. To highlight this better, we added the following sentence in the section "Scaling up": "This contrasts with networks of independent neurons in which performance will scale linearly with any change in redundancy for a fixed readout."
Second, we note that we illustrated the response of a trained network model without coordinated redundancy to various perturbations. This network is not robust to our cumulative perturbations, yet our identically sized network with coordinated redundancy is robust (Figure 1C vs 1D). While we only include one example simulation, we emphasize that this simulation is representative (which any reader can see by simulating our code).
Third, we would like to point out that we also show that more redundancy is not always better. Notably, our work suggests that including more neurons makes the network more sensitive to certain kinds of perturbations (e.g., synaptic noise in Figure 6GH, or delays in Figure 7G).
7. Intuitively one might expect that the performance is fragile to changes in the recurrent coupling, since the decoding vector is tied to the feedforward and recurrent weights (Equation 10). Because of this the claimed robustness to perturbations in synaptic weight seems surprising. Is the 'synapse scaling factor' in Figure 6G,H the parameter \δ_{\Omega} in Equation 17? The way it is defined it is the maximal change, however most synapses are perturbed to a much smaller extent. Also, for \δ approaching 0.2 the performance only drops by 20% (Figure 6G). However, for such large \δ's the firing rates in the network seem unreasonably high (Figure 6H), suggesting it no longer performs efficiently. Overall it is hard to judge from the current presentation how robust the network is to synaptic perturbation. Perhaps the firing rates could be bounded, or the performance penalized for using very high rates?
The reviewer is correct that the synapse scaling factor in Figure 6G,H is the parameter δ_{Ω} in Equation 17, and that this corresponds to the maximal possible scaling rather than the actual scaling of the synaptic strengths.
To explain the robustness of the network against changes in the synapses, we have to first note the asymmetry of perturbations: a synapse that is smaller than its ideal value will lead to a temporary shift of the postsynaptic neuron’s threshold away from the bounding box when the presynaptic neuron spikes (Figure 6D). Just as with inhibitory perturbations, this perturbations is harmless as it does not really change the overall shape of the box in redundant systems. A synapse that is larger than its ideal value will lead to a temporary shift into the bounding box which can be harmful, depending on how long it lasts (or how fast the synaptic input decays).
We have now sought to clarify this section by writing:
"We again note an asymmetry: a synapse with decreased strength leads to an outward move of the postsynaptic neuron’s threshold, which is generally harmless. Random synaptic failures, which cause temporary decreases in synaptic strength, do therefore not influence the bounding box functionality. However, a synapse with increased strength leads to an inward move, which could be a temporarily harmful perturbation."
We also agree with the reviewer that synaptic scaling can make these networks less inefficient in terms of number of spikes (Figure 6H). We now write:
"Overall, we find that more redundant networks (with consequently more synapses) are typically more vulnerable to these perturbations, and that synaptic scaling can lead to highly inefficient networks in terms of spike rate, regardless of the network redundancy."
Relation to plasticity:
8. Prior work (e.g. Bourdoukan 2012) proposed a Hebbian learning process that could tune up such a network automatically. If this plasticity is "always on" how would it respond to the various perturbations being considered? How might that contribute to robustness on a longer time scale? How does the size/shape of the bounding box depend on the parameters of the learning rule? Can one get looser and tighter boxes, as needed here for certain types of robustness?
The learning rules we have previously developed (Bourdoukan et al., 2012, Neurips; Brendel et al., 2020, PLOS CB), will learn an optimally arranged bounding box with a predefined width. We have not systematically studied how perturbations would affect the learning—a very interesting question, but we believe beyond the confines of the current study.
Still, a few answers we already know. Specifically, perturbations that lower the excitability of neurons are unlikely to affect learning. For instance, when eliminating neurons, the target connectivity of the learning rules does not change. Moreover, the learning rules we developed will still push synapses in the remaining neurons towards this target connectivity. In machinelearning language, the networks are fully capable of dealing with the dropout of neurons, and dropout may even help in speeding up learning (although we have not full studied this effect). More generally, inhibiting neurons does not affect learning because of the network’s compensatory properties. Perturbations that excite neurons (e.g. synaptic noise) are different. We would speculate that they are likely to strongly perturb the learning process.
https://doi.org/10.7554/eLife.73276.sa2Article and author information
Author details
Funding
Fundação para a Ciência e a Tecnologia (FCTPTDC/BIAOUT/32077/2017IC&DTLISBOA010145FEDER)
 Christian K Machens
Simons Foundation (543009)
 Christian K Machens
Fundação para a Ciência e a Tecnologia (SFRH / BD / 52217 / 2013)
 Nuno Calaim
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Alfonso Renart for helpful discussions and comments on the manuscript. This work was supported by the Fundação para a Ciência e a Tecnologia (project FCTPTDC/BIAOUT/32077/2017IC&DTLISBOA010145FEDER) and by the Simons Foundation (Simons Collaboration on the Global Brain #543009).
Senior Editor
 Michael J Frank, Brown University, United States
Reviewing Editor
 Markus Meister, California Institute of Technology, United States
Reviewers
 Brent Doiron, University of Pittsburgh, United States
 Markus Meister, California Institute of Technology, United States
Publication history
 Preprint posted: June 15, 2020 (view preprint)
 Received: August 23, 2021
 Accepted: May 22, 2022
 Accepted Manuscript published: May 30, 2022 (version 1)
 Version of Record published: July 22, 2022 (version 2)
Copyright
© 2022, Calaim, Dehmelt, Gonçalves et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,506
 Page views

 634
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
Improving muscle function has great potential to improve the quality of life. To identify novel regulators of skeletal muscle metabolism and function, we performed a proteomic analysis of gastrocnemius muscle from 73 genetically distinct inbred mouse strains, and integrated the data with previously acquired genomics and >300 molecular/phenotypic traits via quantitative trait loci mapping and correlation network analysis. These data identified thousands of associations between protein abundance and phenotypes and can be accessed online (https://muscle.coffeeprot.com/) to identify regulators of muscle function. We used this resource to prioritize targets for a functional genomic screen in human bioengineered skeletal muscle. This identified several negative regulators of muscle function including UFC1, an E2 ligase for protein UFMylation. We show UFMylation is upregulated in a mouse model of amyotrophic lateral sclerosis, a disease that involves muscle atrophy. Furthermore, in vivo knockdown of UFMylation increased contraction force, implicating its role as a negative regulator of skeletal muscle function.

 Computational and Systems Biology
Intracellular states probed by gene expression profiles and metabolic activities are intrinsically noisy, causing phenotypic variations among cellular lineages. Understanding the adaptive and evolutionary roles of such variations requires clarifying their linkage to population growth rates. Extending a cell lineage statistics framework, here we show that a population’s growth rate can be expanded by the cumulants of a fitness landscape that characterize how fitness distributes in a population. The expansion enables quantifying the contribution of each cumulant, such as variance and skewness, to population growth. We introduce a function that contains all the essential information of cell lineage statistics, including mean lineage fitness and selection strength. We reveal a relation between fitness heterogeneity and population growth rate response to perturbation. We apply the framework to experimental cell lineage data from bacteria to mammalian cells, revealing distinct levels of growth rate gain from fitness heterogeneity across environments and organisms. Furthermore, third or higher order cumulants’ contributions are negligible under constant growth conditions but could be significant in regrowing processes from growtharrested conditions. We identify cellular populations in which selection leads to an increase of fitness variance among lineages in retrospective statistics compared to chronological statistics. The framework assumes no particular growth models or environmental conditions, and is thus applicable to various biological phenomena for which phenotypic heterogeneity and cellular proliferation are important.