Cusps Enable Line Attractors for Neural Computation

TL;DR

This paper links line attractors in neuronal networks to pulse-gating mechanisms, showing how they enable robust graded information propagation through a cusp catastrophe in high-dimensional, stochastic neural models.

Contribution

It reveals the dynamical origin of line attractors via cusp catastrophes in pulse-gated neural networks, connecting theoretical analysis with biological plausibility.

Findings

Line attractors are associated with cusp catastrophes in neural dynamics.

Ghost dynamics near bifurcations support robustness of graded information.

Pulse gating, synaptic noise, and stochasticity enable stable line attractors.

Abstract

Line attractors in neuronal networks have been suggested to be the basis of many brain functions, such as working memory, oculomotor control, head movement, locomotion, and sensory processing. In this paper, we make the connection between line attractors and pulse-gating in feedforward neuronal networks. In this context, because of their neutral stability along a one-dimensional manifold, line attractors are associated with a time-translational invariance that allows graded information to be propagated from one neuronal population to the next. To understand how pulse-gating manifests itself in a high-dimensional, non-linear, feedforward integrate-and-fire network, we use a Fokker-Planck approach to analyze system dynamics. We make a connection between pulse-gated propagation in the Fokker-Planck and population-averaged mean-field (firing rate) models, then identify an approximate line attractor in state space as the essential structure underlying graded information propagation. An analysis of the line attractor shows that it consists of three fixed points: a central saddle with an unstable manifold along the line and stable manifolds orthogonal to the line, which is surrounded on either side by stable fixed points. Along the manifold defined by the fixed points, slow dynamics give rise to a ghost. We show that this line attractor arises at a cusp catastrophe, where a fold bifurcation develops as a function of synaptic noise; and that the ghost dynamics near the fold of the cusp underly the robustness of the line attractor. Understanding the dynamical aspects of this cusp catastrophe allows us to show how line attractors can persist in biologically realistic neuronal networks and how the interplay of pulse gating, synaptic coupling and neuronal stochasticity can be used to enable attracting one-dimensional manifolds and thus, dynamically control the processing of graded information.

Figures

A), and graded current and firing rate (Fig. 1B) may be propagated via pulse-gating between layers in the Fokker-Planck model. In Fig. 1B inset, we compare currents between a time-averaged firing rate model of pulse-gated graded information propagation and the Fokker-Planck model. For a given value of $T$ in a mean-field rate model, we can find an exact $S$. Because of the slow onset of firing (blue, inset) relative to the firing rate model (dashed red, inset), the value of synaptic coupling that gave rise to graded transfer, $S=S_{graded}$, was larger than the mean-field prediction, $S_{exact

A,B, we show firing rates and corresponding current amplitudes for a range of input currents in the Fokker-Planck model. For all but the lowest input current, the output current (the current at the end of the gating pulse) is very close to the input current. Plotting input versus output current for ranges of gating currents, $\bar{I}^{g}$ (Fig. 2C), and synaptic couplings, $S$ (Fig. 2D), we find that changes in gating currents translate the input-output function upwards and changes in synaptic coupling change the slope of the input-output function. By varying $\bar{I}^{g}$ and $S$, we can find

A,B, we compare numerical and approximate (Gaussian) solutions. Note the similar advection of the population density with obvious differences in curvature and reset populations (which may be seen in Fig. 5A as a linear wall of density along the reset value, [math], of the membrane potential). In Fig. 5C, we show the fold bifurcation, found numerically in Fig. 4, for a range of values of $S$ and $\bar{I}^{\mathrm{g}}$. In Fig. 5D, we demonstrate the existence of a cusp catastrophe as a function of the width of the potential distribution $\sigma$. In Fig. 5E,F, we show how fast or slow (ghost) d

Equations

$$\frac{d}{dt}{V}_{i,j}$$

$$\tau \frac{d}{dt}{I}_{i,j}^{\mathrm{ff}}$$

$$\frac{\partial }{\partial t}{\rho }_j\left(V,t\right)$$

$$\tau \frac{d}{dt}{\overline{I}}_j^{\mathrm{ff}}$$

$$J_j\left(V,t\right)=\left(\left[-g_L\left(V-V_{\mathrm{R}}\right)+{\bar{I}}^{\mathrm{g}}+{\overline{I}}_j^{\mathrm{ff}}\right]-D_j\frac{\partial }{\partial V}\right){\rho }_j\left(V,t\right).$$

$$D_j={\sigma }_0^2+\frac{1}{2}\frac{S^2}{pN}{m}_{j-1}\left(t\right).$$

$$m_j\left(t\right)=J_j\left(V_{\mathrm{Th}},t\right).$$

$$\mu (t)={\mu }_0{e}^{-g_{L}t}+\frac{{\bar{I}}^{\mathrm{g}}}{g_L}(1-e^{-g_{L}t})+\frac{A}{\frac{1}{\tau }-g_L}(e^{-g_{L}t}-e^{-t/\tau })$$

$$m(t)={\left[(-g_L\mu (t)+{\bar{I}}^{\mathrm{g}}+A{e}^{-t/\tau })\frac{1}{P}{e}^{-(1-\mu (t){)}^2/2{\sigma }^2}\right]}^{+},$$

$${\bar{I}}_d^{\mathrm{ff}}=S{e}^{-T/\tau }{\int }_0^T{e}^{t/\tau }m(t)\frac{dt}{\tau }.$$

$$\tau {\frac{d{\bar{I}}_j}{dt}}^{ff}$$

$$m_j$$

$$S_{exact}=\frac{\tau }{T}{e}^{T/\tau }$$

$$I^{ff}_{j}(t)=\left\{\begin{array}[]{ll}A\left(\frac{t-(j-1)T}{\tau}\right)e^{-(t-(j-1)T)/\tau},&(j-1)T\leq t\leq jT\\ A\left(\frac{T}{\tau}e^{-T/\tau}\right)e^{-(t-jT)/\tau},&jT<t<\infty\end{array}\right.$$

$$m_{j}(t)=\left\{\begin{array}[]{ll}0,&0<t<(j-1)T\\ A\left(\frac{T}{\tau}e^{-T/\tau}\right)e^{-(t-jT)/\tau},&(j-1)T\leq t\leq jT\\ 0,&jT<t<\infty\end{array}\right.$$

Full text

I Introduction

Attractor neural networks have been appealed to by theoretical neuroscientists to explain working memory (Wang, 2001), oculomotor control (Seung, 1996), head movement (Zhang, 1996), locomotion (Bruno et al., 2015), sensory processing (Ben-Yishai et al., 1995) and many other experimentally observed brain functions. Some of these neuronal functions have been rationalized using the persistent activity of an attractor (e.g., working memory), while other functions have been modeled by the gradedness (i.e., order preserving property (Machens et al., 2005)) of an attractor (e.g., oculomotor control and sensory processing). Because of these dynamical features, attractor networks have been particularly useful in modeling the encoding of continuous external stimuli (Mante et al., 2013) or the maintenance of internal neuronal representations (Machens et al., 2005). Furthermore, it has been shown that attractor network dynamics have the capacity to perform optimal computations (Latham et al., 2003) and thus, can implement Bayesian inference (Wu et al., 2003; Deneve et al., 2007).

Pulse-gating is a mechanism capable of transferring packet-based spiking activity from one neuronal population to another (Sornborger et al., 2015; Wang et al., 2016). For instance, the synfire-gated synfire chain (SGSC) (Wang et al., 2016) is a mechanism consisting of two sets of neural populations, one representing rate-coded (graded) information and another that gates the flow of that information. This separation into information-carrying and information-control populations provides a basis for understanding how some brain areas can be responsible for the control of information transfer (Gisiger and Boukadoum, 2011), such as the mediodorsal nucleus of the thalamus (Floresco and Grace, 2003), and others for information processing, such as the somatosensory areas (Hernandez et al., 2000).

The original information propagation model that demonstrated the concept of pulse-gated firing rate propagation was based on a mean-field firing rate model of a network of current-based integrate-and-fire (I&F) neurons that was coarse-grained in time. There was good correspondence between the mean-field model and mean spiking rates in a network of I&F neurons (Sornborger et al., 2015; Wang et al., 2016). However, in general, mean-field models make critical use of the relationship between firing rate and current (f-I curves) (or other input-output functions), which correspond to the steady-state responses of a population of neurons driven by a constant input. There is no a priori reason that a temporally-averaged mean-field model should correspond with an I&F network in a transiently driven (pulse-gating) context. Furthermore, the original mean-field solution that demonstrated pulse-gated information propagation suffered from parameter fine-tuning (Seung, 1996; Goldman, 2008; MacNeil and Eliasmith, 2011), even though I&F simulations appeared robust to parameter variability.

In this paper, in order to better understand the robustness of pulse-gated, graded information transfer in a network of spiking neurons, we construct Fokker-Planck equations describing the membrane potential probability density function in a feedforward network and study solutions in the state space of a dimensionally-reduced iterative dynamical system. After outlining our methods, we show that pulse-gating in feedforward networks gives rise to approximately time-translationally-invariant spiking probabilities that are propagated from layer to layer in a feedforward network. We then examine synaptic current input-output relations allowing for the construction of an effective population firing rate model and show that this model, averaged across neuronal populations, is very similar to a mean-field model based on rectified linear input-output functions. We show how the input-output function depends on synaptic coupling, gating pulse amplitude, and synaptic noise, and that there is a sizeable region of parameter space within which graded propagation exists.

We then demonstrate the dynamical convergence of membrane potential density to a one-dimensional manifold in parameter space. The Fokker-Planck system reveals that this one-dimensional manifold arises due to a saddle node bifurcation creating an unstable fixed point in the center of a one-dimensional manifold and two stable fixed points on either side of the saddle node. Dynamics orthogonal to the one-dimensional manifold give rise to rapid convergence to the one-dimensional manifold. Ghost (slow) dynamics along the one-dimensional manifold allow the unstable manifold of the saddle to be viewed as an attracting one-dimensional manifold.

We finally show that this approximate line attractor is robust and generic. By using a reduced analytical model of gating induced transients in the Fokker-Planck system, we show that the propagation of firing amplitudes can be mapped to a cusp catastrophe, and that a bundle of one-dimensional ghost manifolds exists in the region surrounding the fold of the cusp. Our results reveal how the coordination of pulse-gating, synaptic coupling, and membrane potential dynamics enables approximate line attractors in feedforward networks and demonstrate the robustness of graded information propagation when pulse-gating is incorporated.

II Methods

We study a feedforward network of $j=1,\dots,M$ populations of $i=1,\dots,N$ excitatory, current-based, integrate-and-fire (I&F) neurons whose membrane potential, $V$, and synaptic current are described by

[TABLE]

where $V_{\mathrm{R}}$ is the rest voltage (also the reset voltage), $\tau$ is the synaptic timescale, $S$ is the synaptic coupling strength, $p_{jk}$ is a Bernoulli distributed random variable and $p=\langle p_{jk}\rangle$ is the mean synaptic coupling probability. The $l$’th spike time of the $k$’th neuron in layer $j-1$ is determined by $V_{i,j}(t^{l}_{j-1,k})=V_{\mathrm{Th}}$, i.e. when the neuron reaches threshold (after which $V_{i,j}$ is immediately reset to $V_{R}$). The gating current, $I^{g}_{j}$, is a white noise process with a square pulse envelope, $\theta(t-(j-1)T)-\theta(t-jT)$, where $\theta$ is a Heaviside theta function and $T$ is the pulse length (Sornborger et al., 2015) of pulse height $\bar{I}^{\mathrm{g}}$ and variance $\sigma_{0}^{2}$. Note that an exponentially decaying current is injected in population $1$ providing synchronized activity that will subsequently propagate downstream through populations $j=2,\dots,M$.

Feedforward networks of this type admit discrete-time-translationally-invariant solutions describing graded packet transfer Sornborger et al. (2015); Wang et al. (2016). Such solutions exist for gating pulses that are either temporally sequential or overlapping Wang et al. (2016).

To understand this phenomenon, we used a Fokker-Planck analysis. Assuming the spike trains in Eq. (1d) to be Poisson distributed, the collective behavior of this feedforward network may be described by the (Fokker-Planck) equations

[TABLE]

These equations describe the evolution of the probability density function, $\rho_{j}(V,t)$, in terms of the probability density flux, $J_{j}(V,t)$, the mean feedforward synaptic current, $\bar{I}^{\mathrm{ff}}_{j}$, and the population firing rate, $m_{j}$. For each layer, $j$, the probability density function gives the probability of finding a neuron with membrane potential $V\in(-\infty,V_{\mathrm{Th}}]$ at time $t$.

The probability density flux is given by

[TABLE]

The effective diffusivity is

[TABLE]

(Below, we take $N\rightarrow\infty$ and neglect the second term on the right in (3).) The population firing rate is the flux of the probability density function at threshold,

[TABLE]

The boundary conditions for the Fokker-Planck equations are $J_{j}\left({V_{\mathrm{R}}^{+},t}\right)=J_{j}\left({{V_{\mathrm{Th}}},t}\right)+J_{j}\left({V_{\mathrm{R}}^{-},t}\right)$, ${\rho_{j}}\left({V_{\mathrm{R}}^{+},t}\right)={\rho_{j}}\left({{V_{\mathrm{Th}}},t}\right)+{\rho_{j}}\left({V_{\mathrm{R}}^{-},t}\right)$, and ${\rho_{j}}\left({V=-\infty,t}\right)=0$.

To improve the efficiency of exploring the bifurcation structure of this system in a large state space, we also investigate an approximate model in which the initial distribution is assumed to be Gaussian, $\rho_{0}(V,t)=(1/P)\exp{(-(V-\mu(t))^{2}/2\sigma^{2})}$, with width $\sigma$ and mean $\mu(t)$, where $P=\int_{-\infty}^{V_{\mathrm{Th}}}\rho_{0}(V,0)$ is a normalization factor accounting for the truncation of the Gaussian at threshold, $V_{\mathrm{Th}}$. As the gating current turns on, the distribution is uniformly advected toward the voltage threshold, $V_{\mathrm{Th}}$, and the population begins to fire. Uniform advection neglects a small amount of firing due to a diffusive flux across the firing threshold, thus the fold bifurcation (see Results) occurs at a slightly larger value of synaptic coupling, $S$, for this approximation relative to numerical simulations. Since the timescale of the pulse is fast, neurons only have enough time to fire once (approximately). Thus, we neglect the re-emergent population at $V_{\mathrm{R}}$, which contributes negligibly to firing during the transient pulse.

With this approximation, Eq. (2a) gives rise to $\dot{\mu}=-g_{L}(\mu-V_{\mathrm{R}})+\bar{I}^{\mathrm{g}}+\bar{I}_{u}^{\mathrm{ff}}$, where $\sigma^{2}=\sigma_{0}^{2}/g_{L}$. With upstream current $\bar{I}_{u}^{\mathrm{ff}}=Ae^{-t/\tau}$. Setting $V_{\mathrm{Th}}=1$, this integrates to

[TABLE]

and from Eq. (4), we have

[TABLE]

which, from Eq. (2d), results in a downstream synaptic current at $t=T$

[TABLE]

After the end of the pulse, the current decays exponentially. This decaying current feeds forward and is integrated by the next layer. Thus, for an exact transfer, $\bar{I}_{d}^{\mathrm{ff}}(S,\bar{I}^{\mathrm{g}},A,T)=A$.

III Results

In Fig. 1, we show that population density (Fig. 1A), and graded current and firing rate (Fig. 1B) may be propagated via pulse-gating between layers in the Fokker-Planck model. In Fig. 1B inset, we compare currents between a time-averaged firing rate model of pulse-gated graded information propagation and the Fokker-Planck model. For a given value of $T$ in a mean-field rate model, we can find an exact $S$. Because of the slow onset of firing (blue, inset) relative to the firing rate model (dashed red, inset), the value of synaptic coupling that gave rise to graded transfer, $S=S_{graded}$, was larger than the mean-field prediction, $S_{exact}$ (see Appendix), by a factor of $1.07$ for $T$ near $5$ ms. In Fig. 1C, we demonstrate stable graded propagation across many layers.

In Fig. 2A,B, we show firing rates and corresponding current amplitudes for a range of input currents in the Fokker-Planck model. For all but the lowest input current, the output current (the current at the end of the gating pulse) is very close to the input current. Plotting input versus output current for ranges of gating currents, $\bar{I}^{g}$ (Fig. 2C), and synaptic couplings, $S$ (Fig. 2D), we find that changes in gating currents translate the input-output function upwards and changes in synaptic coupling change the slope of the input-output function. By varying $\bar{I}^{g}$ and $S$, we can find an optimal input-output function very close to the diagonal (Fig. 2E), hence giving very accurate propagation of graded information. A nearby basin of parameters also gives good graded propagation (Fig. 2F) as measured by the distance of the input-output function from the diagonal. The input-output function in Fig. 2E may be used in a population firing rate model for large pulse-gated systems for which spiking networks are impractical.

In the graded transfer regime, the probability density may be described by a few of its low-lying moments, viz. $M^{1}_{j}(t)=\int_{-\infty}^{V_{\mathrm{Th}}}V\rho_{j}(V,t)dV$ and $M^{2}_{j}(t)=\int_{-\infty}^{V_{\mathrm{Th}}}V^{2}\rho_{j}(V,t)dV$. Changes in the moments are related to changes in the shape of the density. Higher moments rapidly converge to the one-dimensional manifold (not shown). To identify time-translationally invariant solutions admitting graded propagation, we search for parameters, $S$, $A$, etc. such that for a given $T$, $\bar{I}_{j}^{\mathrm{ff}}(t+T)\approx\bar{I}_{j-1}^{\mathrm{ff}}(t)$.

In Fig. 3, we plot the map $(m_{j}(jT),M^{1}_{j}(jT),M^{2}_{j}(jT))$ for successive values of $j$ (i.e. mapping feedforward propagation from layer to layer). As may be seen in Fig. 3, for parameters admitting graded information propagation, initial conditions rapidly approach a one-dimensional manifold. On the manifold itself, there are three fixed points, two stable fixed points at the extremes and one unstable saddle in between. For these parameters, the dynamics along the unstable direction of the saddle are slow. On the other hand, the dynamics along the stable direction of the saddle (orthogonal to the one-dimensional manifold) are fast. Thus, after an initial transient taking the packet to the one-dimensional manifold, the amplitude and waveform of each packet remain relatively unchanged from transfer to transfer. Furthermore, the relative ordering of the amplitudes is retained even as the amplitude slowly changes.

Fixed points of the map $(m_{j}(jT),M^{1}_{j}(jT),M^{2}_{j}(jT))$ are shown in Fig. 4, showing a fold bifurcation (along the unstable manifold of the saddle) as a function of $S$.

As $S$ moves away from the region of graded transfer, the dynamics along the unstable manifold become fast. Thus, all amplitudes rapidly approach the fixed points at the extremes. This leads to a mechanism for binary information transfer (in distinction from graded information transfer), where pulse gating can only transfer low and high amplitudes. As $S$ increases, the location of the high fixed point increases giving more distinction between low and high fixed points.

This behavior may be understood across a large region of parameter space with the use of the approximate, Gaussian model (see Methods).

In Fig. 5A,B, we compare numerical and approximate (Gaussian) solutions. Note the similar advection of the population density with obvious differences in curvature and reset populations (which may be seen in Fig. 5A as a linear wall of density along the reset value, [math], of the membrane potential). In Fig. 5C, we show the fold bifurcation, found numerically in Fig. 4, for a range of values of $S$ and $\bar{I}^{\mathrm{g}}$. In Fig. 5D, we demonstrate the existence of a cusp catastrophe as a function of the width of the potential distribution $\sigma$. In Fig. 5E,F, we show how fast or slow (ghost) dynamics occur near the fold of the cusp depending on $\sigma$ and $S$. In Fig. 5E,F, blue shades indicate locations where feedforward amplitudes are slowly changing, green and yellow shades where amplitudes change rapidly. Isoclines (the same in both panels) are plotted denoting the fixed point (central isocline) and nearby locations of small, but non-zero change. Note how in Fig. 5E,F, dynamics are fast and amplitudes rapidly approach the attractors for large synaptic coupling, but dynamics remain slow in the region near the fold bifurcation. This region gives rise to an approximate line attractor that becomes more robust, i.e. covering a large range of amplitudes ($A/g_{L}$), as $\sigma$ increases (compare Fig. 5E and Fig. 5F).

IV Discussion

To understand pulse-gated information propagation, we used a Fokker-Planck analysis to derive input-output curves that may be used in population firing rate models of pulse-gated propagation and to reveal an approximate line attractor in the network of spiking neurons. We further showed that the line attractor is associated with ghost dynamics occurring along the front fold (relatively small synaptic coupling) of a cusp catastrophe as $\sigma$ varies in state space. Locally, the nearly one-dimensional attracting manifold is the result of a fold bifurcation where the stable directions of the saddle and attracting fixed points are strongly attracting and the unstable manifold exhibits ghost dynamics. In terms of input-output functions, this reflects that the fold of this cusp separates a region of sigmoidal f-I curves, with a relatively large linear interval, from a region of bistability. Since this region occupies a sizeable volume in parameter space, the graded propagation is robust. Furthermore, the fact that a common cusp catastrophe underlies the fast timescale dynamics of our system indicates that this type of line attractor is generic and will persist in feedforward networks of other types of spiking neurons.

One of the major problems of neuroscience is to understand how complex neural functions emerge from the collective dynamics of neuronal networks. Towards this goal, researchers have tried to construct models using mean-field firing rate theories and large-scale numerical simulations. However, except for a few examples (e.g. (Gerstner, 1995; Nykamp and Tranchina, 2000; Cai et al., 2006; Ostojic and Brunel, 2011; Montbrio et al., 2015; Fasoli et al., 2016)), the precise correspondence between the underlying microscopic spiking neurons and the macroscopic coherent dynamics of the neuronal populations has not been established.

Our analysis provides a significant step toward understanding how macroscopic attracting manifolds can emerge from the dynamic interactions of microscopic spiking neurons via the coordination of pulse-gating, synaptic weights, and intrinsic and extrinsic noise (i.e. the distribution of the membrane potential across the population), and offers possible order parameters for which macroscopic descriptions can be derived from the underlying microscopic dynamical model.

Furthermore, transfer mechanisms, such as those found in the graded and binary transfer parameter regimes shown in Fig. 5. (i.e. for $S$ near the fold (graded) or $S$ to the right of the fold (binary)), provide a novel means of understanding dynamic network interactions such as the detailed measurements of population activity underlying complex neural tasks provided by modern experimental techniques. Already, pulse-gated transfer mechanisms have been shown to be capable of implementing dynamic modules representing complex neural functions, such as short term memory, decision making, and the control of neural circuits (Sornborger et al., 2015; Wang et al., 2016).

V Appendix: Mean Field Model

A mean-field firing rate model of Eq. (1) is given by

[TABLE]

where $[]^{+}$ denotes a thresholded linear function, with threshold $g_{0}$, for the input-output relation of an I&F neuron.

In (Sornborger et al., 2015), we showed that when the gating pulse cancels the threshold ($\bar{I}^{g}=g_{0}$), and the feedforward synaptic coupling strength was

[TABLE]

we get exact, graded propagation, where the mean synaptic current and firing rates were

[TABLE]

and

[TABLE]

Acknowledgements.

L.T. thanks the UC Davis Mathematics Department for its hospitality. A.T.S. would like to thank Liping Wei and the Center for Bioinformatics at the School of Life Sciences at Peking University for their hospitality. This work was supported by the Natural Science Foundation of China grants 91232715 (Z.X, L.T.), 31771147 (L.T.), 91430216 (J.Z.), and U1530401 (J.Z.), by the Open Research Fund of the State Key Laboratory of Cognitive Neuroscience and Learning grant CNLZD1404 (Z.X., L.T.), by the Beijing Municipal Science and Technology Commission under contract Z151100000915070 (Z.X., L.T.), and by the Undergraduate Honors Research Program of the School of Life Sciences at Peking University (Z.X., L.T.).