Statistical mechanics of phase-space partitioning in large-scale spiking neuron circuits

TL;DR

This paper develops a statistical theory to understand how the phase space of large-scale spiking neural circuits is partitioned into attractor basins, revealing how connectivity and dynamics influence neural computation.

Contribution

It introduces a novel analytical framework linking spike-time collision events to phase space partitioning in spiking neural networks, advancing understanding of neural circuit dynamics.

Findings

Basin boundaries are pre-images of spike-time collision events.

Basin diameter increases with inhibitory coupling strength.

Basin size shrinks with higher spike event rates.

Abstract

Synaptic interactions structure the phase space of the dynamics of neural circuits and constrain neural computation. Understanding how requires methods that handle those discrete interactions, yet few exist. Recently, it was discovered that even random networks exhibit dynamics that partitions the phase space into numerous attractor basins. Here we utilize this phenomenon to develop theory for the geometry of phase space partitioning in spiking neural circuits. We find basin boundaries structuring the phase space are pre-images of spike-time collision events. Formulating a statistical theory of spike-time collision events, we derive expressions for the rate of divergence of neighboring basins and for their size distribution. This theory reveals that the typical basin diameter grows with inhibitory coupling strength and shrinks with the rate of spike events. Our study provides an analytical and generalizable approach for dissecting how connectivity, coupling strength, single neuron dynamics and population activity shape the phase space geometry of spiking circuits.

Figures

Equations

$$\tau {\dot{V}}_n(t)=-V_n(t)+I_n\left(t\right),$$

$$I_n\left(t\right)=I_{\mathrm{Ext}}+\tau J\sum _s{A}_{n{n}_s}\delta (t-t_s),$$

$$\bar{\nu }=\frac{I_0}{J_0\tau }+\mathcal{O}\left(\frac{1}{\sqrt{K}}\right).$$

$${\varphi }_n(t)=\frac{\tau }{T_{\mathrm{free}}}\ln \left[\frac{I_{\mathrm{Ext}}-V_R}{I_{\mathrm{Ext}}-V_n(t)}\right],$$

$${\dot{\varphi }}_n\left(t\right)=T_{free}^{-1},$$

$$T_{\mathrm{free}}$$

$$Z\left(\varphi \right)$$

$$\text{with}d:$$

$${\lambda }_{\mathrm{inh}}\approx -Kd.$$

$${ϵ}^{\ast }({\vec{\varphi }}_0,\vec{\xi }):=\sup \left\{ϵ|\underset{t\to \infty }{\lim }{D}_t\left(ϵ\right)=0\right\}.$$

$$D_t\left(ϵ\right):=\frac{1}{N}\sum _{n=1}^N\left|{\varphi }_{n,t}-{\varphi }_{n,t}\left(ϵ\right)\right|,$$

$$\Delta t_{\mathrm{jump}}$$

$${\lambda }_{\mathrm{sus}}=p/\overline{\Delta t}=K\bar{\nu },$$

$$\hat{S}\left(ϵ\right)={\left[1_{\mathrm{FT}}\left(ϵ\right)\right]}_{\rho \left({\vec{\varphi }}_0,\vec{\xi }\right)},$$

$$\Delta t_s\left(ϵ\right)=t_s(ϵ)-t_{s-1}(ϵ)=\Delta t_s\left(0\right)+\delta t_s\left(ϵ\right)-\delta t_{s-1}\left(ϵ\right),s\ge 2.$$

$$\delta t_s\left(ϵ\right)\approx -C{a}_sϵ,$$

$$\Delta t_s\left(ϵ\right)\approx \Delta t_s-C\Delta a_sϵ,$$

$$S\left(ϵ\right)={\left[\hat{S}\left(ϵ\right)\right]}_{P_A\left(\left(A_{mn}\right)\right)}.$$

$${\rho }_T=\rho (\left\{\Delta a_s\right\},\left\{\Delta t_s\right\},M,{\vec{\varphi }}_0|\vec{\xi },\left(A_{mn}\right))\rho (\vec{\xi })P_A\left(\left(A_{mn}\right)\right),$$

$${\rho }_T\sim P_{A_{mn}}\left(A_{mn}\right)P_T(M)\prod _{s=2}^M{\rho }_t\left(\Delta t\right)2\Theta (\Delta a_s){\rho }_a\left(\Delta a_s\right),$$

$$S\left(ϵ\right)$$

$$S_s(ϵ)$$

$$S\left(ϵ\right)\approx$$

$$S\left(ϵ\right)$$

$$\overline{{ϵ}^{\ast }}$$

$$=\sqrt{N}\overline{\Delta t}\frac{\left|J\right|/I_{\mathrm{Ext}}}{\tau /I_{\mathrm{Ext}}}=\sqrt{N}\overline{\Delta t}\frac{\left|J\right|}{\tau }=\frac{\sqrt{N}}{\tau }\frac{\left|-J_0/\sqrt{K}\right|}{N\bar{\nu }}$$

$$\overline{{ϵ}^{\ast }}$$

Full text

{affiliations}

Max Planck Institute for Dynamics and Self-Organization, Göttingen, Germany

Bernstein Center for Computational Neuroscience, Göttingen, Germany

Laboratoire de Physique Théorique, ENS-PSL Research University, Paris, France

Kavli Institute for Theoretical Physics, University of California Santa Barbara, Santa Barbara, USA

{frontpage}

Maximilian Puelma Touzel

Laboratoire de Physique Théorique

Ecole Normale Supérieure

24 rue Lhomond

75231 Paris Cedex 05

Paris, France

Ph. +49 551 51 76 420

e-mail:  puelma@lpt.ens.fr

Key words: neuronal circuits, dynamics of networks, disordered systems, basins of attraction, high dimensional systems, pulse-coupled systems, sequence generation

{introduction}

Computing devices, whether natural or artificial, perform their function by finely orchestrated state changes of internal dynamical variables. In nervous systems these dynamical variables are physico-chemical states of nerve cells and synapses that connect them into complex networks called neural circuits. The causal dependencies arising from the synaptic interactions between cells greatly extend the space of functions computable by the circuit, beyond that of single neurons.

Mathematical models of neural circuits have been formulated in two fundamentally distinct ways1. Most synaptic interactions in the brain are driven by sparsely-fired nerve impulses, called spikes, each lasting only a millisecond. In spiking neural network models this fundamental granularity of neuronal interactions is explicitly represented: all interactions depend on a discrete set of spike event times. Alternatively, continuous variable models for neural circuit dynamics are formulated by assuming that a frequency of nerve impulse generation, the firing rate, represents the information-encoding variable causally relevant for neural circuit computation. Firing rate models have been commonly used to model neural circuits 2, theoretically study their dynamics 3; 4 and learning 5; 6; 7; 8, and are the basis of spectacular advances in artificial computing systems 9. Statistical physics has played a role in this development, e.g. in clarifying the disordered phase space organization 10.

From a dynamical systems perspective, attractor states and their basins of attraction play a fundamental role in theories of neural computation. While analogous in some cases 11; 12, however, rate models are not equivalent to temporally coarse-grained versions of spiking neural networks, even if they are closely matched in structure 13. Moreover, low firing rates (not much more than 1 Hz) in the cerebral cortex 14 make it hard to imagine how continuous rate variables associated to single neurons could provide a causally accurate description on behavioral time scales (hundreds of milliseconds). Developing theory for spiking networks may well require a dedicated approach. The absence of relevant averages and even a tractable ensemble of spiking trajectories, however, has thus far limited statistical approaches. Methods to design them 15 or to theoretically dissect the associated phase space organization are only starting to emerge.

Recently it has been discovered that, with dominant inhibition, even randomly wired networks partition their phase space into a complex set of basins of attraction, termed flux tubes 16; 17. Here we utilize this setting to develop a statistical theory of phase space partitioning in spiking neural circuits. We first present a simulation study of flux tubes, uncovering their shape and revealing it is structured by a spike time collision event. Formulating these events, we then derive the conditions for and rate of the mutual divergence of neighboring tubes. Our main calculation is the derivation of the distribution of flux tube sizes, which we obtain from statistics of these events by leveraging the random connectivity to average over the disorder. Our analytical approach provides a transparent method to determine how coupling strength, connectivity, single neuron dynamics and population activity combine to shape the phase space geometry of spiking neural circuits.

Methods

We study a tractable instance of the inhibition-dominated regime of neural circuits. $N$ neurons are connected by an Erdős-Rényi graph with adjacency matrix $A=\left(A_{mn}\right)$. $A_{mn}=1$ denotes a connection from neuron $n$ to $m$, realized with probability, $p=K/N$. The neurons’ membrane potentials, $V_{n}\in\left(-\infty,V_{T}\right]$, are governed by Leaky Integrate-and-Fire (LIF) dynamics,

[TABLE]

for $n\in\left\{1,\dots,N\right\}$. Here, $\tau$ is the membrane time constant and $I_{n}\left(t\right)$ the synaptic current received by neuron $n$; when $V_{n}$ reaches a threshold, $V_{T}=0$, neuron $n$ ‘spikes’, and $V_{n}$ is reset to $V_{R}=-1$. At the spike time, $t_{s}$, the spiking neuron, $n_{s}$, delivers a current pulse of strength $J$ to its $\mathcal{O}(K)$ postsynaptic neurons, $\left\{m|A_{mn_{s}}=1\right\}$, ($s$ indexes the spikes in the observation window). The total synaptic current is

[TABLE]

where $I_{\mathrm{Ext}}>0$ is a constant external current and $J<0$ is the recurrent coupling strength. An $\mathcal{O}\left(1/\sqrt{K}\right)$-scaling of $J$ is chosen to maintain finite current fluctuations at large $K$ and implies that the external drive is balanced by the recurrent input. As a consequence, firing in this network is robustly asynchronous and irregular 18; 19; 20; 21. Setting $I_{\mathrm{Ext}}=\sqrt{K}I_{0}$, with $I_{0}>0$, and $J=-J_{0}/\sqrt{K}$ with $J_{0}>0$, the corresponding stationary mean-field equation for the population-averaged firing rate, $\bar{\nu}$, is 17

[TABLE]

It is convenient to map the voltage dynamics to a pseudophase representation 22; 17, $\vec{\phi}\left(t\right)$, with

[TABLE]

where $T_{free}$ is the oscillation period of a neuron driven only by $I_{\mathrm{Ext}}$. $\phi_{n}\left(t\right)$ evolves linearly in time,

[TABLE]

between spike events, i.e. $t\notin\left\{t_{s}\right\}$, and undergoes shifts given by the phase response curve, $Z(\phi)$, across input spike times where $\phi$ is the state at spike reception. In the large-$K$ limit, $T_{\mathrm{free}}$ and $Z\left(\phi\right)$ simplify to

[TABLE]

respectively (see Supplemental Methods for details). The differential phase response, $\tfrac{\mathrm{d}}{\mathrm{d}\phi}Z=-d$, is essential for the strongly dissipative nature of the collective dynamics. For $J=0$, the dynamics (equation (5)) would preserve phase space volume. This volume, however, is strongly contracted by spikes received in the post-synaptic neurons. Consider trajectories from a small ball of initial conditions as they emit the same future spike. The ball of phases at this spike contracts by a factor $1-d$ along each of the $K$ dimensions of the subspace spanned by the post-synaptic neurons. The volume thus contracts by $(1-d)^{K}\to e^{\lambda_{\mathrm{inh}}}$ per spike, for $K\gg 1$, with exponential rate,

[TABLE]

$\lambda_{\mathrm{inh}}<0$ is responsible for the linear stability of the dynamics given by equations (1) and (2), first shown in Refs. 22; 16.

Phase-space partitioning

The phase space volume taken up by an ensemble of nearby trajectories at a given spike is contracted at the spike’s reception. Larger phase space volumes, however, are not uniformly contracted but torn apart, with the pieces individually contracted and overall dispersed across the entire traversed phase space volume. The elementary phenomenon is illustrated in Fig. 1.

We define the critical perturbation strength, $\epsilon^{*}$, as the flux tube’s extent out from a given state $\vec{\phi}_{0}$ on the equilibriated trajectory, $\vec{\phi}_{t}$, and in a given orthogonal perturbation direction, $\vec{\xi}$,

[TABLE]

Here, $D_{t}\left(\epsilon\right)$ is the 1-norm distance,

[TABLE]

between $\vec{\phi}_{t}$ and the perturbed trajectory, $\vec{\phi}_{t}(\epsilon)$, evolving freely from the perturbed state, $\vec{\phi}_{0}\left(\epsilon\right):=\vec{\phi}_{0}+\epsilon\vec{\xi}$ (reference time $t=0$ and $||\vec{\xi}||=1$; see Supplemental Methods for details). $\epsilon^{*}$ is the largest value below which $D_{t}\left(\epsilon\right)$ vanishes in time. $D_{t}$ initially decays exponentially, but for a supercritical perturbation, $\epsilon>\epsilon^{*}$, there exists a divergence event time, $t^{*}>0$, defined and obtained as the time at which a sustained divergence in $D_{t}$ begins (see Fig. 1a).

A 2D cross-section of the phase space around $\vec{\phi}_{0}$ (Fig. 1b) reveals that the locations of these critical perturbations form lines which intersect to form polygon-shaped basin boundaries. Before developing a theory for this phase space organization (caricatured in Fig. 1c), we first analyze two main features of the geometry of a flux tube: the punctuated exponential decay of its cross-sectional volume and the exponential separation of neighboring tubes.

Punctuated geometry of flux tubes

As expected from the typical phase space volume contraction (see equation (8)), we find along a simulated trajectory that the orthogonal phase space volume enclosed by the local flux tube exhibits exponential decay. This decay, however, is punctuated by blowup events. Figure 2a displays the spiking activity produced by the typical trajectory, $\vec{\phi}_{t}$. The neighborhood around $\vec{\phi}_{t}$ over a time window is visualized in a folded representation using a fixed, 2D projection of the phase space (Fig. 2b and Supplementary Video; see Supplemental Information for construction details). The basin of attraction surrounding $\vec{\phi}_{t}$ (Fig. 2b) consists of lines which remain fixed between spike times. Across spike times, new lines appear and existing lines disappear. At irregular intervals breaking up time windows of exponential contraction, large abrupt blowup events take the boundary away from the center trajectory (Fig. 2b, c), producing jumps in the area enclosed by the boundary. It is important to note that these events do not mean that the evolving phase space volume from an ensemble of states contained in the tube would expand. Such volumes only contract and converge to the same asymptotic trajectory. The basin of attraction itself, however, does not exclusively contract with time. In fact, it should on average maintain a typical size.

The blowup events typically coincide with a divergence event time, $t^{*}$ (Fig. 1a), in some perturbation direction. Two such coincidences are visible in Fig. 2c,d. We conclude that the local basin at any time extends out in phase space until the perturbed trajectory approaches the pre-image of a divergence event occurring at a future time. Flux tube shape is then determined by the statistics of such events.

Tube boundary and divergence

We analyzed a set of divergence events from simulations. We find that a collision of a pair of spikes constitutes the elementary event triggering the divergence of the perturbed trajectory. These pairs, hereon called susceptible spike pairs, were generated by connected pairs of neurons. Moreover, a perturbation-induced collision of a susceptible spike pair generated an abrupt spike time shift in one or both of the connected neurons’ spike times. We found that the nature of the spike time shift depends on the motif by which the two neurons connect. We denote the backward-connected pair motif $n_{s^{*}}\leftarrow n_{s^{\prime}{}}$, where $s^{*}$, the decorrelation event index, is the spike index of the earlier of the pair (note that $t^{*}\equiv t_{s^{*}}$), and $s^{\prime}{}>s^{*}$ labels the later spike in the pair.

For $\epsilon\lesssim\epsilon^{*}$, the presynaptic spike time, $t_{s^{\prime}{}=s^{*}+1}$, is advanced with increasing $\epsilon$ relative to the postsynaptic spike time, $t_{s^{*}}$, until the two spikes collide (see Fig. 3). At collision ($\epsilon=\epsilon^{*}$), the pulsed inhibition and the rate of approach to voltage threshold cause an abrupt delay of $t_{s^{*}}$ by $\Delta t_{\mathrm{jump}}$. Using equation (3), we obtain

[TABLE]

for $d\ll 1$. Further details and the other two motifs (forward-connected and symmetrical) are discussed in the Supplementary Notes.

For each spike in the network sequence, the rate of its susceptible spike partners is

[TABLE]

where $\mkern 1.5mu\overline{\mkern-1.5mu\Delta t\mkern-1.5mu}\mkern 1.5mu=\left(N\bar{\nu}\right)^{-1}$ is the average distance between successive spikes. Since $\Delta t_{\mathrm{jump}}\approx\lambda_{\mathrm{sus}}^{-1}$, the spike time of neuron $n_{s^{*}}$ is shifted forward typically as far as its next nearest susceptible partner spike. Thus, one collision event will typically induce another in at least one of the $\mathcal{O}\left(K\right)$ neurons to which the involved pair of neurons are presynaptic. A cascade of collision events then follows with near certainty (see Supplemental Notes for details).

The shift in $t_{s^{*}}$ by $\Delta t_{\mathrm{jump}}$ is carried forward to all future spike times of $n_{s^{*}}$, so that $n_{s^{*}}$ becomes a source of collision events. The total collision rate is then $\lambda_{\mathrm{sus}}$ multiplied by the number of source neurons, which approximately increments with each collision in the cascade. Averaging over realizations of the cascade (reference time $t^{*}=0$), the average number of collisions, $\bar{m}$, grows as $\tfrac{\mathrm{d}}{\mathrm{d}t}\bar{m}=\lambda_{\mathrm{sus}}\bar{m}$. Finally, since each collision produces a jump in distance of equal size, we obtain the pseudoLyapunov exponent, $\lambda_{p}=\lambda_{\mathrm{sus}}$ from its implicit definition, $\bar{D}_{t}=\bar{D}_{0}\exp\left[\lambda_{p}t\right]$ (see Supplemental Notes), as the exponential rate at which flux tubes diverge.

Statistical theory of flux tube diameter

The geometry of a flux tube is captured by the flux tube indicator function, $\mathbb{1}_{\mathrm{FT}}\left(\epsilon\right)=\Theta\left(\epsilon^{*}(\vec{\phi}_{0},\vec{\xi})-\epsilon\right)$, evaluated across network states, $\vec{\phi}_{0}$, of its contained attracting trajectory and perturbation directions, $\vec{\xi}$. Using the Heaviside function, $\Theta(x)$, $\mathbb{1}_{\mathrm{FT}}\left(\epsilon\right)=1$ for perturbations remaining in the tube ($\epsilon<\epsilon^{*}$), and 0 otherwise. The average of $\mathbb{1}_{\mathrm{FT}}\left(\epsilon\right)$ over $\vec{\phi}_{0}$ and $\vec{\xi}$ ,

[TABLE]

is the survival function: the probability that an $\epsilon$-sized perturbation does not lead to a divergence event later in the perturbed trajectory, i.e. $\epsilon<\epsilon^{*}$, and is formally defined as $\hat{S}\left(\epsilon\right):=1-\int_{0}^{\epsilon}\rho\left(\epsilon^{*}\right)\mbox{d}\epsilon^{*}$, with $\rho\left(\epsilon^{*}\right)$ the transformed density over $\epsilon^{*}$. $\hat{S}\left(0\right)=1$ and decays to 0 as $\epsilon\to\infty$. The scale of this decay defines the typical flux tube size. Calculating $\hat{S}\left(\epsilon\right)$ requires two steps: firstly, establishing a tractable representation of $\epsilon^{*}(\vec{\phi}_{0},\vec{\xi})$ and secondly, performing the average in equation (13). Both of these in general pose intricate problems. However, as we will see next, both substantially simplify when generic properties of the asynchronous, irregular state are taken into account.

Perturbed spike intervals are obtained using the spike time deviations, $\delta t_{s}\left(\epsilon\right):=t_{s}\left(\epsilon\right)-t_{s}\left(0\right)$, $s=1,2,\dots$,

[TABLE]

In a linear approximation we find,

[TABLE]

where $C=\frac{T_{\mathrm{free}}}{\sqrt{N}}$ converts network phase deviation to spike time deviation and $a_{s}$ is a dimensionless susceptibility that depends on the adjacency matrix, $A=\left(A_{mn}\right)$, derivatives of the phase response curve evaluated at the network states at past spike times, $\{\vec{\phi}_{s^{\prime}}=\vec{\phi}_{t_{s^{\prime}}}\}$ for $s^{\prime}<s$, and the perturbation direction $\vec{\xi}$ (see Supplemental Notes for its derivation). Substituting equation (15) into equation (14) gives

[TABLE]

with $\Delta t_{s}=\Delta t_{s}\left(0\right)$. Note that $\Delta t_{s}\left(\epsilon\right)$ can have a zero, i.e. a spike time collision only when $\Delta a_{s}=a_{s}-a_{s-1}>0$.

To obtain the scaling behavior of the flux tube geometry it is sufficient to examine the statistics of flux tube borders using the corresponding divergence events generated by collisions of backward-connected susceptible spike pairs in the perturbed trajectory (Fig. 3). In these cases, the perturbation strength $\epsilon\to\epsilon^{*}$ as the network spike interval $\Delta t_{s^{*}}\left(\epsilon\right)\to 0$ for $A_{n_{s^{*}}n_{s^{*}+1}}=1$. In fact, the latter condition serves in these cases as an implicit definition of $\epsilon^{*}$ and $s^{*}$.

According to equation (16), $\hat{S}\left(\epsilon\right)$ in principle depends on the adjacency matrix, $A=\left(A_{mn}\right)$, of the network realization. Removing this dependence by averaging over the ensemble of graphs, $P_{A}\left(\left(A_{mn}\right)\right)$, simplifies the calculation of the survival function,

[TABLE]

Evaluating the right-hand side of equation (17) using the perturbed spike intervals, linearized in $\epsilon$, requires knowledge of the joint probability density of all variables present in equation (16),

[TABLE]

where we have chosen the perturbation direction, $\vec{\xi}$, to be statistically independent of the state, $\vec{\phi}_{0}$, being perturbed at $t=0$. Here, the unperturbed spike pattern is represented by two random variables: $M$, the number of spikes in the time interval $\left[0,T\right]$ after the perturbation, and $\left\{\Delta t_{s}\right\}$, the set of all $M-1$ inter-spike intervals in this window. It is well understood that in the large-system limit in a sparse graph, $1\ll K\ll N$, the currents driving individual neurons in the network converge to independent, stationary Gaussian random functions 23. For low average firing rates, this implies that the pattern of network spikes $\left(M,\left\{\Delta t_{s}\right\}\right)$ resembles a Poisson process 24. Furthermore, the susceptibility becomes state-independent in this limit. Neglecting the weak dependence between the distribution of network spike patterns and $A=\left(A_{mn}\right)$, the full density, $\rho_{T}$ (equation (18)), approximately factorizes,

[TABLE]

with distribution of a single adjacency matrix element, $P_{A_{mn}}\left(A_{mn}=1\right)=p$, $P_{A_{mn}}(A_{mn}=0)=1-p$, count distribution of spikes in the observation window, $P_{T}(M)$, and distribution of single inter-spike interval $\rho_{t}\left(\Delta t\right)$. The latter is exponential with rate $\mkern 1.5mu\overline{\mkern-1.5mu\Delta t\mkern-1.5mu}\mkern 1.5mu$. With these approximations (see Supplementary Notes for details), all dependencies on the distribution of perturbation direction are mediated by the susceptibilities, $\left\{\Delta a_{s}\right\}$. For any isotropic $\rho(\vec{\xi})$ having finite-variance, $\rho_{a}\left(\Delta a_{s}\right)$ has zero mean and standard deviation proportional to $\exp\left[\frac{\lambda_{\mathrm{inh}}}{N}s\right]$, with the average contraction rate per neuron, $\frac{\lambda_{\mathrm{inh}}}{N}=-\frac{Kd}{N}=-pd$, due to the inhibition. The factor $2\Theta(\Delta a_{s})$ places support only positive values of $\Delta a_{s}$ as required.

As $\rho_{T}$ factorizes, so does $S\left(\epsilon\right)$,

[TABLE]

where $S_{s}\left(\epsilon\right)$ is the probability that a perturbation of strength $\epsilon$ does not lead to a collision event involving the $s^{\mathrm{th}}$ spike. With the above simplifications,

[TABLE]

Evaluating equation (21) (see Supplementary Notes for the derivation), we find

[TABLE]

where $x_{s}=\frac{C}{\mkern 1.5mu\overline{\mkern-1.5mu\Delta t\mkern-1.5mu}\mkern 1.5mu}e^{\frac{\lambda_{\mathrm{inh}}}{N}s}\epsilon\leq\epsilon/\sqrt{p}$, and $\mathrm{Erfcx}\left[x\right]=e^{x^{2}}\left(1-\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-y^{2}}\mathrm{d}y\right)$ is the scaled complementary error function. $\mathrm{Erfcx}\left[x_{s}\right]-1\approx-x_{s}$ for $\epsilon/\sqrt{p}\ll 1$, so that finally

[TABLE]

where we have identified $\lambda_{\mathrm{sus}}=\frac{p}{\mkern 1.5mu\overline{\mkern-1.5mu\Delta t\mkern-1.5mu}\mkern 1.5mu}$. Employing the logarithm and $C\lambda_{\mathrm{sus}}\epsilon\propto\sqrt{p}\epsilon\ll 1$,

[TABLE]

with

[TABLE]

where we have used equations (6) and (7) in the second line and note the cancellation of $p$ and $I_{\mathrm{Ext}}$. Equation (23) shows for $1\ll K\ll N$ that the basin diameter, $\epsilon^{*}$, is exponentially distributed and so completely determined by its characteristic scale, $\mkern 1.5mu\overline{\mkern-1.5mu\epsilon^{*}\mkern-1.5mu}\mkern 1.5mu$ (equation (25)), that is smaller for larger network size, higher average in-degree, higher population activity, and larger membrane time constant, $\tau$. $\mkern 1.5mu\overline{\mkern-1.5mu\epsilon^{*}\mkern-1.5mu}\mkern 1.5mu$ grows, however, with the synaptic coupling strength, $J_{0}$. In Fig. 4b, we see quantitative agreement in simulations between the definition of $\hat{S}\left(\epsilon\right)$ (equation (13) using the definition of $\epsilon^{*}$, equation (9)) and its approximate microstate parametrization (equations (20), (21)). These also confirm the exponential form of our reduced expression (equations (23), (25)) and a scaling dependence on $J_{0}$ (Fig. 4c). The latter holds until $J$ is no longer of size $\mathcal{O}\left(1/\sqrt{K}\right)$. The other scalings were reported in Ref.17. A derivation of only the characteristic scaling of $\mkern 1.5mu\overline{\mkern-1.5mu\epsilon^{*}\mkern-1.5mu}\mkern 1.5mu$, but not depending on the Poisson spiking assumption, is given in the Supplemental Notes.

The geometry of phase space partitioning

Figure 5 presents the phase space organization of these spiking circuits as we have revealed it, replacing the caricature of Fig. 1c. For a perturbation made to a stable trajectory, the geometry of the determining collision event is shown in Figure 5a, in a folded representation. The pre-images of this event determine the flux tube boundary back to the perturbation. Our results also provide a global, i.e. non-folded geometry of the partitioning (Fig. 5b(left)). Susceptible spike collisions are edges of the $N$-dimensional unit hypercube of phases where the corresponding voltages of two connected neurons both approach threshold. The Poincare section obtained by projecting the dynamics orthogonal to the trajectory (since no motion exists orthogonal to this subspace) then reveals the intrinsic partition. Here, the polygon basin boundaries arise as the pre-images of the projections of susceptible edges lying nearby the trajectory at future spike times (Fig. 5b(right)).

{discussion}

We have developed a theory of phase space partitioning in spiking neural circuits, exemplified using the phenomenon of flux tubes. Importantly, the approach yields the dependence on various control parameters. We find the flux tube diameter contracts with the rate of volume contraction per neuron, $\lambda_{\mathrm{inh}}/N=(N\bar{\nu}\tau)^{-1}$, due to the inhibition received across the post-synaptic subspace of each spike. This contraction is punctuated, however, by collision events between susceptible spikes, i.e. those from pairs of connected neurons, occurring at rate $\lambda_{\mathrm{sus}}=(K\bar{\nu})^{-1}$ and across which the basin volume expands out to a pre-image of the next collision event. For some neighboring tube, this collision event sets off a cascade of such events with exponential rate, $\lambda_{\mathrm{sus}}$ that is responsible for their mutual divergence. Using these collision events to identify the spiking trajectories lying on flux tube boundaries, we were able calculate the size distribution of these basins. The average size is controlled by the ratio of these two exponential rates. Leaving out a factor converting shifts in spike time to shifts in state,

[TABLE]

The final scaling, $\mkern 1.5mu\overline{\mkern-1.5mu\epsilon^{*}\mkern-1.5mu}\mkern 1.5mu=J_{0}/\left(\sqrt{KN}\bar{\nu}\tau\right)$, thus combines the contraction from the single neuron dynamics responsible for the dissipative dynamics, with the overall rate of spikes, which appears since each spike can be involved in a destabilizing collision event. Both contracting and expanding rates scale with the probability of connection, $p$, so we intuitively expect $p$ to appear in $\mkern 1.5mu\overline{\mkern-1.5mu\epsilon^{*}\mkern-1.5mu}\mkern 1.5mu$ only implicitly through $J$ and, reassuringly, $p$ indeed cancels out.

Our framework motivates a variety of extensions. Our calculations can be performed for different disordered connectivity ensembles (e.g. correlated entries from annealed dilution processes 25 and structured second-order statistics 26), different activity regimes (e.g. non-Markovian spike interval processes 27), and different single neuron models (e.g. any threshold neuron with known phase response curve). We have applied the theory to an instability caused by abrupt changes in spike time due to an inhibitory input near voltage threshold, a scenario that can also be analyzed in neuron models with smooth thresholds (e.g. the rapid theta-neuron 28 that has the LIF neuron as a limit). The theory may also apply to other, as yet unknown instabilities involving spike collision events. Finally, while the linear stability of the dynamics precludes finite, asymptotic (Kolmogorov-Sinai) entropy production, the partition refinement picture we provide in Fig. 5b suggests a transient production of information about the perturbation on timescales of the order of the divergence event time, $t^{*}$. Making this connection to ergodic theory more precise is an interesting direction for future research.

Applying our approach in a relatively idealized context allowed for a tractable assessment of phase space organization. Despite its simplicity, however, the LIF neuron accurately captures many properties of cortical neurons, such as their dynamic response 29. We have also neglected heterogeneity in many properties. For instance, in contrast to the locally stable regime studied here, mixed networks of excitatory and inhibitory neurons can instead be conventionally chaotic 30. This chaos can nevertheless be suppressed in the ubiquitous presence of fluctuating external drive 31; 32 or with spatially-structured connectivity 33, suggesting a generality to locally stable dynamics and phase space partitioning in neural computation. Our approach, in particular the way we have quantified the ensemble of perturbed spiking trajectories, can inform formulations of local stability in these more elaborate contexts. Of particular interest are extensions where a macroscopic fraction of tubes remain large enough to realize encoding schemes tolerant of intrinsic and stimulus noise. For example, extensions to random dynamical systems 34; 35 could provide theoretical control over spiking dynamic variants of rate network-based learning schemes to generate stable, input-specific trajectories 7.

Recent advances in experimental neuroscience have allowed for probes of the finite-size stability properties of cortical circuit dynamics call for in vivo. For example, simultaneous intra- and extra-cellular recordings in the whisker motion-sensing system of the rat reveal that the addition of a single spike makes a measurable impact on the underlying spiking dynamics of the local cortical area 36. Indeed, rats can be trained to detect perturbations to single spikes emitted in this area 37. Representative toy theories, such as the one we provide, can guide this work by highlighting the features of spiking neural circuits that contribute to phase space partitioning. The combined effort promises to elucidate the dynamical substrate for neural computation at the level at which the neuronal interactions actually operate.

Acknowledgements

M.P.T. would like to acknowledge discussions with Michael Monteforte, Sven Jahnke and Rainer Engelken. This work was supported by BMBF (01GQ07113, 01GQ0811, 01GQ0922, 01GQ1005B), GIF (906-17.1/2006), DFG (SFB 889), VW-Stiftung (ZN2632), and the Max Planck Society.

Author Contributions

M.P.T. conceived the project, developed the concepts, and wrote the manuscript. F.W. supervised the project, discussed the results and edited the manuscript.

Additional Information

The authors declare no competing financial interests.