Contrasting the effects of adaptation and synaptic filtering on the timescales of dynamics in recurrent networks

TL;DR

This study compares how adaptation and synaptic filtering influence the timescales of neural network dynamics, revealing that synaptic filtering prolongs network activity timescales while adaptation does not.

Contribution

It provides a detailed analysis showing that synaptic filtering and adaptation have distinct impacts on network timescales, with only synaptic filtering extending the overall activity duration.

Findings

Synaptic filtering increases network timescales proportionally.

Adaptation does not lead to larger network timescales.

Network dynamics are more influenced by synaptic filtering than adaptation.

Abstract

Neural activity exhibits a vast range of timescales that can be several fold larger than the membrane time constant of individual neurons. Two types of mechanisms have been proposed to explain this conundrum. One possibility is that large timescales are generated by a network mechanism based on positive feedback, but this hypothesis requires fine-tuning of the synaptic connections. A second possibility is that large timescales in the neural dynamics are inherited from large timescales of underlying biophysical processes, two prominent candidates being adaptive ionic currents and synaptic transmission. How the timescales of these processes influence the timescale of the network dynamics has however not been fully explored. To address this question, we analyze large networks of randomly connected excitatory and inhibitory units with additional degrees of freedom that correspond to adaptation or synaptic filtering. We determine the fixed points of the systems, their stability to perturbations and the corresponding dynamical timescales. Furthermore, we apply dynamical mean field theory to study the temporal statistics of the activity in the fluctuating regime, and examine how the adaptation and synaptic timescales transfer from individual units to the whole population. Our overarching finding is that synaptic filtering and adaptation in single neurons have very different effects at the network level. Unexpectedly, the macroscopic network dynamics do not inherit the large timescale present in adaptive currents. In contrast, the timescales of network activity increase proportionally to the time constant of the synaptic filter. Altogether, our study demonstrates that the timescales of different biophysical processes have different effects on the network level, so that the slow processes within individual neurons do not necessarily induce slow activity in large recurrent neural networks.

Tables

Table 1. Table 1: Parameter values. Parameter values used in the simulations.

Parameter	Value
Number of units $N$	3000
In-degree $C$	100
Excitatory inputs $C_E$	80
Inhibitory inputs $C_I$	20
Ratio I-E coupling strength $g$	4.1
Threshold $\gamma$	-0.5
Maximum firing rate ${\varphi }_{\text{max}}$	2

Equations

$$\varphi \left(x\right)=\{\begin{array}{ll}{\left[x-\gamma \right]}^{+}& \text{if }x-\gamma <{\varphi }_{\text{max}}\\ {\varphi }_{\text{max}}& \text{otherwise},\end{array}$$

$$\{\begin{array}{ll}& {\tau }_m{\dot{x}}_i\left(t\right)=-x_i\left(t\right)+{\sum }_{j=1}^N{J}_{ij}\varphi \left(x_j\left(t\right)\right)-g_w{w}_i\left(t\right)+I_i\left(t\right)\\ & {\tau }_w{\dot{w}}_i\left(t\right)=-w_i\left(t\right)+x_i\left(t\right)-\gamma ,\end{array}$$

$$\{\begin{array}{ll}& {\tau }_m{\dot{x}}_i\left(t\right)=-x_i\left(t\right)+s_i\left(t\right)\\ & {\tau }_s{\dot{s}}_i\left(t\right)=-s_i\left(t\right)+{\sum }_{j=1}^N{J}_{ij}\varphi \left(x_j\left(t\right)\right)+I_i\left(t\right).\end{array}$$

$$h_s\left(t\right)=\frac{1}{{\tau }_s-{\tau }_m}\left(e^{-\frac{t}{{\tau }_s}}-e^{-\frac{t}{{\tau }_m}}\right)\Theta \left(t\right),$$

$${\tau }^{\pm }=\frac{2{\tau }_m{\tau }_w}{{\tau }_w+{\tau }_m}{\left(1\pm \sqrt{1-\frac{4{\tau }_m{\tau }_w\left(1+g_w\right)}{{\left({\tau }_m+{\tau }_w\right)}^2}}\right)}^{-1}.$$

$$h_w\left(t\right)=\left(-\frac{g_w}{{\tau }_w}{e}^{-\left(1+g_w\right)\frac{t}{{\tau }_w}}+\frac{1}{{\tau }_m}{e}^{-\frac{t}{{\tau }_m}}\right)\Theta \left(t\right).$$

$$\left(1+g_w\right)x_0=J\left(C_E-g{C}_I\right)\varphi \left(x_0\right)+g_w\gamma ,$$

$$x_0=J\left(C_E-g{C}_I\right)\varphi \left(x_0\right),$$

$${\tau }^{+}$$

$${\tau }^{-}$$

$${\tau }^{+}$$

$${\tau }^{-}$$

$${\tau }_w>\frac{{\tau }_m}{g_w+\sqrt{2g_w\left(g_w+1\right)}},$$

$$\{\begin{array}{ll}& {\tau }_m{\dot{x}}_i\left(t\right)=-x_i\left(t\right)-g_w{w}_i\left(t\right)+I_i\left(t\right)\\ & {\tau }_w{\dot{w}}_i\left(t\right)=-w_i\left(t\right)+\varphi \left(x_i\left(t\right)\right),\end{array}$$

$$\varphi \left(x\right)=\{\begin{array}{ll}{\left[x-\gamma \right]}^{+}& \text{if }x-\gamma <{\varphi }_{\text{max}}\\ {\varphi }_{\text{max}}& \text{otherwise}.\end{array}$$

$$\{\begin{array}{ll}& {\tau }_m{\dot{x}}_i\left(t\right)=-x_i\left(t\right)-g_w{w}_i\left(t\right)+I_i\left(t\right)\\ & {\tau }_w{\dot{w}}_i\left(t\right)=-w_i\left(t\right)+x_i\left(t\right)-\gamma ,\end{array}$$

$$\{\begin{array}{ll}& {\tau }_m{\dot{x}}_i\left(t\right)=-x_i\left(t\right)+s_i\left(t\right)+I_i\left(t\right)\\ & {\tau }_s{\dot{s}}_i\left(t\right)=-s_i\left(t\right)+I_i\left(t\right).\end{array}$$

$$I_i\left(t\right)=I_{\text{syn},i}=\sum _j{J}_{ij}\varphi \left(x_j\left(t\right)\right).$$

$$J_{\text{eff}}:=\sum _j{J}_{ij}=J\left(C_E-g{C}_I\right).$$

$$x\left(t\right)=\left(h\ast I\right)\left(t\right)={\int }_{-\infty }^{+\infty }d{t}^{\prime }h\left(t^{\prime }\right)I\left(t-t^{\prime }\right).$$

$$z\left(t\right)={\int }_{-\infty }^{\infty }d{t}^{\prime }{e}^{A{t}^{\prime }}\Theta \left(t^{\prime }\right)b\left(t-t^{\prime }\right),$$

$$A_{s}=\left(\begin{array}[]{c c}-\tau_{m}^{-1}&\tau_{m}^{-1}\\ 0&-\tau_{s}^{-1}\end{array}\right).$$

$$h_s\left(t\right)=\Theta \left(t\right)\frac{1}{{\tau }_m-{\tau }_s}\left(e^{-\frac{t}{{\tau }_m}}-e^{-\frac{t}{{\tau }_s}}\right).$$

$$\displaystyle A_{w}=\left(\begin{array}[]{c c}-\tau_{m}^{-1}&-g_{w}\tau_{m}^{-1}\\ \tau_{w}^{-1}&-\tau_{w}^{-1}\end{array}\right).$$

$${\lambda }_w^{\pm }=\frac{1}{2}\left(-{\tau }_m^{-1}-{\tau }_w^{-1}\pm \sqrt{{\left({\tau }_m^{-1}+{\tau }_w^{-1}\right)}^2-4\left(1+g_w\right){\tau }_m^{-1}{\tau }_w^{-1}}\right).$$

$${\xi }^{\pm }={\left(\frac{g_w}{{\tau }_m},\frac{1}{2}\left(-\frac{1}{{\tau }_m}+\frac{1}{{\tau }_w}\mp \sqrt{{\left(\frac{1}{{\tau }_m}-\frac{1}{{\tau }_w}\right)}^2-4\frac{g_w}{{\tau }_m{\tau }_w}}\right)\right)}^T.$$

$${\lambda }_w^{+}=$$

$${\lambda }_w^{-}=$$

$$h_w\left(t\right)=\frac{1}{{\xi }_1^{+}{\xi }_2^{-}-{\xi }_1^{-}{\xi }_2^{+}}\left({\xi }_1^{+}{\xi }_2^{-}{e}^{{\lambda }^{+}t}-{\xi }_1^{-}{\xi }_2^{+}{e}^{{\lambda }^{-}t}\right),$$

$${\xi }_{-}$$

$${\xi }_{+}$$

$$h_w\left(t\right)=\frac{g_w}{{\tau }_m\left(2g_w+1\right)-{\tau }_w}{e}^{-\frac{1+g_w}{{\tau }_w}t}+\frac{1}{{\tau }_m}\frac{1-\left(1+g_w\right)\frac{{\tau }_m}{{\tau }_w}}{1-\left(1+2g_w\right)\frac{{\tau }_m}{{\tau }_w}}{e}^{-\left(\frac{1}{{\tau }_m}-\frac{g_w}{{\tau }_w}\right)t}.$$

Full text

Contrasting the effects of adaptation and synaptic filtering on the timescales of dynamics in recurrent networks

Manuel Beiran1*, Srdjan Ostojic1*

1 Group for Neural Theory, Laboratoire de Neurosciences Cognitives Computationnelles, Département d’Études Cognitives, École Normale Supérieure, INSERM U960, PSL University, Paris, France

• manuel.beiran@ens.fr (MB); srdjan.ostojic@ens.fr (SO)

Abstract

Neural activity in awake behaving animals exhibits a vast range of timescales that can be several fold larger than the membrane time constant of individual neurons. Two types of mechanisms have been proposed to explain this conundrum. One possibility is that large timescales are generated by a network mechanism based on positive feedback, but this hypothesis requires fine-tuning of the strength or structure of the synaptic connections. A second possibility is that large timescales in the neural dynamics are inherited from large timescales of underlying biophysical processes, two prominent candidates being intrinsic adaptive ionic currents and synaptic transmission. How the timescales of adaptation or synaptic transmission influence the timescale of the network dynamics has however not been fully explored. To address this question, here we analyze large networks of randomly connected excitatory and inhibitory units with additional degrees of freedom that correspond to adaptation or synaptic filtering. We determine the fixed points of the systems, their stability to perturbations and the corresponding dynamical timescales. Furthermore, we apply dynamical mean field theory to study the temporal statistics of the activity in the fluctuating regime, and examine how the adaptation and synaptic timescales transfer from individual units to the whole population. Our overarching finding is that synaptic filtering and adaptation in single neurons have very different effects at the network level. Unexpectedly, the macroscopic network dynamics do not inherit the large timescale present in adaptive currents. In contrast, the timescales of network activity increase proportionally to the time constant of the synaptic filter. Altogether, our study demonstrates that the timescales of different biophysical processes have different effects on the network level, so that the slow processes within individual neurons do not necessarily induce slow activity in large recurrent neural networks.

Introduction

Adaptive behaviour requires processing information over a vast span of timescales [1], ranging from micro-seconds for acoustic localisation [2], milliseconds for detecting changes in the visual field [3], seconds for evidence integration [4] and working memory [5], to hours, days or years in the case of long-term memory. Neural activity in the brain is matched to the computational requirements imposed by behaviour, and consequently displays dynamics over a similarly vast range of timescales [6, 7, 8]. Since the membrane time constant of an isolated neuron is of the order of tens of milliseconds, the origin of the long timescales observed in the neural activity has been an outstanding puzzle.

Two broad classes of mechanisms have been proposed to account for the existence of long timescales in the neural activity. The first class relies on non-linear collective dynamics that emerge from synaptic interactions between neurons in the local network. Such mechanisms have been proposed to model a variety of phenomena that include working memory [9], decision-making [10] and slow variability in the cortex [11]. In those models, long timescales emerge close to bifurcations between different types of dynamical states, and therefore typically rely on the fine tuning of some parameter [12]. An alternative class of mechanisms posits that long timescales are directly inherited from long time constants that exist within individual neurons, at the level of hidden internal states [13]. Indeed biophysical processes at the cellular and synaptic level display a rich repertoire of timescales. These include short-term plasticity that functions at the range of hundreds of milliseconds [14, 15], a variety of synaptic channels with timescales from tens to hundreds of milliseconds [16, 17, 18, 19], ion channel kinetics implementing adaptive phenomena [20], calcium dynamics [21] or shifts in ionic reversal potentials [22]. How the timescales of these internal processes affect the timescales of activity at the network level has however not been fully explored.

In this study, we focus on adaptative ion-channel currents, which are known to exhibit timescales over several orders of magnitude [23, 24, 25]. We contrast their effects on recurrent network dynamics with the effect of the temporal filtering of inputs through synaptic currents, which also expands over a large range of timescales [26]. To this end, we extend classical rate models [27, 28, 29, 30] of randomly connected recurrent networks by including for each individual unit a hidden variable that corresponds to either the adapting of the synaptic current. We systematically determine the types of collective activity that emerge in such networks. We then compare the timescales on the level of individual units with the activity within the network.

Results

We consider $N$ coupled inhibitory and excitatory units whose dynamics are given by two variables: the input current $x_{i}$ and a slow variable $s_{i}$ or $w_{i}$ that accounts for the synaptic filtering or adaptation current respectively. The instantaneous firing rate of each neuron is obtained by applying a static non-linearity $\phi\left(x\right)$ to the input current at every point in time. For simplicity, we use a positive and bounded threshold-linear transfer function

[TABLE]

where $[\cdot]^{+}$ indicates the positive part, $\gamma$ is the activation threshold and $\phi_{\text{max}}$ the maximum firing rate.

Single neuron adaptation is described by the variable $w\left(t\right)$ that low-pass filters the linearized firing rate with a timescale $\tau_{w}$, slower than the membrane time constant $\tau_{m}$, and feeds it back with opposite sign into the input current dynamics (see Methods). The dynamics of the $i$-th adaptive neuron are given by

[TABLE]

where $I_{i}\left(t\right)$ is the external input current to neuron $i$.

Synaptic filtering consists in low-pass filtering the synaptic input received by a cell with time constant $\tau_{s}$, before it contributes to the input current. The dynamics of the $i$-th neuron in a network with synaptic filtering are

[TABLE]

The matrix element $J_{ij}$ corresponds to the synaptic coupling strength from neuron $j$ onto neuron $j$. In this study we focus on neuronal populations of inhibitory and excitatory units, whose connectivity is sparse, random, with constant in-degree: all neurons receive exactly the same number of excitatory and inhibitory connections, $C_{E}$ and $C_{I}$, as in [31, 32, 33]. All excitatory synapses have equal strength $J$ and all inhibitory neurons $-gJ$. Furthermore, we consider the large network limit where the number of synaptic neurons $N$ is large while keeping the excitatory and inhibitory inputs $C_{E}$ and $C_{I}$ fixed.

Single unit: timescales of dynamics

In the models studied here the input current of individual neurons is described by a linear system. Thus, their activity is fully characterized by the response $h\left(t\right)$ to a brief impulse signal, i.e. the linear filter. When such neurons are stimulated with a time-varying input $I\left(t\right)$, the response is the convolution of the filter with the input, $x\left(t\right)=\left(h\ast I\right)\left(t\right)$. These filters can be determined analytically for both neurons with adaptation or synaptic filtering and directly depend on the parameters of these processes. Analyzing the differences that these two slow processes produce in the linear filters is useful for studying the differences in the response of adaptive and synaptic filtering neurons to temporal stimuli (Fig 1 A), and will serve as a reference for comparison to the effects that emerge at the network level.

In particular, the filter of a neuron with synaptic filtering, $h_{s}\left(t\right)$, is the sum of two exponentially decaying filters of opposite signs and equal amplitude, with time constants $\tau_{s}$ and $\tau_{m}$:

[TABLE]

where $\Theta\left(t\right)$ is the Heaviside function (see Methods). Thus, the current response of a neuron to an input pulse received from an excitatory presynaptic neuron is positive and determined by two different timescales. The response first grows with timescale $\tau_{m}$, so that the neuron cannot respond to any abrupt changes in the synaptic input faster than this timescale, and then decreases back to zero with timescale $\tau_{s}$ (grey curves, Fig 1 B).

The adaptation filter is given as well by the linear combination of two exponential functions. In contrast to the synaptic filter, since the input in the adaptive neuron model affects directly the current variable $x_{i}\left(t\right)$, there is an instantaneous change in the firing rate to an input delta-function (red curves, Fig 1 B). The timescales of the two exponentials can be calculated as

[TABLE]

When the argument of the square root in Eq (5) is negative, i.e. for an adaptation time constant $\tau_{w}$ smaller than $\tau_{m}\left(4\left(1+g_{w}\right)\right)^{-1}$, the two timescales correspond to a pair of complex conjugate numbers, so that the filter is an oscillatory function whose amplitude decreases monotonically to zero at a single timescale. If the argument of the square root is positive, for slow enough adaptation, the two timescales are real numbers and correspond to exponential functions of opposing signs of decaying amplitude. However, the amplitudes of these two exponentials are different (see Methods). To illustrate this, we focus on the limit of large adaptation time constants with respect to the membrane time constant, where the two exponential functions evolve with timescales that decouple the contribution of the membrane time constant and the adaptation current. In that limit, the adaptive filter reads

[TABLE]

The amplitude of the slow exponential is inversely related to its timescale so that the integral of this mode is fixed, and independent of the adaptation time constant. This implies that a severalfold increase of the adaptation time constant does not lead to strong changes in the single neuron activity for time-varying signals (Fig 1A).

Furthermore, we can characterize the timescale of the single neuron response as the sum of the exponential decay timescales weighed by their relative amplitude, and study how this characteristic timescale evolves as a function of the time constants of either the synaptic or the adaptive current (Fig 1C). For adaptive neurons, the activity timescale is bounded as a consequence of the decreasing amplitude of the slow mode, i.e. increasing the adaptation time constant beyond a certain value will not lead to a slower response. In contrast, the activity of an individual neuron with synaptic filtering scales proportionally to the synaptic filter time, since the relative amplitudes of the two decaying exponentials are independent of the time constants.

When any of the two neuron types are stimulated with white Gaussian noise, the variance in the response is always smaller than the input variance, due to the low pass filtering properties of the neurons. However, this gain in the variance of the input currents is modulated by the different neuron parameters (Fig 1D). For a neuron with synaptic filtering, the gain is inversely proportional to the time constant $\tau_{s}$. In contrast, for a neuron with adaptation, increasing the adaptation time constant has the opposite effect of increasing the variance of the current response. This is because when the adaptation time constant increases, the amplitude of the slow exponential decreases accordingly, and the low-pass filtering produced by this slow component is weaker. Following the same reasoning, increasing the adaptation coupling corresponds to strengthening the low-pass filtering performed by adaptation, so that the variance decreases (Fig 1D, dashed vs full red curves).

Population-averaged dynamics

In the absence of any external input, a non-trivial equilibrium for the population averaged activity emerges due to the recurrent connectivity of the network. The equilibrium firing rate is identical across network units, since all units are statistically equivalent. We can write the input current $x_{0}$ at the fixed point as the solution to the transcendental equation

[TABLE]

for the network with adaptation, and to

[TABLE]

for synaptic filtering (see Methods). Based on Eq (7), we find that the adaptation coupling $g_{w}$ reduces the mean firing rate of the network, independently of whether the network is dominated by inhibition or excitation (Fig 2A). Synaptic filtering instead does not play any role in determining the equilibrium activity of the neurons, since Eq (8) is independent of the synaptic filtering parameter $\tau_{s}$.

We next study the stability and dynamics of the equilibrium firing rate in response to a small perturbation uniform across the network, $x_{i}\left(t\right)=x_{0}+\delta x\left(t\right)$. Because of the fixed in-degree of the connectivity matrix, the linearized dynamics of each neuron are identical, so that the analysis of the homogeneous perturbation on the network reduces to the study of a two-dimensional deterministic system of differential equations which corresponds to the dynamics of the population-averaged response (see Methods). The stability and timescales around equilibrium depend on the two eigenvalues of this linear 2D-system. More specifically, the fixed point is stable to a homogeneous perturbation if the two eigenvalues of the dynamic system have negative real part, in which case the inverse of the unsigned real part of the eigenvalues determines the timescales of the response. For both the network with synaptic filtering and the network with adaptive neurons, the order parameter of the connectivity that determines the stability of the fixed point is the effective recurrent coupling $J\left(C_{E}-gC_{I}\right)$ each neuron receives, resulting from the sum of all input synaptic connections. A positive (negative) effective coupling corresponds to a network where recurrent excitation (inhibition) dominates and the recurrent input provides positive (negative) feedback [32, 33].

For networks with synaptic filtering, we find that the synaptic time constant does not alter the stability of the equilibrium state, so that the effective coupling alone determines the stability of the population-averaged activity. As the effective input coupling strength is increased, the system undergoes a saddle-node bifurcation when the effective input is $J\left(C_{E}-gC_{I}\right)=1$ (Fig 2C). In other words, the strong positive feedback loop generated by the excitatory recurrent connections destabilizes the system.

To analyze the timescales elicited by homogeneous perturbations, we calculate the eigenvalues and eigenvectors of the linearized dynamic system (see Methods). We find that for inhibition-dominated networks ($J\left(C_{E}-gC_{I}\right)<0$), the network shows population-averaged activity at timescales that interpolate between the membrane time constant and the synaptic time constant. As the effective coupling is increased, the slow timescale at the network level can be made arbitrarily slow by tuning the effective synaptic coupling close to the bifurcation value, a well-known network mechanism to achieve slow neural activity [12].

In the limit of very slow synaptic timescale, the two timescales of the population-averaged activity are

[TABLE]

so that the timescale $\tau^{-}$ is proportional to the membrane time constant and $\tau^{+}$ is proportional to the slow synaptic time constant, effectively decoupling the two timescales. The relative contribution of these two timescales is the same, independently of the time constant $\tau_{s}$, as we found in the single neuron analysis.

The network with adaptation shows different effects on the population-averaged activity. First, the presence of adaptation modifies the region of stability: the system is stable when the effective recurrent input $J\left(C_{E}-gC_{I}\right)$ is less than the minimum of $1+g_{w}$ and $1+\frac{\tau_{m}}{\tau_{w}}$ (see Methods). Therefore, the stability region is larger than for the network with synaptic filtering (Fig 2B vs Fig 2C). In other words, the effective excitatory feedback required to destabilize the network is larger due to the counterbalance provided by adaptation. Moreover, adaptation allows the network to undergo two different types of bifurcations as the effective input strength increases, depending on the adaptation parameters. One possibility is a saddle-node bifurcation, as in the synaptic case, which takes place when $J\left(C_{E}-gC_{I}\right)=1+g_{w}$. Beyond the instability all neurons in the network saturate. The other possible bifurcation, which happens if $\frac{\tau_{m}}{\tau_{w}}<g_{w}$, at an effective coupling strength $J\left(C_{E}-gC_{I}\right)=1+\frac{\tau_{m}}{\tau_{w}}$, is a Hopf bifurcation: the fixed point of network becomes unstable, leading in general to oscillating dynamics of the population-averaged response. Note that in the limit of very slow adaptation, the system can only undergo a Hopf bifurcation (Fig 2B).

The two timescales of the population-averaged activity in the stable regime for the adaptive network decouple the two single neuron time constants when adaptation is much slower than the membrane time constant. In this limit, up to first order of the adaptive time ratio $\frac{\tau_{m}}{\tau_{w}}$, the two activity timescales are

[TABLE]

Similar to the single neuron dynamics, the amplitude of the slow mode, corresponding to $\tau^{-}$, decreases as $\tau_{w}$ is increased, so that the contribution of the slow timescale is effectively reduced when $\tau_{w}$ is very large. On the contrary, the mode corresponding to $\tau^{+}$, proportional to the membrane time constant can be tuned to reach arbitrarily large values. This network mechanism to obtain slow dynamics does not depend on the adaptation properties.

Heterogeneous activity

Linear stability analysis

Previous studies have shown that random connectivity can lead to heterogeneous dynamics where the activity of each unit fluctuates strongly in time [29, 34, 35, 33]. To assess the effects of additional hidden degrees of freedom on the emergence and timescales of such fluctuating activity, we examine the dynamics when each unit is perturbed independently away from the equilibrium, $x_{i}\left(t\right)=x_{0}+\delta x_{i}\left(t\right)$. By linearizing the full $2N$-dimensional dynamics around the fixed point, we can study the stability and timescales of the activity characterized by the set of eigenvalues of the linearized system, $\lambda_{s}$ and $\lambda_{w}$ for the network with synaptic filtering neurons and adaptation, respectively. These sets of eigenvalues are determined by a direct mapping to the eigenvalues of the connectivity matrix, $\lambda_{J}$ (see Methods). The eigenvalues $\lambda_{J}$ of the connectivity matrices considered are known in the limit of large networks [36, 33]: they are enclosed in a circle of radius $J\sqrt{C_{E}+g^{2}C_{I}}$, except for an outlier that corresponds to the population-averaged dynamics, studied in the previous section. Therefore, we can map the circle that encloses the eigenspectrum $\lambda_{J}$ into a different shape in the space of eigenvalues $\lambda_{s/w}$ (insets Fig. 3). In order to determine the stability of the response to the perturbation, we assess whether the real part of the eigenspectrum $\lambda_{s/w}$ is negative at all possible points. Furthermore, the type of bifurcation is determined by whether the curve enclosing the eigenvalues $\lambda_{s,w}$ crosses the imaginary axis at zero frequency or at a finite frequency when the synaptic coupling strength is increased, leading respectively to a zero-frequency or to a Hopf bifurcation [37].

The order parameter of the connectivity that affects the stability and dynamics of the network is now the radius of the circle of eigenvalues $\lambda_{J}$, i.e. $J\sqrt{\left(C_{E}+g^{2}C_{I}\right)}$. This parameter is the standard deviation of the synaptic input weights of a neuron (see Methods), which contrasts with the order parameter of the population-averaged response, that depends on the mean of the synaptic input weights. The mean and standard deviation of the synaptic connectivity can be chosen independently, so that while the population-averaged activity remains stable, the individual neurons might not display stable dynamics. To analyze solely the heterogeneous response of the network to the perturbation, we focus in the following on network connectivities whose population-averaged activity is stable, i.e. the effective synaptic coupling is inhibitory or weakly excitatory.

We find that in the network with synaptic filtering, the eigenspectrum $\lambda_{s}$ always crosses the stability bound through the real axis, which takes place when the spectral radius of the connectivity is one, $J\sqrt{C_{E}+g^{2}C_{I}}=1$. Thus the system undergoes a zero-frequency bifurcation similar to randomly connected networks without hidden variables [29, 35, 38, 33], leading to strong fluctuations at the single neuron level that are self-sustained by the network connectivity (Fig 3 Bii-Biv). The critical coupling at which the equilibrium firing rate loses stability is independent of the synaptic time constant, i.e. synaptic filtering does not affect the stability of heterogeneous responses (Fig 4 A). However, the synaptic time constant $\tau_{s}$ affects the timescales at which the system returns to equilibrium after a perturbation, because the eigenvalues $\lambda_{s}$ (see Eq (73) in Methods) depend explicitly on $\tau_{s}$.

For a network with adaptive neurons, we calculate the eigenspectrum $\lambda_{w}$ and find that the transition to instability $\text{Re}\left(\lambda_{w}\right)=0$ can happen either at zero frequency or at a finite frequency (see Methods), leading to a Hopf bifurcation (as in inset Fig 3 Aiii). In particular, the network dynamics undergo a Hopf bifurcation when

[TABLE]

so that strong adaptation coupling and slow adaptation time constants lead to a finite frequency bifurcation. In particular, if the coupling $g_{w}$ is larger than $\sqrt{5}-2\approx 0.236$, only the Hopf bifurcation is possible, since by construction $\frac{\tau_{m}}{\tau_{w}}<1$. We can also calculate the frequency of oscillations at the Hopf bifurcation. We find that, for slow adaptive currents, the Hopf frequency is inversely related to the adaptation time constant (Fig 4B), so that slower adaptation currents produce slower oscillations at the bifurcation.

Adaptation also increases the stability of the equilibrium firing rate to a heterogeneous perturbation, in comparison to a network with synaptic filtering (Fig 4 C). This can be intuitively explained in geometrical terms by analyzing how adaptation modifies the shape of the eigenspectrum $\lambda_{w}$ with respect to the circular eigenspectrum of the connectivity matrix $\lambda_{J}$.

The Hopf bifurcation leads to the emergence of a new dynamical regime in the network (Fig 3 Aiv), which is studied in the following section. Right at the Hopf bifurcation, the system shows marginal oscillations at a single frequency that can be reproduced in finite-size simulations whenever only one pair of complex conjugate eigenvalues is unstable (Fig 3 Aiii).

Fluctuating activity: dynamical mean field theory

The classical tools of linear stability theory applied so far can only describe the dynamics of the system up to the bifurcation. To study the fluctuating regime, we take a different approach and focus on the temporal statistics of the activity, averaged over different connectivity matrices: we determine the mean and autocorrelation function of the single neuron firing rate, and characterize the timescale of the fluctuating dynamics [29, 34, 39, 40, 35, 38, 33]. For large networks, the dynamics can be statistically described by applying dynamical mean field theory (DMFT), which approximates the deterministic input to each unit by an independent Gaussian noise process. The full network is then reduced to a two-dimensional stochastic differential equation, where the first and second moments of the noise must be calculated self-consistently. We solve the self-consistent equations using a numerical iterative procedure, similar to the schemes followed in [41, 42, 43, 44, 34] (see Methods for an explanation of the iterative algorithm and its practical limitations).

For the network with synaptic filtering, we find that the autocorrelation function of the firing rates in the fluctuating regime corresponds to a monotonically decreasing function (Fig 5 A), qualitatively similar to the correlation obtained in absence of synaptic filtering [33]. This fluctuating state has often been referred to as rate chaos and shows non-periodical heterogeneous activity which is intrinsically generated by the network connectivity. The main effect of synaptic filtering is on the timescale of these fluctuations. When the synaptic time constant is much larger than the membrane time constant, the timescale of the network activity is proportional to the synaptic time constant $\tau_{s}$, as indicated by the linear dependence between the half-width of the autocorrelation function and the synaptic timescale $\tau_{s}$, when all other network parameters are fixed (Fig 5 B).

For the network with adaptation, we focus on large adaptation time constant $\tau_{w}$, where the network dynamics always undergo a Hopf bifurcation. The autocorrelation function in such a case displays damped oscillations (Fig 5 C). The decay in the envelope of the autocorrelation function is due to the chaotic-like fluctuations of the firing rate activity.

We define the time lag at which the envelope of the autocorrelation function decreases as the timescale of the network dynamics (see Methods). The timescale of the activity increases as the adaptation timescale is increased, when all the other parameters are fixed (Fig 5 D). However, this activity timescale saturates for large values of the adaptation timescale: the presence of very slow adaptive currents, beyond a certain value, will not slow down strongly the network activity. This saturation value depends on the connectivity strength.

Effects of noise

The networks studied so far, for a fixed connectivity matrix, are completely deterministic. We next study the effects of additional white noise inputs to each neuron, as a proxy towards understanding recurrent networks of spiking neurons with adaptation and synaptic filtering. On the mean-field level, such noise is equivalent to studying a recurrent network whose neurons fire action potentials as a Poisson process with instantaneous firing rate $\phi\left(x_{i}\left(t\right)\right)$ [45, 35].

Numerical simulations show that in the stable regime the additive external noise generates weak, fast stationary dynamics around the fixed point (Fig 6 Ai, Bi). The timescale of these fluctuations and their amplitude depend on the distance of the eigenspectrum to the stability line, so that the stable fluctuations for weak synaptic coupling strength (Fig 6Ai) are smaller in amplitude than those for larger coupling strength (Fig 6Aii), whose eigenspectrum is closer to the stability boundary. For adaptation, in the fluctuating regime beyond the Hopf bifurcation, the network activity shows again a combination of fluctuating activity and oscillations.

We further extend the DMFT analysis to account for the additional variance of the external white noise sources (see Methods). The autocorrelation function of the firing rates, as predicted by DMFT, does not vary drastically when weak noise is added to the network, except for very short time lags, at which white noise introduces fast fluctuations (see Fig 7). For the network with adaptation, the autocorrelation function of the firing rates still shows damped oscillations (Fig 7 A), while for the network with synaptic filtering, similarly, weak noise does not affect much the decay of the autocorrelation function (Fig 7 D). Very strong external noise on the other hand will reduce the effect of the underlying recurrent dynamics of the rate network, since the signal to noise ratio in the synaptic input of all neurons is low.

For a fixed external noise intensity, reducing the adaptation coupling or increasing the adaptation time constant increases the variance of the firing rate (Fig 7B), which resembles the dependence of the variance gain for individual neurons (Fig 1D). Conversely, slower synaptic filtering reduces the variance of the neuron’s firing rates. This is because in the network with synaptic filtering the noise is also filtered at the synapses –in the limit of very large $\tau_{s}$, the whole white noise is filtered out– whereas in the network with adaptation the noise affects directly the input current, without being first processed by the adaptation variable.

However, the timescale of the activity is nonetheless drastically affected by strong noise. External noise adds fast fluctuations on top of the intrinsically generated dynamics of the heterogeneous network with adaptation or synaptic filtering. If the noise is too strong, the effective timescale of the activity takes into account mostly this fast component. In that limit, the timescale of the activity is almost independent of the synaptic or adaptive time constants (Fig 7 C and F, largest noise intensity).

Discussion

We examined dynamics of excitatory-inhibitory networks in which each unit had a hidden degree of freedom that represented either firing-rate adaptation or synaptic filtering. The core difference between adaptation and synaptic filtering was how external inputs reached the single-unit activation variable that represents the membrane potential. In the case of adaptation, the inputs directly entered the activation variable, which was then filtered by the hidden, adaptive variable through a negative feedback loop. In the case of synaptic filtering, the external inputs instead reached first the hidden, synaptic variable and were therefore low-pass filtered before being propagated in a feed-forward fashion to the activation variable. While both mechanisms introduce a second timescale in addition to the membrane time constant, our main finding is that the interplay between those two timescales is very different in the two situations. Surprisingly, in presence of adaptation, the membrane timescale remains the dominant one in the dynamics, while the contribution of the adaptation timescale appears to be weak. In contrast, in a network with synaptic filtering, the dominant timescale of the dynamics is directly set by the synaptic variable, and the overall dynamics are essentially equivalent to a network in which the membrane time-constant is replaced with the synaptic one.

We used a highly abstracted model, in which each neuron is represented by membrane current that is directly transformed into a firing-rate through a non-linear transfer function. This class of models has been popular for dissecting dynamics in excitatory-inhibitory [27, 28, 46, 47, 48] or randomly-connected networks [29, 30, 33], and for implementing computations [49, 50]. Effects of adaptation in this framework have to our knowledge not been examined so far. We therefore extended the standard rate networks by introducing adaptation in an equally abstract fashion[24], as a hidden variable specified solely by a time constant and a coupling strength. Different values of those parameters can be interpreted as corresponding to different specific membrane conductances that implement adaptation, e.g. the calcium dependent potassium $I_{ahp}$ current or the slow voltage-dependent potassium current $I_{m}$, which are known to exhibit timescales over several orders of magnitude [51, 52]. To cover the large range of adaptation timescales observed in experiments [23], it would be straightforward to superpose several hidden variables with different time constants. Our approach could also be easily extended to include simultaneously adaptation and synaptic filtering.

A number of previous works have studied the effects of adaptation within more biologically constrained, integrate-and-fire models. These works have in particular examined the effects of adaptation on the spiking statistics [53, 54, 55], firing-rate response [56, 57], synchronisation [25, 58, 59, 55, 60], perceptual bistability [61] or single-neuron coding [62, 63]. In contrast, we have focused here on the relation between the timescales of adaptation and those of network dynamics. While our results rely on a simplified firing-rate model, we expect that they can be directly related to networks of spiking neurons by exploiting quantitative techniques for mapping adaptive integrate-and-fire models to effective firing rate descriptions[64].

A side result of our analysis is the finding that strong coupling in random recurrent networks with adaptation generically leads to a novel dynamical state, in which individual units exhibit a mixture of oscillatory and strong temporal fluctuations. The characteristic signature of this dynamical state is a damped oscillation found in the auto-correlation function of single-unit activity. In contrast, classical randomly connected networks lead to a fluctuating, chaotic state in which the auto-correlation function decays monotonically [29, 34, 35, 33]. Note that the oscillatory activity of different units is totally out of phase, so that no oscillation is seen at the level of population activity. This dynamical phenomenon is analogous to heterogeneous oscillations in anti-symmetrically connected networks with delays [37]. In both cases, the oscillatory dynamics emerge through a bifurcation in which a continuum of eigenvalues crosses the instability line at a finite-frequency. Similar dynamics can be also found in networks in which the connectivity is a superposition of a random and a rank two structured part [33]. In that situation, the heterogeneous oscillations however originate from a Hopf bifurcation due to an isolated pair of eigenvalues that correspond to the structured part of the connectivity.

Our main aim here was to determine how hidden variables could induce long timescales in randomly-connected networks. Long timescales could alternatively emerge from non-random connectivity structure. As extensively investigated in earlier works, one general class of mechanism relies on setting the connectivity parameters close to a bifurcation that induces arbitrarily long timescales [29, 12]. Another possibility is that non-random features of the connectivity, such as the over-representation of reciprocal connections [65, 66] slow down the dynamics away from any bifurcation. A recent study [67] has indeed found such a slowing-down. Weak connectivity structure of low-rank type provides yet another mechanism for the emergence of long timescales. Indeed, rank-two networks can generate slow manifolds corresponding to ring attractors provided a weak amount of symmetry is present [68].

Ultimately, the main reason for looking for long timescales in the dynamics is their potential role in computations performed by recurrent networks [69, 70]. Recent works have proposed that adaptive currents may help implement computations in spiking networks by either introducing slow timescales or reducing the amount of noise due to spiking [71, 72]. Our results suggest that synaptic filtering is a much more efficient mechanism to this end than adaptation. Identifying a clear computational role for adaptation in recurrent networks therefore remains an open and puzzling question.

Methods

Network model

We compare the dynamics of two different models: a recurrent network with adaptive neurons, and a recurrent network with synaptic filtering. Each model is defined as a set of $2N$ coupled differential equations. The state of the $i$-th neuron is determined by two different variables, the input current $x_{i}\left(t\right)$ and the adaptation (synaptic) variable $w_{i}\left(t\right)$ ($s_{i}\left(t\right)$).

Adaptation

The dynamics of the recurrent network with adaptive neurons are given by

[TABLE]

where $\phi\left(x\right)$ is a monotonically increasing non-linear function that transforms the input current into firing rate. In this study, we use a threshold-linear transfer function with saturation:

[TABLE]

In Eq (14) adaptation in single neuron rate models is defined as a low-pass filtered version with timescale $\tau_{w}$ of the neuron’s firing rate $\phi\left(x_{i}\left(t\right)\right)$, and is fed back negatively into the input current, with a strength that we call the adaptation coupling $g_{w}$. For the sake of mathematical tractability, we linearize the dynamics of the adaptation variable by linearizing the transfer function (Eq 15), $\phi\left(x_{i}\left(t\right)\right)\approx x_{i}\left(t\right)-\gamma$. Therefore, the dynamics of the network model with adaptation studied here read

[TABLE]

Note that this approximation allows for adaptation to increase the input current of a neuron, when the neuron’s current is below the activation threshold $\gamma$.

Synaptic filtering

For the recurrent network with synaptic filtering, the dynamics are

[TABLE]

In Eqs (14), (16), and (17), $I\left(t\right)$ represents the total external input received by the neuron. In general, we are interested in the internally generated dynamical regimes of the network, so that the input is given by the synaptic inputs

[TABLE]

The matrix element $J_{ij}$ indicates the coupling strength of the $j$-th neuron onto the $i$-th neuron. The connectivity matrix is sparse and random, with constant in-degree [32, 73, 33]: all neurons receive the same number of input connections $C$, from which $C_{E}$ are excitatory and $C_{I}$ inhibitory. All excitatory synapses have coupling strength $J$ while the strength of all inhibitory synapses is $-gJ$. Moreover, each neuron can only either excite or inhibit the rest of the units in the network, following Dale’s principle. Therefore, the total effective input coupling strength, which is the same for all neurons, is

[TABLE]

Single neuron dynamics

The dynamics of each individual neuron are described by a two-dimensional linear system, which implies that the input current response $x\left(t\right)$ to a time-dependent input $I\left(t\right)$ is the convolution of the input with a linear filter $h\left(\tau\right)$ that depends on the parameters of the linear system:

[TABLE]

In general, for any linear dynamic system $\dot{z}\left(t\right)=Az+b\left(t\right)$, where $A$ is a square matrix in $\mathbb{R}^{N\times N}$ and $b\left(t\right)$ is a $N$-dimensional vector, the dynamics are given by

[TABLE]

where $\Theta\left(t\right)$ is the Heaviside function. Thus, comparing Eqs (21) and (20), the linear filter is determined by the elements of the so-called propagator matrix $P\left(t\right)=e^{At}\Theta\left(t\right)$.

Synaptic filtering

For a single neuron wit synaptic filtering, the dynamics are given by Equation (17), where the input $I_{i}\left(t\right)$ represents the external current. We write the response in its vector form $\left(x\left(t\right),s\left(t\right)\right)^{T}$ and the input as $\left(0,I\left(t\right)\right)^{T}$. The dynamic matrix is

[TABLE]

The linear filter, $h_{s}\left(t^{\prime}\right)$, is given by the entries of the propagator matrix that links the input $I\left(t\right)$ to the output element $x\left(t\right)$, which are in this case only the entry in row one and column two: $h_{s}\left(t^{\prime}\right)=\left[P\left(t^{\prime}\right)\right]_{12}$. To compute the required entry of the propagator, we diagonalize the dynamic matrix $A=VDV^{-1}$. The matrix $D$ is a diagonal matrix with the eigenvalues of matrix $A$ in the diagonal entries, and $V$ is a matrix whose columns are the corresponding eigenvectors. Applying the identity $e^{tVDV^{-1}}=Ve^{tD}V^{-1}$ and the definition of propagator we obtain that

[TABLE]

The two timescales of the activity are defined by the inverse of the eigenvalues of the system, which coincide with $\tau_{m}$ and $\tau_{s}$. Every time a pulse is given to the neuron, both modes get activated with equal amplitude and opposing signs, as indicated by Eq 23. This means that there is a fast ascending phase after a pulse, at a temporal scale $\tau_{m}$, and a decay towards zero with timescale $\tau_{s}$.

Adaptation

The dynamics of a single adaptive neuron are determined by Equation (16), where $I_{i}\left(t\right)$ is the external input to the neuron. We apply the same procedure to determine the timescales of the response of an adaptive neuron to time-dependent perturbations. The dynamic matrix for an adaptive neuron reads

[TABLE]

Its eigenvalues are

[TABLE]

and the eigenvectors

[TABLE]

The eigenvalues are complex if and only if $g_{w}>\left(4\tau_{m}\tau_{w}\right)^{-1}\left(\tau_{w}-\tau_{m}\right)^{2}$, and in that case their real part is $\frac{1}{2\tau_{m}\tau_{w}}\left(\tau_{m}+\tau_{w}\right)$. As the adaptive time constant becomes slower, at a certain critical adaptation time constant both eigenvalues become real. We are interested in the behavior when the adaptation time constant is large. The absolute value of the inverse of the eigenvalues determines the time constants of the dynamics. Therefore, for large $\tau_{w}$ we can calculate the two real eigenvalues to first order of $\tau_{w}^{-1}$

[TABLE]

In this limit of slow adaptation, the time constant of one eigenmode is proportional to $\tau_{w}$, whereas the second mode scales with $\tau_{m}$. We are interested in the amplitude of each mode with respect to the other.

By explicitly calculating the first entry of the propagator matrix we obtain the adaptive filter in terms of the eigenvectors and eigenvalues,

[TABLE]

where we use the notation $\xi_{1}^{+}$ to indicate the first component of the eigenvector associated to the eigenvalue $\lambda^{+}$. Approximating to leading order of $\tau_{w}^{-1}$ the eigenvectors in Eq (28), we obtain the eigenvectors

[TABLE]

Then, using Eqs (31), (32) and (33), we determine the linear filter:

[TABLE]

Interestingly, in contrast with synaptic filtering, the amplitude of the two modes are not equal. The amplitude of the slow mode (first term in Eq 34), whose timescale is proportional to $\tau_{w}$, decays proportionally to $\tau_{w}^{-1}$ with respect to the fast mode, when $\tau_{w}\ll\tau_{m}\left(2g_{w}+1\right)$. Therefore, the area under the linear filter corresponding to this mode is independent of $\tau_{w}$ for very large adaptation time constants:

[TABLE]

It follows that, if the adaptation timescale is increased, its relative contribution to the activity will decrease by the same factor, so that very slow adaptive currents will effectively be masked by the fast mode.

Equilibrium activity

The two systems possess a non-trivial equilibrium state at which the input current of all units stays constant. Since all units are statistically equivalent, the equilibrium activity is the same for all units. For synaptic filtering, the input current at equilibrium is given by a transcendental equation, that is obtained by setting to zero the left hand side of Eq (17):

[TABLE]

This equilibrium coincides with the fixed point of the system without synaptic filtering.

For adaption, instead, from Eq (16) we obtain that the equilibrium is determined by

[TABLE]

We further assume unless otherwise specified that the fixed point of the system is in the linear regime of the transfer function, so that $\phi\left(x\right)=x-\gamma$. In that case $x_{0}=\left(J\left(C_{E}-gC_{I}\right)-g_{w}\right)\left(x_{0}-\gamma\right)$, so that larger adaptation coupling corresponds to weaker input currents, i.e. decreasing stationary firing rate. The adaptation time constant does not affect the fixed point.

Dynamics of homogeneous perturbations

We study the neuronal dynamics in response to a small perturbation uniform across the network

[TABLE]

Synaptic filtering

Linearizing Eq 17 we obtain

[TABLE]

where we use the notation $\phi_{0}^{\prime}:=\left.\frac{d\phi\left(x\right)}{dx}\right\rvert_{x_{0}}$. Because the perturbation $\delta x$ in Eq (39) is independent of $j$, using Eq (19) the dynamics for all units are equivalent to the population-averaged dynamics and are given by

[TABLE]

From Eq (40) we can define the dynamic matrix

[TABLE]

The only difference in the linearized dynamics of the population-averaged current with respect to the single neuron dynamics (Eq 22) is the non-diagonal entry $\phi_{0}^{\prime}J\left(C_{E}-gC_{I}\right)$. When either the derivative at the fixed point cancels, or when the total effective input is zero, the population dynamics equals the dynamics of a single neuron. The eigenvalues of the population-averaged dynamics are

[TABLE]

and the eigenvectors

[TABLE]

For very large synaptic time constants, the eigenvalues are approximated to leading order as

[TABLE]

Approximating as well the eigenvectors to leading order, we obtain

[TABLE]

the filter of the linear response to weak homogeneous perturbations reads:

[TABLE]

Note that the amplitude of the two exponential terms is the same, independently of the effective coupling and time constants.

Adaptation

For the system with adaptive neurons, the linearized system reads

[TABLE]

As for the network with synaptic filtering, the dynamics of the perturbation are equivalent for each unit, so that we can write down the dynamic matrix for the population-averaged response to homogeneous perturbations

[TABLE]

The difference with respect to the linear single neuron dynamics (Eq 50) is that the effective recurrent coupling appears now in the first diagonal entry of the dynamic matrix.

When the fixed point is located within the linear range of the transfer function, the derivative is one, so that we do not further specify the factor $\phi_{0}^{\prime}$ in the following equations. Consequently, the dynamics of the system to small perturbations do not depend on the exact value of the fixed point, which does not hold for more general transfer functions.

The eigenvalues of the system read

[TABLE]

with eigenvectors

[TABLE]

In the limit of very slow adaptation, given that the two eigenvalues are real, they can be approximated to leading order as

[TABLE]

and the corresponding eigenvectors read

[TABLE]

Therefore, if the perturbation is stable (see next section) we can write down the corresponding linear filter as

[TABLE]

The area under the slow mode is again independent of the adaptation time constant in this limit,

[TABLE]

Stability of homogeneous perturbations

The equilibrium point is stable when the real part of all eigenvalues is negative. Equivalently, in a two dimensional system –as it is the case for the population-averaged dynamics–, the dynamics are stable when the trace of the dynamic matrix is negative and the determinant positive.

Synaptic filtering

In the system with synaptic filtering, the trace and determinant are

[TABLE]

The trace is therefore always negative. The determinant is positive, and therefore the population-averaged dynamics are stable, when the effective coupling $J\left(C_{E}-gC_{I}\right)$ is smaller than unity. In contrast, if the effective coupling is larger than unity, i.e. if positive feedback is too strong, the equilibrium firing rate is unstable, so that any small perturbation to the equilibrium firing rate will lead the system to a different state. Right at the critical effective coupling, one eigenvalues is zero and the other one equals $\text{Tr}_{s}$, implying that the population-averaged dynamics undergo a saddle-node bifurcation. Beyond the bifurcation, the network reaches a state where the firing rates of all neurons saturate.

Adaptation

In the adaptive population dynamics, the recurrent connectivity has a different effect on the stability of the adaptive population dynamics. The trace and determinant of the dynamic matrix are

[TABLE]

Both the timescale $\tau_{w}$ and the strength $g_{w}$ of adaptation affect the trace and determinant of the dynamic matrix, and therefore the stability. The system is unstable if the determinant is negative (one positive and one negative real eigenvalue) or if the determinant is positive and the trace is positive. The determinant is negative, and therefore the system becomes unstable through a saddle-node bifurcation, when $J\left(C_{E}-gC_{I}\right)>1+g_{w}$. Note that the adaptation strength increases the stability of the system: a stronger positive feedback loop is required to destabilize the fixed point, in comparison to the network with synaptic filtering. The determinant and trace are positive if $J\left(C_{E}-gC_{I}\right)<1+g_{w}$ but $J\left(C_{E}-gC_{I}\right)>1+\frac{\tau_{m}}{\tau_{w}}$, respectively, leading to a Hopf bifurcation: the system produces sustained marginal oscillations at the bifurcation in response to small perturbations around the fixed point. Beyond the Hopf bifurcation, the oscillations are maintained in time, unless the system shows a fixed point when all neurons saturate ($x_{0}=\frac{1}{1-g_{w}}\left(J\left(C_{E}-gC_{I}\right)\phi_{max}+g_{w}\gamma\right)$). This fixed point exists if $x_{0}>\phi_{max}+\gamma$.

Heterogeneous activity

We next study the network dynamics beyond the population-averaged activity, along modes where different units have different amplitudes. We study perturbations of the type

[TABLE]

We define the $2N$-dimensional vector $\textbf{x}=\left(\delta x_{1},...,\delta x_{N}^{1},\delta w_{1}^{1},...,\delta w_{N}^{1}\right)^{T}$. Since the dynamics of each unit is now different, the dynamic matrix of the linearized system, $A$, is described by a squared matrix of dimensionality $2N$. Therefore, the perturbations generate dynamics along $2N$ different modes whose timescales are determined by the eigenvalues of the matrix $A$. The eigenvalues are determined by the characteristic equation $\left|A-\lambda I\right|=0$. In order to calculate these eigenvalues, we make use of the following identity which holds for any block matrix $Z=A-\lambda I$, that is composed by the four square matrices P,Q, R, and S and the block S is invertible:

[TABLE]

Consequently, if we set Eq (67) to zero, since we assumed that $\left\lvert S\right\rvert\neq 0$, we obtain

[TABLE]

The identity in Eq (67) can be shown by using the decomposition

[TABLE]

together with the fact that when a non-diagonal block is zero. The determinant of such a matrix is the product of determinants of the diagonal blocks.

Synaptic filtering

The dynamical matrix for the network with synaptic filtering, obtained by linearizing Eqs (17), is

[TABLE]

The matrix J is the connectivity matrix. Again, we assume in the following that the fixed point is located in the linear range of the transfer function, so that $\phi_{0}^{\prime}=1$.

The characteristic equation, obtained by combining Eqs (68) and (70), reads

[TABLE]

where $\lambda_{J}$ are the eigenvalues of the connectivity matrix. Solving for $\lambda_{J}$ we obtain the equation which maps the eigenvalues of the synaptic filtering network dynamics $\lambda_{s}$ onto the eigenvalues of the connectivity matrix $\lambda_{J}$,

[TABLE]

In contrast, solving for the eigenvalues of the dynamic matrix $\lambda_{s}$ we obtain the inverse mapping

[TABLE]

In other words, Eqs 73 and  72 constitute two different approaches to assessing the stability of the system [37]. One approach is to examine whether the domain of eigenvalues $\lambda_{s}$ resulting from Eq (73) intersect the line $\text{Re}\left(\lambda_{s}\right)=0$ (Fig 3, insets in B). The eigenvalues $\lambda_{J}$ of the connectivity matrix are distributed within a circle in the complex plane, whose radius is proportional to the synaptic strength, $\lambda_{J}<J\sqrt{C_{E}+g^{2}C_{I}}$ plus an outlier real eigenvalue at $J\left(C_{E}-gC_{I}\right)$ that corresponds to the homogeneous perturbations studied above (see [36]). We focus in this section on the bulk of eigenvalues that corresponds to modes of activity with different amplitudes for different units. We can therefore parametrize the eigenvalues $\lambda_{J}$ as

[TABLE]

and introduce the parametrization into Eq (73) to obtain an explicit expression for the curve that encloses the eigenspectrum $\lambda_{s}$. Note that in an abuse of notation, we denote the limits of the eigenspectrum as $\lambda$ and not the eigenvalues themselves that constitute the eigenspectrum.

The alternative approach is to use the inverse mapping from the eigenvalues $\lambda_{s}$ to the eigenvalues of the connectivity $\lambda_{J}$, by mapping the line $\text{Re}\left(\lambda_{s}\right)=0$ into the space of eigenvalues $\lambda_{J}$ (Fig 8). More specifically, the line $\text{Re}\left(\lambda_{s}\right)=0$ can be parametrized as

[TABLE]

and introduced into Eq (72). In this case, the stability is assessed by whether the eigenspectrum of the connectivity matrix J crosses the stability boundary or not (insets Fig 8). This alternative approach is useful for some calculations due to the simple geometry of the connectivity eigenspectrum $\lambda_{J}$.

Taking the alternative approach, introducing Eq (75) into Eq (72), we obtain the stability bound in the complex plane of eigenvalues $\lambda_{J}$:

[TABLE]

The first point of the stability curve $\lambda_{J}^{sb}\left(\omega\right)$ intersecting with a circle of increasing radius centered at the origin is the closest point of the curve to the origin, i.e. the minimum of $\left|\lambda_{J}^{sb}\right|^{2}$ with respect to $\omega$. The squared distance to the origin is

[TABLE]

whose minimum happens trivially at $\omega=0$, $\lambda_{J}=1$ (Fig 8A). In conclusion, the system is unstable if

[TABLE]

Note that this is the same condition as in the case without synaptic filtering, The synaptic filtering system approaches the no-filtering system when $\tau_{s}\to 0$. Although we are considering in this work synaptic timescales that are larger than the membrane time constant, the analysis is valid for arbitrarily fast synaptic time constants. In that limit, the stability curve in Eq (76) approaches the curve $\lambda_{J}^{sb}=1$, retrieving the stability boundary found in [29].

To study the limit of slow synaptic time constant, $\tau_{s}\gg\tau_{m}$, we analyze the direct approach, i.e. study how the parameters of adaptation modify the eigenspectrum of the dynamic matrix $A_{s}$ in the complex plane of eigenvalues $\lambda_{s}$. To this end, we introduce the parametrized connectivity eigenspectrum (Eq 74) into Eq (73), and approximate it to leading order of $\frac{\tau_{m}}{\tau_{s}}$. We obtain that the eigenspectrum of eigenvalues $\lambda_{s}$ are enclosed by the curves

[TABLE]

The equations above approximate the full eigenspectrum by two disjoint circles of radius $\tau_{s}^{-1}J\sqrt{C_{E}+g^{2}C_{I}}$, the one corresponding to the $\lambda_{s}^{+}$ eigenvalues centered at $-\frac{1}{\tau_{s}}$, and the other circle $\lambda^{-}$ centered at $-\frac{1}{\tau_{m}}+\frac{1}{\tau_{s}}$. The circle centered closer to the instability bound, $\lambda_{s}^{+}$ sets the slow timescales of the network, and its associated timescale is proportional to $\tau_{s}$. This gives an intuitive explanation to why the network timescale scales linearly with the synaptic time constant (Fig 5).

Adaptation

For adaptation, we repeat the same procedure as for the synaptic filtering to determine the stability to heterogeneous perturbations. The dynamical matrix reads

[TABLE]

Using Eqs (68) and (81) we can obtain the characteristic equation. Solving for $\lambda_{J}$ we obtain the mapping between the $\lambda_{w}$ eigenvalues and the connectivity eigenvalues

[TABLE]

while solving for $\lambda_{w}$ we obtain the expression for the inverse mapping:

[TABLE]

We first explore the inverse mapping. Inserting the parametrization in Eq (75) into Eq (82), the stability curve in the complex plane of connectivity eigenvalues reads

[TABLE]

The bifurcation parameters can then be found by determining the closest point of the stability boundary to the origin. The corresponding value of $\omega$ determines the oscillatory frequency of the first unstable mode. This value can be zero, corresponding to a zero-frequency bifurcation, which generally leads to slowly fluctuating activity referred to as rate chaos ([29], [35], [40], [33]). Alternatively, when the parameter $\omega$ that minimizes the norm of $\lambda_{J}^{sb}$ is non-zero, the system undergoes a Hopf bifurcation.

It is useful to consider the different geometries of the stability curve in Eq (84) in order to identify the closest point of the curve to the origin. Note that the curve shows symmetry with respect to the real axis, $\lambda_{J}^{sb}\left(-\omega\right)=\lambda_{J}^{sb*}\left(\omega\right)$.

The curve might cross the real axis $\text{Im}\left(\lambda_{J}\right)=0$ either in one or two different values of $\left|\omega\right|$. Solving $\text{Re}\left(\lambda_{J}^{sb}\right)=0$, we find that the curve crosses twice the real axis, when $\tau_{m}<\tau_{w}g_{w}$ (Fig 8 Bi). In that case, one crossing is the point $\tau_{m}\lambda_{J}=1+\frac{\tau_{m}}{\tau_{w}}$ and the other $\tau_{m}\lambda_{J}=1+g_{w}$. This second intersection corresponds to $\omega=0$. Therefore, it is clear that, since the first crossing of the real axis is closer to the origin than the point at $\omega=0$, the bifurcation necessarily occurs at non-zero frequency for $\tau_{m}<\tau_{w}g_{w}$.

When the curve crosses only once the zero axis, the point $\lambda_{J}=1+g_{w}$, corresponding to a zero-frequency, is not necessarily the closest one to the origin (Fig 8 Bii). One approach to determine analytically whether the system undergoes a Hopf or a zero-frequency bifurcation is to look at the curvature at the point $\omega=0$ and compare it to the curvature of a circle with radius $1+g_{w}$. To do so, we approximate both the stability line and the circle by a parabola, and compare their curvatures (dashed curve, Fig8 Bii and Biii). First, we write the stability boundary in its implicit form, $\lambda_{J}^{sb}:=x_{J}^{sb}+iy_{J}^{sb}$, as

[TABLE]

Then, we consider small deviations of the coordinates $x^{sb}=1+g_{w}+\epsilon_{x}$ and $y^{sb}=\epsilon_{y}$. If we approximate up to first order of $\epsilon_{x}$ and second order of $\epsilon_{y}$ we obtain the parabola

[TABLE]

Repeating the same procedure for the circle of eigenvalues, with radius $r=1+g_{w}$ we obtain $\epsilon_{y}^{2}=2\left(1+g_{w}\right)\epsilon_{x}+O\left(\epsilon_{x}^{2}\right)$. By requiring the circle of eigenvalues to be interior to the boundary curve (for the same $\epsilon_{x}$, $\epsilon_{y,\text{circle}}^{2}<\epsilon_{y,\text{sb}}^{2}$), we obtain that the instability parabola is exterior to the circle, therefore the system undergoes a zero-frequency bifurcation (Fig 8C), when

[TABLE]

which simplifies to

[TABLE]

In the limit of the adaptation timescale approaching the membrane time constant, the left side of the inequality above approaches one. Introducing this value in Eq 88, we find that for adaptive couplings stronger than $g_{w}>\sqrt{5}-2$ only a Hopf bifurcation is possible.

Dynamical Mean Field Theory

The linearization of the dynamical system from the previous section is only valid up to the instability boundary. A commonly used method to study the dynamics that arise beyond the bifurcation is dynamical mean field theory (DMFT) [29, 34, 41, 40, 39, 35, 33]. DMFT approximates the deterministic input to each element of the system by a Gaussian stochastic process, whose first and second moment are determined self-consistently.

The dynamics of the $i$-th neuron in the synaptic and adaptive network are approximated as

[TABLE]

[TABLE]

where $\xi_{i}\left(t\right)$ is a Gaussian variable. In the thermodynamic limit, the noise sources are independent between neurons, so that for $i\neq j$ $\left[\xi_{i}\left(t\right)\xi_{j}\left(t^{\prime}\right)\right]=0$.

The next step is to determine the self-consistent equations, that links the distribution of $\xi_{i}$ to the statistics of the original system in Eqs (16) and (17). First, we relate the statistics of the noise, currents $x_{i}$ and rates $\phi\left(x_{i}\right)$ based on the dynamics. Then, we close the equations by explicitly assuring that the transfer function relates the currents and the rates.

To determine the first moment of the noise, we apply that $\xi_{i}\left(t\right)=\sum_{j}J_{ij}\phi\left(x_{j}\left(t\right)\right)$ and average over the population, as in [33]. The first moment of the noise then obeys

[TABLE]

We calculate next the relation for the second moment of the noise, which again is the same as in [33]:

[TABLE]

where $C\left(\tau\right)=\left\langle\phi_{i}\left(t\right)\phi_{i}\left(t+\tau\right)\right\rangle$ .

These equations show that the first and second moment of the Gaussian sources do not depend on the identity of neuron $i$, so that all neurons are statistically equivalent. Thus, we can reduce the full $2N$-deterministic system to a two-variable stochastic system, describing a prototypical neuron in the network.

The equations (91) and (92) describe how the noise is related to the properties of the connectivity and the statistics of the rates $\phi\left(x\right)$. The next step is to calculate how the first and second moment of the noise are related to the statistics of the input current, which we write as $\mu:=\left[x_{i}\right]$ for the first moment and $\Delta\left(\tau\right):=\left[x_{i}\left(t\right)x_{i}\left(t+\tau\right)\right]-\mu^{2}$ for the second moment.

For the mean of the input current, averaging over units Eqs (89) and (90) and introducing the result in (91) for the synaptic and adaptive system respectively, we obtain

[TABLE]

By differentiating twice $\Delta\left(\tau\right)$ with respect to the lag $\tau$ and using Eqs (89) and (92), as in [29, 33] we obtain:

[TABLE]

where $Q_{s}\left(\tau\right):=\int_{-\infty}^{+\infty}dth_{s}\left(t\right)h_{s}\left(t+\tau\right)$ is the autocorrelation function of the single neuron filter $h_{s}$ (Eq.23). Equivalently, for the adaptive system, using Eq (90) and (92) we obtain

[TABLE]

where we define in relation to Eq (6) $h_{w}^{sym}\left(\tau\right)=h_{w}\left(\left|\tau\right|\right)$, and the autocorrelation function of the adaptive filter $Q_{w}:=\int_{-\infty}^{+\infty}dth_{w}\left(t\right)h_{w}\left(t+\tau\right)$.

Secondly, in order to close the self-consistent description, we can link the statistics of the rates $\phi_{i}\left(t\right)$ with the statistics of the currents $x_{i}\left(t\right)$ by writing the input currents explicitly as Gaussian variables. We can write down the input current at time $t$ and $t+\tau$ explicitly as (see [34]):

[TABLE]

This explicit construction in terms of Gaussian variables $z_{1}$, $z_{2}$ and $z_{3}$ realizes the constraints $\left[x^{2}\left(t\right)\right]-\mu^{2}=\Delta\left(0\right)$, $\left[x^{2}\left(t+\tau\right)\right]-\mu^{2}=\Delta\left(0\right)$ and $\left[x\left(t\right)x\left(t+\tau\right)\right]-\mu^{2}=\Delta\left(\tau\right)$. Now, explicitly calculating the first moment of the rates by replacing the average for a Gaussian integral and using Eq (97) we obtain

[TABLE]

where we use the short-hand notation $\int Dz=\int_{-\infty}^{+\infty}\frac{e^{-\frac{z^{2}}{2}}}{2}dz$.

For the second moment, introducing Eqs (97) and (98) into the definition of autocorrelation function of the rate, we get

[TABLE]

Therefore, in order to determine the self-consistent solution, we need to find a mean and autocorrelation function for the currents that satisfy both Eqs (99) and  (100) and Eqs (93) and (95) (for the synaptic system) and Eqs (94) and (96) (for the adaptive system). Once the statistics of the currents and rates are known, it is straight-forward to obtain the statistics of the noise, using Eqs (91) and (92).

In previous works [29, 35, 40, 38, 33] it was possible to further simplify the self-consistent equations because the resulting analogous equation to Eqs (95) and Eq (96) was a conservative system. However, in the networks studied here, synaptic filtering and adaptation add the convolutional terms in Eqs (95) and Eq (96) that make the system non-conservative. Therefore, we followed an alternative approach and found the solutions to the self-consistent equations using an iterative scheme, that circumvents solving directly the integral equations.

Iterative scheme

We solve the self-consistent equations numerically following a single-unit iterative scheme, as in [42, 43, 44, 41]:

• •

First, we simulate the dynamics in Eqs (89) and (90) assuming white Gaussian noise with a certain mean $\left[\xi\right]^{\left(0\right)}$ and autocorrelation function $\left[\xi\left(t\right)\xi\left(t+\tau\right)\right]=\left(\sigma_{\xi}^{\left(0\right)}\right)^{2}\delta\left(\tau\right)$.

• •

We calculate the autocorrelation functions of the firing rate and input currents empirically, $\mu^{\left(0\right)}$, $\Delta^{\left(0\right)}$, $\left\langle\phi\right\rangle^{\left(0\right)}$ and $C^{\left(0\right)}\left(\tau\right)$.

• •

We simulate in the new iteration $k+1$ the noise following the self-consistent statistics obtained in the previous iteration, as indicated by Eqs (91) and (92)

[TABLE]

In order to numerically generate a Gaussian variable with autocorrelation function $G\left(\tau\right)$, we first generate the noise in the Fourier domain, where each frequency component of the noise is given by

[TABLE]

where $\tilde{G}\left(\omega\right)$ denotes the Fourier transform of the target autocorrelation function, and $\psi$ is a random variable with uniform probability density in the range $\left[-\pi,\pi\right]$.

• •

We repeat the previous step until the values $\mu^{\left(k\right)}$, $\Delta^{\left(k\right)}$, $\left\langle\phi\right\rangle^{\left(k\right)}$ and $C^{\left(k\right)}\left(\tau\right)$ do not vary beyond a certain tolerance for new iterations.

We find that such an iterative method applied to the systems studied here converges to a solution for the parameters of the noise after a few iterations, independently of the noise properties used in the initial step.

The drawbacks of this iterative scheme are that the two-dimensional system needs to be simulated several times at each iteration in order to determine the first and second order statistics of the input current and the firing rate, which is in general a computationally costly operation. We also find that the method converges more robustly to the solution (given the fact that both the trial length in the simulation and the number of trials are finite), at the expense of initial speed convergence, when the first and second moments of the noise are only partially updated at each iteration, so that

[TABLE]

and similarly for the second-moment equation, where $\alpha$ is a parameter between zero and one. In this work, we used $\alpha=0.6$.

This method is inefficient for very large adaptation and synaptic time constants, since it requires simulating with both a fine temporal resolution (faster than the membrane time constant) over very large intervals (much larger than the slow adaptive/synaptic timescale). Another drawback of the iterative method is that its convergence is based on the assumption that smooth changes in the noise statistics lead to smooth changes in the statistics of the firing rates. In general, close to a bifurcation, this requirement may not hold.

Dynamics with intrinsic noise

We next study how white Gaussian noise, independent between neurons and intrinsic to each unit in the network, affects the dynamics of the system. On the mean-field level, this is equivalent to studying a network where each neuron spikes at a Poisson process whose rate varies in time as $\phi\left(x_{i}\left(t\right)\right)$ [35]. The additional input to each neuron, whose dynamics are given in Eqs (2) and (3), is now

[TABLE]

where $\left[\eta_{i}\right]=0$, and $\left[\eta_{i}\left(t\right)\eta_{j}\left(t+\tau\right)\right]=\delta_{ij}\frac{\sigma_{\eta}^{2}}{2}\delta\left(\tau\right)$, and Gaussian distributed. The DMF equations are derived following the same steps as in the absence of intrinsic noise. The stochastic variable $\xi\left(t\right)$ is the sum of the recurrent input and the intrinsic noise. Its first moment remains unchanged:

[TABLE]

which is the same result as Eq (91). The second moment of the stochastic process is the sum of the variance generated by the recurrent connections and the variance of the intrinsic noise

[TABLE]

Accordingly, the iterative scheme now takes into account the equation above, so that the equation for the second moment of the self-consistent relation (Eq 102) reads when there is intrinsic noise

[TABLE]

Adding white noise produces a discontinuity in the derivative of the autocorrelation function of the firing rates at zero lag (Fig 7 A and D). This can be shown by integrating explicitly both sides of the Eqs (95) and (96) around zero when the external noise is added. It results in the condition

[TABLE]

Since the autocorrelation function is a symmetric function, $\dot{\Delta}\left(0^{+}\right)=-\dot{\Delta}\left(0^{+}\right)$, leading to

[TABLE]

Thus, the autocorrelation function of the input current decays linearly at zero time lag with a slope proportional to the external noise intensity, which also extends to the autocorrelation function of the firing rate.

Definition of the timescale of the activity

The activity of multivariable dynamical systems ranges over several timescales. In particular, for stable linear systems, the timescales of the activity are given by the inverse of the absolute values of the real part of the eigenvalues. As we showed before, for single adaptive or synaptic neurons, the activity consists of two modes that evolve at two different timescales. However, the relative contribution of each of the excited modes can make one timescale more predominant than the other, as it happens for slow adaptation time constant, which becomes effectively undetectable in the single neuron dynamics.

In this work, we calculate the timescale of the activity for linear systems as the average of the timescales of the activated input current modes, weighed by their contribution (Fig 1). For a linear system with filter $h\left(t\right)=\sum_{k}a_{k}e^{-\frac{t}{\tau_{k}}}$, the correlation time is

[TABLE]

For large networks, which are high-dimensional non-linear systems, we define the main timescale of the activity as the time lag at which the autocorrelation function has decayed to a fraction $e^{-\frac{1}{2}}$ of its maximum (Figs 5 and 7):

[TABLE]

where $E[C\left(\tau\right)]$ is the envelope of the autocorrelation function, calculated as the norm of its analytic signal, computed using the Hilbert transform. This corresponds to the width of the envelope at which the autocorrelation decays to $e^{-0.5}$ of its value. For an exponentially decaying correlation function, this measure corresponds to the decay time constant. For a Gaussian envelope, this measure would correspond to two times its standard deviation, $2\sigma$.

Acknowledgments

This work was funded by the Ecole des Neurosciences de Paris Ile-de-France, the Programme Emergences of City of Paris, Agence Nationale de la Recherche grants ANR-16-CE37-0016-01 and ANR-17-ERC2-0005-01, and the program “Investissements d’Avenir” launched by the French Government and implemented by the ANR, with the references ANR-10-LABX-0087 IEC and ANR11-IDEX-0001-02 PSL* Research University. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.