How single neuron properties shape chaotic dynamics and signal transmission in random neural networks

TL;DR

This paper develops a dynamical mean-field theory to understand how intrinsic neuron properties influence chaotic dynamics and signal transmission in large, randomly connected neural networks, revealing resonance phenomena and enhanced signal processing.

Contribution

It introduces a novel theoretical framework linking single neuron dynamics to network chaos and signal transmission, applicable to realistic neural models.

Findings

Power spectrum of chaotic activity emerges from nonlinear sharpening of single-unit responses.

Networks with adaptation exhibit resonant chaos with narrow-band oscillations.

Signal transmission is enhanced by adaptation in chaotic regimes.

Abstract

While most models of randomly connected networks assume nodes with simple dynamics, nodes in realistic highly connected networks, such as neurons in the brain, exhibit intrinsic dynamics over multiple timescales. We analyze how the dynamical properties of nodes (such as single neurons) and recurrent connections interact to shape the effective dynamics in large randomly connected networks. A novel dynamical mean-field theory for strongly connected networks of multi-dimensional rate units shows that the power spectrum of the network activity in the chaotic phase emerges from a nonlinear sharpening of the frequency response function of single units. For the case of two-dimensional rate units with strong adaptation, we find that the network exhibits a state of "resonant chaos", characterized by robust, narrow-band stochastic oscillations. The coherence of stochastic oscillations is maximal at the onset of chaos and their correlation time scales with the adaptation timescale of single units. Surprisingly, the resonance frequency can be predicted from the properties of isolated units, even in the presence of heterogeneity in the adaptation parameters. In the presence of these internally-generated chaotic fluctuations, the transmission of weak, low-frequency signals is strongly enhanced by adaptation, whereas signal transmission is not influenced by adaptation in the non-chaotic regime. Our theoretical framework can be applied to other mechanisms at the level of single nodes, such as synaptic filtering, refractoriness or spike synchronization. These results advance our understanding of the interaction between the dynamics of single units and recurrent connectivity, which is a fundamental step toward the description of biologically realistic network models in the brain, or, more generally, networks of other physical or man-made complex dynamical units.

Tables

Table 1. Table 1: Parameters of the models in Fig. 8 . Matrix A A \mathrm{A} defining the rate model for the different examples in Fig. 8 .

Fig. 8a,b,c	Fig. 8d,e,f

Equations

$${\dot{x}}_i^{\alpha }(t)=\sum _{\beta =1}^D{\mathrm{A}}^{\alpha \beta }{x}_i^{\beta }(t),$$

$${\dot{x}}_i^{\alpha }(t)=\sum _{\beta =1}^D{\mathrm{A}}^{\alpha \beta }{x}_i^{\beta }(t)+{\delta }^{\alpha 1}\left(\sum _{j=1}^N{J}_{ij}\varphi (x_j^1(t))+I_i(t)\right)$$

$$J_{ij}\sim \mathcal{N}\left(0,g^2/N\right)$$

$${\dot{x}}_i(t)$$

$${\dot{a}}_i(t)$$

$$g_c(\gamma ,\beta )=\{\begin{array}{ll}\sqrt{1-\gamma (\gamma +2\beta )+2\sqrt{{\gamma }^2\beta (2\gamma +2\beta +2)}},& \beta >{\beta }_H(\gamma )\\ 1+\beta ,& \beta \le {\beta }_H(\gamma )\end{array}$$

$$f_m=\frac{1}{2\pi }\sqrt{-{\gamma }^2+\sqrt{\beta {\gamma }^2(\beta +2\gamma +2)}}.$$

$${\varphi }_{PL}(x)=\{\begin{array}{l}-1\text{for}x<-1\\ x\text{for}-1<x<1\\ 1\text{for}x>1\end{array}$$

$${\dot{x}}^{\alpha }(t)=\sum _{\beta =1}^D{\mathrm{A}}^{\alpha \beta }{x}^{\beta }(t)+{\delta }^{\alpha 1}\left(\eta (t)+I(t)\right),$$

$$⟨\eta (t)\eta (s)⟩=g^2⟨\varphi (x^1(t))\varphi (x^1(s))⟩.$$

$${\tilde{x}}^1(f)={\tilde{\chi }}_0(f)\left(\tilde{\eta }(f)+\tilde{I}(f)\right),$$

$${\tilde{\chi }}_0(f)={\left[{\left(2\pi if{\mathrm{I}}_D-\mathrm{A}\right)}^{-1}\right]}^{1,1},$$

$$S_x(f)=\tilde{G}(f)\left(S_{\eta }(f)+S_I(f)\right),$$

$$S_{\eta }(f)=g^2{S}_{\varphi (x^1)}(f).$$

$$\tilde{G}(f)=\frac{{\left|{\left[\text{adj}(2\pi if{\mathrm{I}}_D-\mathrm{A})\right]}^{1,1}\right|}^2}{{\prod }_{i=1}^D|2\pi if-{\lambda }_{\mathrm{A}}^i{|}^2},$$

$$g_c^2\underset{f}{\max }\tilde{G}(f)=1.$$

$$\dot{x}(t)=$$

$$\dot{a}(t)=$$

$$S_x(f)=\tilde{G}(f)\left(g^2{S}_{\varphi (x)}(f)+S_I(f)\right).$$

$$\tilde{G}(f)=\frac{{\gamma }^2+{\omega }^2}{{\omega }^4+(1+{\gamma }^2-2\beta \gamma ){\omega }^2+{\gamma }^2(1+\beta {)}^2},$$

$$\tilde{\chi }(f)=\frac{\tilde{G}(f)}{1-g^2\tilde{G}(f)}.$$

$$Q=\frac{f_p}{\Delta f_{\text{HM}}},$$

$$t_c=\frac{{\int }_0^{\infty }\tau \left|C_x(\tau )\right|d\tau }{{\int }_0^{\infty }\left|C_x(\tau )\right|d\tau }.$$

$$I_i(t)=A_I\cos \left(2\pi f_{I}t+{\theta }_i\right),$$

$$S_x(f)=S_{\text{bkg}}(f)+S_{\text{osc}}(f)\cdot =S_{\text{bkg}}+\sum _{k=1}^{\infty }{b}_k\left(\delta (f-k{f}_I)+\delta (f+k{f}_I)\right),$$

$$S_x(f_I)=A_{\text{bkg}}+A_{\text{osc}}\cdot =\frac{S_x(f_I-\Delta f)+S_x(f_I-\Delta f)}{2}+\frac{b_1}{\Delta f},$$

$$\text{SNR}(f_I)=\frac{A_{\text{osc}}}{A_{\text{bkg}}}.$$

$$\text{Var}(x)={\int }_{-\infty }^{\infty }{S}_{\text{bkg}}(f)df+2\sum _{k=1}^{\infty }{b}_k=:P_{\text{bkg}}+P_{\text{osc}}.$$

$$|{\tilde{\chi }}_{\beta }(f){|}^2=\frac{\tilde{G}(f)}{1-g^2\tilde{G}(f)},$$

$$S_x(f)=|{\tilde{\chi }}_{\beta }(f){|}^2\left(S_I(f)+S_{\eta }(f)\right).$$

$$\mathrm{SNR}(f_I)=\frac{A_{\text{osc}}}{A_{\text{bkg}}}=\frac{|{\tilde{\chi }}_{\beta }(f){|}^2{S}_I(f_I)}{|{\tilde{\chi }}_{\beta }(f){|}^2{S}_{\eta }(f_I)}=\frac{S_I(f_I)}{S_{\eta }(f_I)},$$

$$\mathrm{B}=\left(\begin{array}[]{c c c c}\mathrm{A}^{11}\mathrm{I}_{N}+\phi^{\prime}(0)\mathrm{J}&\mathrm{A}^{12}\mathrm{I}_{N}&\dots&\mathrm{A}^{1D}\mathrm{I}_{N}\\ \mathrm{A}^{21}\mathrm{I}_{N}&\mathrm{A}^{22}\mathrm{I}_{N}&\dots&\mathrm{A}^{2D}\mathrm{I}_{N}\\ \dots&\dots&\dots&\dots\\ \mathrm{A}^{D1}\mathrm{I}_{N}&\mathrm{A}^{D2}\mathrm{I}_{N}&\dots&\mathrm{A}^{DD}\mathrm{I}_{N}\end{array}\right)\quad,$$

Full text

I Introduction

Random network models of interacting dynamical elements (units, or nodes), are widely used across different scientific fields Dorogovtsev and Mendes (2002). Examples from biology include neural network models Sompolinsky et al. (1988), metabolic networks Barkai and Leibler (1997) and protein regulatory networks Kauffman (1993); Jeong et al. (2001); Pomerance et al. (2009). The dynamics of random networks also plays an important role in the study of epidemic outbreaks Pastor-Satorras and Vespignani (2001), social networks Newman et al. (2002), power grids Nishikawa and Motter (2015) and transportation networks Banavar et al. (1999), as well as in abstract physics systems such as soft spin models Sompolinsky and Zippelius (1981) and networks of oscillators Rodrigues et al. (2016). In all the above examples, the collective dynamics result from the interplay between network connectivity and the dynamics of single units.

A typical property of large random networks of nonlinear dynamical elements is the existence of a chaotic phase. In the context of neural networks, the rich dynamics of large random networks of neuron-like elements at the edge of chaos has been exploited to learn complex tasks involving generation of temporal patterns Maass et al. (2002); Jaeger and Haas (2004); Sussillo and Abbott (2009); Laje and Buonomano (2013); Nicola and Clopath (2017); DePasquale et al. (2018); Mastrogiuseppe and Ostojic (2018). In these and other related approaches, the chaotic behavior of the network mainly arises from the random interactions, whereas the dynamics of single elements are typically given by first-order differential equations. The simplicity of single elements allows to quantitatively determine the chaotic phase of the coupled elements using dynamical mean-field theory (DMFT) Sompolinsky et al. (1988), even in networks with more realistic connectivity structure Rajan and Abbott (2006); Kadmon and Sompolinsky (2015); Mastrogiuseppe and Ostojic (2017, 2018).

A fascinating question is what kind of activity emerges in brain-like networks, that are subject to additional biological constraints. Individual neurons exhibit rich multi-dimensional internal dynamics Benda and Herz (2003); Lundstrom et al. (2008); La Camera et al. (2006); Pozzorini et al. (2013) that are inconsistent with first-order equations. However, a theoretical understanding of the emergent activity patterns in networks of more realistic multi-dimensional dynamical elements is largely lacking. In particular, beyond one-dimensional and simple oscillator models Sompolinsky et al. (1988); Wieland et al. (2015); van Meegen and Lindner (2018), a self-consistent mean-field theory for fluctuations such as the autocorrelation function or the power spectrum of the network activity is still an unsolved theoretical problem. Here, we develop a theoretical framework that extends DMFT to multi-dimensional rate neurons. Using this framework, we show that the power spectrum of the network activity in the nonlinear, strongly coupled regime, emerges from a sharpening of the single-neuron frequency response function due to strong recurrent connections.

Our theory uses firing rate models with two or more variables per unit. While rate-based models Wilson and Cowan (1972); Deco et al. (2011) discard information on the exact spike-timing of single neurons, they have the advantage of being accessible to an analytical characterization of their dynamics. However, commonly-used one-dimensional rate models cannot fully capture the dynamics of the mean activity of a population of spiking neurons, such as the synchronization of neurons in response to a stimulus onset Mainen and Sejnowski (1995); Bair and Koch (1996); Devalle et al. (2017), an effect that is readily observed in integrate-and-fire model Knight (1972); Konig et al. (1996); Gerstner (2000); Brette and Guignon (2003); Schwalger et al. (2017). To capture rapid synchronization after stimulus onset in rate models, it is necessary to consider at least two equations per rate neuron Mattia and Del Giudice (2002); Schaffer et al. (2013); Montbrió et al. (2015). Multi-dimensional models also account for additional cellular mechanisms such as refractoriness Berry and Meister (1998), spike-frequency adaptation (SFA) Naud and Gerstner (2012); Schwalger and Lindner (2013); Deger et al. (2014); Schwalger et al. (2017), synaptic filtering Fourcaud and Brunel (2002); Schwalger and Schimansky-Geier (2008), subthreshold resonance Richardson et al. (2003) or for the effect of dendritic compartments Ostojic et al. (2015); Doose et al. (2016).

To be specific, we focus on SFA, the decrease of a neuron’s firing rate in response to a sustained stimulus, but our theory can also be applied to other phenomena. SFA is present in neurons at all stages of sensory processing, and is believed to play a crucial role for efficient coding of external stimuli Benda and Herz (2003). Moreover, SFA over multiple timescales represents an efficient solution for information transmission of sensory signals whose statistics change dynamically Fairhall et al. (2001); Lundstrom et al. (2008); Pozzorini et al. (2013). It is therefore of great interest to understand how adaptation and recurrent connections interact to shape network dynamics and signal transmission Mar et al. (1999); Ávila Åkerberg and Chacron (2009). If connections and adaptation are weak, the network dynamics can be largely understood within linear response theory. In particular, in the presence of signals and noise, linear response theory predicts that adaptation shapes signal and noise in precisely the same manner Deger et al. (2014), canceling the noise-shaping effect of adaptation Mar et al. (1999); Ávila Åkerberg and Chacron (2009); Lindner (2016). In contrast, in strongly coupled networks generating chaotic fluctuations Sompolinsky et al. (1988), linear response theory is not applicable and the effect of adaptation on the signal transmission in this case remains poorly understood. Here, we show that introducing adaptation into a strongly-coupled network of rate units shifts the network to a state of “resonant” chaos that is qualitatively different from the chaotic behavior of the network without adaptation. In this state, the network generates a stable rhythm corresponding to a narrow-band peak in the power spectrum which is robust against quenched disorder in adaptation parameters (heterogeneity). We show that in this new regime the network has two interesting functional properties: first, the correlation time increases with the adaptation timescale; second, the low-frequency power of the chaotic activity is strongly decreased, enabling a better transmission of slow signals.

This paper is organized as follows: In section II we present the microscopic network model and discuss its different dynamical regimes. In sections III and IV, we present the mean-field theory and analyze the chaotic regime. In each of these sections, we first present the general formalism and then apply the result to the case of adaptation with a single auxiliary variable. Finally, we study the functional properties of the resonant chaotic state, focusing on the correlation time (section V) and on the response to external input (section VI). Detailed derivations and an example of a higher-dimensional rate model representing multi-timescale adaptation are provided in the appendix.

II Microscopic model

We are interested in studying the dynamics of a randomly connected recurrent network of multi-dimensional firing-rate units where each unit is described by a set of $D$ variables $x_{i}^{1},\dotsc,x_{i}^{D}$. The first variable $x_{i}^{1}$ is an activation variable that defines the output rate $y$ via a nonlinear function $\phi$, i.e. $y_{i}(t)=\phi(x_{i}^{1}(t))$. More precisely, $\phi(x_{i}^{1}(t))$ should be interpreted as the deviation of the firing rate from some reference rate. Therefore, $\phi(x_{i}^{1}(t))$ can take both positive and negative values.

The remaining $D-1$ variables are auxiliary variables. In isolation, each unit obeys a system of $D$ first-order linear differential equations

[TABLE]

where the dot denotes the temporal derivative. In what follows, subscripts (in Latin letters) indicate the index of the unit in the network and run from 1 to $N$, while superscripts (in Greek letters) indicate the index of the variable in the rate model and run from 1 to $D$. The matrix $\mathrm{A}$ is assumed to be non-singular and to have eigenvalues with negative real parts. We assume that the rate $\phi(x_{j}^{1}(t))$ is the only signal that unit $j$ uses to communicate with other units. Conversely, the signals coming from other units only influence the variable $x_{i}^{1}$, i.e. the rate of unit $j$ is directly coupled only to the first variable of unit $i$. The choice of having the same variable sending and receiving signals is dictated by simplicity and is not necessary for the development of the theory. Unit $i$ receives input from all the other units, via a set of random connections $J_{ij}$, sampled i.i.d. from a Gaussian distribution with mean zero and variance $g^{2}/N$. When incorporating these assumptions, the network equations read

[TABLE]

where $\delta^{\alpha\beta}$ is the Kronecker delta symbol. The external input $I_{i}(t)$ is assumed to have stationary statistics and zero mean.

Network with adaptation.

As a biologically relevant two-dimensional example, we consider rate units that undergo firing rate adaptation. We measure time in units of the timescale of the first variable such that $\mathrm{A}_{11}=-1$ (i.e. time is considered dimensionless), and work with a single auxiliary variable that mediates adaptation. To ease our notation we drop the superscripts in the variables and write $x_{i}$ instead of $x_{i}^{1}$ and $a_{i}$ instead of $x_{i}^{2}$. The adaptation variable $a_{i}(t)$ of neuron $i$ is driven by the neuron activation variable $x_{i}(t)$ and provides negative feedback onto $x_{i}$. The network equations (Eq. 2) for the adaptation case become

[TABLE]

where the parameter $\gamma>0$ can be interpreted as the ratio of the timescales of the two variables $x$ and $a$, while $\beta>0$ is a parameter that controls the strength of adaptation.

Numerical simulations of the network with adaptation show that for low connection strength $g$, the network exhibits transient dynamics before it settles to a fixed point in which all $x_{i}$ are zero (Fig. 1a,b). By analyzing the stability of this fixed point (see appendix A for details), we find that, in the $N\rightarrow\infty$ limit, the critical value of $g$ at which stability is lost depends on the adaptation parameters via

[TABLE]

where $\beta_{H}(\gamma)=-1-\gamma+\sqrt{2\gamma^{2}+2\gamma+1}$. Notice that $g_{c}(\gamma,\beta)>1$ for all $\gamma,\beta>0$, i.e. adaptation stabilizes the dynamics since in the case without adaptation we have $g_{c}(0,0)=1$ Sompolinsky et al. (1988). Interestingly, the two different cases in Eq. 5 correspond to two different bifurcation types: for $\beta<\beta_{H}(\gamma)$ we find that at $g=g_{c}(\gamma,\beta)$ we have a saddle-node bifurcation. In contrast, if $\beta>\beta_{H}(\gamma)$, the system undergoes a Hopf bifurcation (see appendix A), whose characteristic frequency is given by

[TABLE]

Examples of eigenvalue spectra of the randomly coupled network close to these two different bifurcations are shown in the insets of Fig. 1: For $\beta>\beta_{H}(\gamma)$, the eigenvalue spectrum is deformed such that the eigenvalues with the largest real part are complex (Fig. 1a,c). In contrast, for $\beta<\beta_{H}(\gamma)$ the spectrum is convex, exhibiting a rightmost eigenvalue that is real (Fig. 1b,d). Above the bifurcation, i.e. for $g>g_{c}(\gamma,\beta)$, the network exhibits self-sustained, irregular fluctuations (Fig. 1c,d) that we will characterize in the next sections.

In all the simulations and numerical integrations, we choose $\phi(x)$ as a piecewise-linear function given by

[TABLE]

unless stated otherwise.

III Mean-field theory

General theory.

The dynamics of the $ND$-dimensional dynamical system in Eq. 2 for large $N$ is too high-dimensional to be studied at the microscopic level. In contrast, using dynamical mean-field theory Sompolinsky et al. (1988), we can find properties of the network dynamics that are independent of the specific connectivity realization. In what follows, we will assume that the external input $I_{i}(t)$ to each unit is an independent realization of the same Gaussian process. Following Sompolinsky et al. (1988), we approximate the network input to a representative unit $i$ with a Gaussian process $\eta$ and substitute the average over time, initial conditions and network realizations with the average over realizations of $\eta$. This approximation is valid in the large-$N$ limit, in which neurons become independent Schücker et al. (2016); Crisanti and Sompolinsky (2018). In the mean-field description, the activity of each individual unit in the network follows a realization of the following system of $D$ stochastic differential equations, to which we refer to as mean-field equations (see appendix C for more details)

[TABLE]

where $\eta(t)$ is a Gaussian process. The mean $\langle\eta(t)\rangle$ vanishes because the averaging over the Gaussian process statistics mimics the average over different neurons and network realizations, and the connections $J_{ij}$ in Eq. 2 are sampled from a Gaussian distribution with mean zero. Thanks to the fact that neurons become independent in the large-$N$ limit Sompolinsky et al. (1988), the average of the network input over network realizations is also zero (see appendix C for more details). On the other hand, the autocorrelation function $\langle\eta(t)\eta(s)\rangle$ needs to be determined self-consistently by imposing (cf. appendix C)

[TABLE]

Thanks to the mean-field approximations, we reduced the $ND$-dimensional, deterministic, nonlinear system of Eq. 2 to the D-dimensional, stochastic system of Eq. 8, which looks linear at first glance. However, the nonlinearity is important and is hidden in the self-consistent match of the second moment, as expressed by Eq. 9. The linear mathematical structure of Eq. 8 allows us to write, in the frequency domain

[TABLE]

where $\tilde{\chi}_{0}(f)$ is the linear response function (susceptibility) of the mean-field system (Eq. 8), which is equal to the linear response function of an uncoupled single neuron in the microscopic description (Eq. 1). For the linear dynamics given by Eq. 8, the linear response function $\tilde{\chi}_{0}(f)$ is given by

[TABLE]

where $\mathrm{I}_{D}$ is the $D$-dimensional identity matrix and the upper indices 1,1 indicate the first element of the first row of the matrix inside the square brackets.

In what follows, we assume that the external input $I(t)$ is stationary and zero-mean, and that the network is in the stationary regime. Therefore, the mean of all variables is equal to zero. The second-order statistics must be determined self-consistently. In the frequency domain, this requires a self-consistent determination of the power spectral density (“power spectrum” for short) $S_{x}(f)$ of the activation variable $x^{1}$, defined as the Fourier transform of the autocorrelation function, $S_{x}(f)=\int_{-\infty}^{\infty}e^{-2\pi if\tau}\langle x^{1}(t+\tau)x^{1}(t)\rangle\,d\tau$. Using the squared modulus of the linear response function \tilde{G}(f)\mathrel{\hbox to0.0pt{\raisebox{1.29167pt}{\cdot}\hss}\raisebox{-1.29167pt}{\cdot}}=|\tilde{\chi}_{0}(f)|^{2}, the power spectrum can be expressed as

[TABLE]

where $S_{\eta}(f)$ and $S_{I}(f)$ denote the power spectral densities of $\eta(t)$ and $I(t)$, respectively. Importantly, from Eq. 9 we have that $S_{\eta}(f)$ depends implicitly on $S_{x}(f)$ through the self-consistency condition

[TABLE]

The factor $\tilde{G}(f)$ can be expressed as a function of the matrix $\mathrm{A}$ as

[TABLE]

where $\text{adj}\left(2\pi if\mathrm{I}_{D}-\mathrm{A}\right)$ is the adjoint matrix of $\left(2\pi if\mathrm{I}_{D}-\mathrm{A}\right)$ and $\lambda_{\mathrm{A}}^{i}$ are the eigenvalues of $\mathrm{A}$. In appendix E, we show that knowing the maximum of $\tilde{G}(f)$ is sufficient to compute the critical value of the coupling $g_{c}$:

[TABLE]

Network with adaptation.

For the network with adaptation defined by Eqs. (3, 4), the mean-field equations read

[TABLE]

with $\langle\eta(t)\rangle=0$ and $\langle\eta(t+\tau)\eta(t)\rangle=g^{2}\langle\phi(x(t+\tau))\phi(x(t))\rangle$. The self-consistent equations (Eqs. (12, 13) ), reduce to

[TABLE]

The factor $\tilde{G}(f)$ can be calculated explicitly, yielding

[TABLE]

with $\omega=2\pi f$.

In the next section, we show how the qualitative features of the dynamics of the network in the fluctuating regime can be predicted by the properties of the single unit linear response function, as summarized in the factor $\tilde{G}(f)$.

IV Resonant chaos in networks with adaptation

The traditional approach in the DMFT literature is to consider the time-domain version of Eq. 12 Sompolinsky et al. (1988). Applying the inverse Fourier transform to Eq. 12 would lead to a differential equation of order $2D$. Unfortunately, by contrast with the case $D=1$, for the multi-dimensional case $D>1$ the dynamics is no longer conservative, which precludes the determination of the initial conditions (see Sompolinsky et al. (1988)). We propose an alternative approach to find a self-consistent solution to Eq. 12 in the Fourier domain. This approach is based on an iterative map, the fixed point of which is the self-consistent solution. Iterative methods have been proposed previously both in the context of spiking Dummer et al. (2014); Wieland et al. (2015) and rate-based networks Stern et al. (2014) using Monte-Carlo methods. Here, we use a semi-analytical iteration method that allows to rapidly solve for the self-consistent power spectrum, and hence to qualitatively understand several features of the network dynamics.

In the frequency domain, the linear transform associated with $\tilde{G}(f)$ is simple, whereas the nonlinearity $\phi(x)$ is difficult to handle. Concretely, we need to express $S_{\phi(x^{1})}$ as a functional of $S_{x}(f)$. This calculation can be performed semi-analytically for the piecewise-linear nonlinearity (a detailed treatment of the nonlinear step is given in appendix F). The idea of our iterative method is to start with an arbitrary initial power spectral density $S_{\phi(x^{1})}^{(0)}(f)$, which we choose to be constant (white noise). We then apply multiple iterations each consisting of a linear step followed by a nonlinear one (Fig. 2f). At each iteration, the linear step is simply a multiplication by $g^{2}\tilde{G}(f)$ and it allows us to compute $(S_{x})^{(n+1)}(f)$. The nonlinear step afterwards transforms $(S_{x})^{(n+1)}(f)$ into $S_{\phi(x^{1})}^{(n+1)}(f)$.

By studying the iterative map that defines the mean-field solution, we conclude that the power spectrum of the network activity emerges from a sharpening of the linear response function $\tilde{G}(f)$ of single units. The sharpening mainly arises from repeated multiplications with the factor $g^{2}\tilde{G}(f)$ in the iteration, which however is balanced by cross-frequency interactions and saturation effects of the nonlinear steps (see appendix B for a detailed discussion). As a result, the network activity exhibits the same frequency bands that are preferred by single neurons, albeit much narrower.

We apply the iterative method to solve the mean-field equations for the network with adaptation (Eq. 18), in the absence of external input ($I=0$). We find that if $g<g_{c}(\gamma,\beta)$, the power spectrum converges to zero, $S_{x}(f)\rightarrow 0$, at all frequencies. Therefore the mean-field variable $x$ is constantly equal to zero. This is consistent with the presence of a stable fixed point at zero and it indicates that, in the thermodynamic limit, the fixed point solution is the only possible one. In this regime we can calculate the mean-field linear response function $\tilde{\chi}(f)$ of the network, which is given by (see appendix E)

[TABLE]

i.e. the network has a sharper linear response function than the single units.

On the other hand, if $g>g_{c}(\gamma,\beta)$, the mean-field network is characterized by a nonzero, continuous power spectral density (Fig. 2). This is an indication that, at the microscopic level, the network is in a chaotic state Schuecker et al. (2018). However, we stress that a more rigorous proof of chaos would require the computation of the maximum Lyapunov exponent of the network, which we will not perform. In contrast to a network without adaptation Sompolinsky et al. (1988), we find that in the presence of adaptation the network can be in two qualitatively different chaotic regimes. For very weak and/or fast adaptation, the chaotic fluctuations are qualitatively the same as for the network without adaptation, i.e. the power spectrum is broad-band with maximum at $f=0$ (Fig. 2b). We refer to this regime as to the non-resonant regime. On the other hand, for strong and/or slow adaptation, the mean-field network settles in a new regime, characterized by an autocorrelation that decays to zero via damped oscillations and, equivalently, by a power spectrum that exhibits a pronounced resonance band around a nonzero resonance frequency $f_{p}$ (Fig. 2a). The decaying autocorrelation function and the continuous power spectral density are an indication that the network is – also in this regime – in a state of microscopic chaos. This new dynamical state, that we refer to as resonant chaos, is qualitatively different from the one of the non-resonant regime and from the one of the non-adaptive network.

Strikingly, whether the network settles in the resonant or in the non-resonant regime can be predicted purely based on the single-unit adaptation properties. More precisely, if $\beta<\beta_{H}(\gamma)$, the function $\tilde{G}(f)$ is monotonically decreasing with the frequency $f$, i.e. it exhibits a low-pass characteristic (Fig. 2b). This low-pass behavior of the single neuron is reflected by a power spectrum of the network that is also dominated by low frequencies, albeit less broad. The network power spectrum corresponds exactly to the non-resonant regime discussed above.

In contrast, if $\beta>\beta_{H}(\gamma)$, the single neuron response amplitude $\tilde{G}(f)$ exhibits a maximum at a nonzero frequency $f_{0}=\frac{1}{2\pi}\sqrt{-\gamma^{2}+\sqrt{\beta\gamma^{2}(\beta+2\gamma+2)}}$. Such a resonance peak is typical of a band-pass filter (Fig. 2a). The frequency $f_{0}$ is identical to $f_{m}$ (Eq. 6), which is derived from the imaginary part of the critical eigenvalue at the Hopf bifurcation (see appendix A). The single-neuron linear response characteristics are qualitatively preserved in the fluctuating activity of the recurrent network, which also exhibits a power spectral density dominated by a nonzero frequency $f_{p}$. This regime corresponds to the resonant regime discussed above. Interestingly, we find numerically that $f_{p}=f_{0}$, i.e. the resonance frequency is not affected by the introduction of recurrent connections (Fig. 2c, tested up to $g=5g_{c}(\gamma,\beta)$). We notice that the non-resonant and resonant regimes are consistent with the fixed point stability analysis of the network in the microscopic description. Indeed, the resonant and non-resonant regimes match the regions in which we observe Hopf or saddle-node bifurcations, respectively (appendix A).

Using simulations of the full microscopic network, we verify that the mean-field description is a good approximation of the system for large but finite $N$. In Fig. 2e we show that the probability density of the activation variable $x$ measured from the microscopic simulations matches the Gaussian distribution predicted by the mean-field theory, with relatively small finite-size effects that increase close to the criticality (see Fig. 2e, $g=1.5g_{c}(\gamma,\beta)$). Moreover, the mean-field solution provides a good description of the system for a wide range of adaptation parameters $\gamma\,,\,\beta$ (Fig. 2c).

Network with heterogeneous adaptation.

The narrow-band oscillations in networks of adapting neurons reported so far have been obtained for networks of identical neurons. The variability of physiological properties of real neurons, however, suggests that adaptation parameters differ among neurons. Heterogeneity of neuronal parameters is known to sensitively influence synchronization in networks of neural oscillators Golomb et al. (2001). In particular, mismatches of oscillation frequencies can impede the formation of neural rhythms. Does heterogeneity have a similar effect in strongly coupled random networks of adapting neurons? To address this question, we introduce a second source of disorder in the system by considering quenched randomness in the adaptation parameters. Specifically, we construct the heterogeneous network by adding Gaussian noise to the parameter $\beta$, i.e. by sampling $\beta\sim\mathcal{N}\left(\bar{\beta},\sigma_{\beta}^{2}\right)$ independently for each neuron. Numerical simulations of the network in the presence of heterogeneous adaptation show that the dynamics of the random network are surprisingly robust to this type of noise (Fig. 2d). Even for relatively high variability ($\sigma_{\beta}/\bar{\beta}$=0.5), the only effect is a barely visible increase of the power spectral density at low-frequencies (Fig. 2d, inset). In appendix D we derive the mean-field equations that correspond to the network with heterogeneous adaptation, and we compute the effective factor $\tilde{G}_{H}(f)$ in this case. The semi-analytical solution of the mean-field theory for heterogeneous adaptation predicts, similar to simulations, a stronger power at low frequencies than in the homogeneous case. However, the deviations predicted by the theory are smaller than the mismatch between theory and simulations for the homogeneous case, so that we could not perform a quantitative verification of the mean-field theory for the heterogeneous case.

We now focus on the resonant chaotic regime, that represents the novel dynamical state that emerges from the introduction of adaptation and study the functional properties of the network in this regime.

V Correlation time and coherence of the oscillations

While the resonance frequency in the resonant regime seems to depend solely on the single-neuron properties, the introduction of recurrent connections increases the coherence of the stochastic oscillations, i.e. decreases the width of the resonance band. The narrower the resonance band, the more coherent the oscillatory behavior will be. To quantify the increase of the oscillation coherence, we measure the quality factor (Q-factor) of the stochastic oscillations, defined as

[TABLE]

where $\Delta f_{\text{HM}}$ is the frequency width of the power spectrum $S_{x}(f)$ at the half-maximum. Intuitively, for a narrow-band oscillation, the quality factor quantifies the number of oscillation cycles during the characteristic decay time of the autocorrelation function. For a single neuron driven by white noise ($\langle\eta(t)\eta(t^{\prime})\rangle=\delta(t-t^{\prime})$), the single-neuron power spectrum of $x$ is proportional to $\tilde{G}(f)$. Compared to this reference shape, we find a higher Q-factor in the recurrent network (Fig. 3a), corresponding to a sharper resonance peak in the power spectrum (see also appendix B). When approaching the criticality from the chaotic phase, $g\rightarrow g_{c}(\gamma,\beta)^{+}$, the quality factor diverges (Fig. 3a), i.e. the dynamics approach regular oscillations.

While the Q-factor measures the decay time constant of the autocorrelation function relative to the mean oscillation period, it is also interesting to consider the absolute correlation time of the activity. As a measure of correlation time of a stochastic process we use the normalized first moment (center of mass) of the absolute value of the autocorrelation function (e.g. Lindenberg et al. (1989)),

[TABLE]

Since the Q-factor diverges when $g\rightarrow g_{c}(\gamma,\beta)$, in this limit the corresponding autocorrelation exhibits sustained oscillations with a diverging correlation time. Due to the increase of the Q-factor, the correlation time also diverges when $g\rightarrow g_{c}(\gamma,\beta)$ (Fig. 3a).

In the regime of slow adaptation, a single unit driven by white noise can have a larger correlation time than a recurrent network (Fig. 3b). This is due to the fact that in this regime the correlation time of the single unit driven by white noise is dominated by the long tail of the autocorrelation. The introduction of recurrent connections increases the oscillatory component, giving a larger “weight” to the short time lags, thus decreasing $t_{c}$. Nevertheless, the correlation time increases with the timescale of adaptation $\tau_{a}$ for both the single unit driven by white noise and the recurrent network (Fig. 3b). Note that the Q-factor goes to zero for very large adaptation timescale ($\gamma\rightarrow 0$), so that the dominant contribution to the correlation time in this regime is the non-oscillatory one.

VI Response to periodic stimuli

In order to go beyond the study of the spontaneous activity of the network, we consider its response to an external oscillatory signal. While signal transmission in linear systems is fully characterized by the frequency response function of the system and by the noise spectrum of the output, the situation is different in the nonlinear neural network that we study here. Similarly to previous approaches Rajan et al. (2010), we provide oscillatory input to each unit in the microscopic network, randomizing the phase (Fig. 4a)

[TABLE]

where $\theta_{i}\sim\mathrm{U}(0,2\pi)$. The corresponding power spectral density of the input is given by $S_{I}(f)=\left(A_{I}^{2}/4\right)\cdot\left(\delta(f-f_{I})+\delta(f+f_{I})\right)$. Thanks to the phase randomization, the network still reaches a stationary state and the mean $\langle x(t)\rangle$ remains at zero. Notice that even if in this case the input is non-Gaussian, the mean-field equation for the power spectrum (Eq. 18) is still valid. However, since $x$ is also not Gaussian anymore, in order to find the mean-field solution we need to modify our iterative scheme by splitting the activation variable $x$ into its Gaussian and its oscillatory part Rajan et al. (2010).

The presence of the input affects the dynamics of the mean-field network, quantified by the power spectral density (Fig. 4b). If the input is given while the network is in the chaotic regime ($g>g_{c}$), sharp peaks at the driving frequency $f_{I}$ and multiples thereof are elicited by the external input, standing out from a background power spectrum that is deformed compared to the case without the external input. For $f_{I}>f_{p}$, as in the example, the bumps of the background spectrum are slightly shifted toward larger values. The opposite happens if $f_{I}<f_{p}$. Notice that both this shift and the shaping of the chaotic activity are nonlinear effects due to the recurrent dynamics. As an additional nonlinear effect, the network activity also exhibits harmonics at the driving frequency of the external input.

To characterize the response to the external stimulus, we split the power spectrum $S_{x}(f)$ into an oscillatory component and a chaotic component that constitutes the background activity

[TABLE]

where $b_{k}$ are positive coefficients and we included the multiples of the driving frequency in order to account for the harmonics. To solve the mean-field equations numerically, we have to consider a finite frequency bin $\Delta f$ (in our numerical results, $\Delta f=0.001$). As a consequence, the heights of the delta peaks in the power spectrum in Eq. 24 are finite and depend on $\Delta f$. First, we will look at the transmission of the oscillatory signal near the driving frequency, i.e. how much of the peak in the power spectrum $S_{x}(f)$ at $f=f_{I}$ is due to the oscillatory drive and how much is due to the background activity. At the driving frequency $f_{I}$ we write (see Fig. 4c)

[TABLE]

i.e. we measure the contribution of chaotic activity to the power spectrum at the driving frequency by interpolating the power spectrum at neighboring frequencies. The signal-to-noise ratio (SNR) at the driving frequency $f_{I}$ is then given by

[TABLE]

Notice the size of the frequency bin $\Delta f$ scales the SNR, but since we are interested in the dependency of the SNR on $f_{I}$ and not in its numerical value, this scaling factor can be neglected. Finally, we have seen in the example in Fig. 4b that the oscillatory input can suppress background activity at frequencies far from $f_{I}$. In order to quantify this chaos-suppression effect, we split the total variance of $x$ into two contributions (Fig. 4c)

[TABLE]

VI.0.1 Weak drive and signal transmission.

If the oscillatory input is weak, chaos is not entirely suppressed and acts as internally-generated noise on the transmission of the oscillatory input. We now study how the network transmits this oscillatory input signal, and how the transmission quality depends on the signal frequency $f_{I}$. It is known from linear response theory that the transmission of weak signals through single homogeneous populations with strong (intrinsic or external) noise does not benefit from adaptation Schwalger (2013); Deger et al. (2014). This is because in the signal-to-noise ratio (SNR) both the signal and the noise are affected in the same way Deger et al. (2014). We wondered whether in a strongly coupled, large random network, adaptation could have a different effect on the oscillatory signal than on the noise, thereby re-shaping the SNR. A particularly interesting question is how signals are transmitted in the presence of purely intrinsically-generated chaotic fluctuations that are shaped by adaptation and recurrent connectivity.

To understand why adaptation cannot shape the SNR in a weakly-coupled network, consider our random network in the non-chaotic regime, with $g\ll g_{c}$ (Fig. 5a). If we drive the non-chaotic network with oscillatory input together with a noise source $\eta$, the typical response of one unit in the network can be approximated using the mean-field linear frequency response function $\tilde{\chi}_{\beta}(f)$, whose square modulus is given by (see appendix E)

[TABLE]

where we add the subscript $\beta$ to stress that $\tilde{\chi}_{\beta}$ depends on the adaptation parameters $\gamma$ and $\beta$ (cf. Eq. 4). If we indicate by $S_{I}(f)$ and $S_{\eta}(f)$ the power spectral density of the oscillatory input and of the external noise respectively, the power spectral density of the output (taken as the network activity $x$) can be approximated by

[TABLE]

This means that both the signal and the noise are shaped by the same factor $|\tilde{\chi}_{\beta}(f)|^{2}$ that characterizes the network (Fig. 5a). The SNR of the output at the driving frequency, defined as in Eq. 26, is given by

[TABLE]

i.e. the parameters of the network, reflected in the linear response function $\chi_{\beta}$, do not influence the SNR (Fig. 5a). Notice that we considered the activation variable $x$ as our output. We verified that considering instead the firing rate $\phi(x)$ as the output yields qualitatively the same results, therefore we will for simplicity continue our analysis for the output $x$. Eq. 30 implies that the SNR depends only on the power spectra of the signal and of the noise. For example, if we consider low-frequency dominated noise, high-frequency signals will be transmitted more easily, but once again the introduction of adaptation will not play any role (Fig. 5a). While this argument is based on a linear response approximation, we verified using the DMFT solution that the linear approximation is quite accurate. Deviations are visible very close to the criticality, but once again the SNR is almost entirely independent of the neuron parameters.

The findings are completely different for a network in the chaotic phase, i.e. $g>g_{c}$. As discussed above, in this regime the network produces internal fluctuations whose power spectrum depends on single neuron parameters (Fig. 5b, see also section IV). For clarity, let us assume that there is no external noise, such that the noise is only internally generated by the network. In this case, the linear response theory framework cannot be applied; in order to predict the effect of the network in shaping both the input and the internally-generated noise, we need to solve the DMFT equations (Eq. 18) iteratively. As in the previous section, the resulting power spectrum can be split into a chaotic component and into an oscillatory component (see Eq. 24). How does the introduction of adaptation shape these two components? We have seen that in the presence of adaptation the network can enter a state of resonant chaos that differs from the traditional chaos of a network without adaptation because of the presence of a dominant frequency band centered at $f_{0}$ (Fig. 5b, see also section IV). The state of resonant chaos survives in the presence of weak input. As the driving frequency changes, the amplitude of the transmitted signal $A_{\text{osc}}$ passes through a maximum at the resonance frequency $f_{0}$ of the network; however, $A_{\text{osc}}$ decreases with the distance from the resonance frequency more slowly than the noise amplitude $A_{\text{bgk}}$ does (Fig. 5b). This difference is related to the nonlinear sharpening effect of the recurrent network, that can be captured by DMFT. As a consequence, the SNR is maximal at very slow driving frequencies, goes through a minimum at the resonance frequency $f_{0}$ before it increases again (Fig. 5b). We can conclude that in the chaotic regime adaptation improves the SNR at low frequencies, whereas in weakly-coupled, non-chaotic networks such an improvement cannot be observed, independently of the choice of the adaptation parameters $\gamma,\beta$. If the strength of the input is increased, the interaction between noise and signal becomes stronger, leading to a deformation of the SNR (Fig. 6a). However, even for strong drive we observe a peak of the SNR at frequencies that are lower than the resonance one.

VI.0.2 Strong drive and chaos suppression.

In the presence of strong input, chaos suppression together with the formation of a sharp peak are indications that at the microscopic level the network is driven towards a limit cycle. Similarly to Rajan et al. (2010), we now study how chaos suppression depends on the driving frequency $f_{I}$. By solving the DMFT equations (Eq. 18) in the presence of external input, we find that both $P_{\text{chaos}}$ and $P_{\text{osc}}$ exhibit a non-monotonic dependence on $f_{I}$ (Fig. 6b). $A_{\text{osc}}$ depends smoothly on $f_{I}$, reaching its largest value around $f_{0}$. On the other hand, $A_{\text{chaos}}$ is zero for input frequencies that are close to $f_{0}$, indicating that the network is driven into a limit cycle. While a network without adaptation also exhibits such a non-monotonic dependence Rajan et al. (2010), in our case this effect is more pronounced due to the resonant power spectrum of the spontaneous activity in the presence of adaptation.

VII Discussion

We studied how the dynamics of a random network of rate neurons are shaped by the properties of single neurons, and in particular by the presence of history-dependent mechanisms such as adaptation. To this end, we generalized DMFT, a well-established theoretical tool Sompolinsky et al. (1988), to the case of multi-dimensional rate units. This allowed us to reduce the high-dimensional, deterministic network model to a low dimensional system of stochastic differential equations. Standard approaches to solving the mean-field equations Sompolinsky et al. (1988) were not fruitful in the multi-dimensional setting. However, the mean-field solution could be found efficiently in a semi-analytical way using an iterative approach. The iterative approach highlights how recurrent connections sharpen the response function of single neurons, i.e. how bands of preferred frequencies become narrower (see also appendix B). Previous studies that considered the role of single neuron properties on random network dynamics focused only on the role of the gain function Kadmon and Sompolinsky (2015); Mastrogiuseppe and Ostojic (2017). To our knowledge, this is the first result that relates the single neuron frequency response function to the spectral properties of random network dynamics.

We studied in detail the important case of neuronal adaptation, using a two-dimensional rate model. We showed that adaptation extends the stability region of a recurrent network of rate units because the transition from a stable fixed point to a fluctuating regime happens at $g=g_{c}>1$, i.e. for higher coupling strength than for the network without adaptation. Crucially, above the criticality and for slow adaptation, the dynamics settle in a state of “resonant chaos” that, unlike the chaotic activity of networks of rate units without adaptation, is dominated by a nonzero resonance frequency. We observed that the resonance frequency can be computed from the single unit properties and it is therefore independent of the connection strength $g$. On the other hand, the presence of recurrent connections increases the coherence of the oscillations and therefore influences the correlation time. The oscillation coherence is maximal at the onset of chaos and decreases with $g$, for $g>g_{c}(\gamma,\beta)$. Indeed, as it is typical of critical behavior, the correlation time in the chaotic phase diverges when approaching the criticality. In the presence of adaptation, this happens because the system approaches a limit cycle.

It is interesting to observe that for slow adaptation there are two separate contributions to the correlation time of the network activity: an oscillatory component, related to the resonance frequency, and a long tail that scales with the adaptation timescale. For finite $\tau_{a}$, the correlation time diverges when $g\rightarrow g_{c}$ due to the oscillatory component. For $\tau_{a}\rightarrow\infty$, the correlation time also diverges, but this is due to the long tail, since both the resonance frequency and the Q-factor go to zero for large $\tau_{a}$, yielding a finite and therefore sub-dominant contribution to the correlation time. Such multi-scale structure of the autocorrelation could be advantageous for network computations that require expressive dynamics over multiple timescales, as it is often the case in motor control. Indeed, adaptation has been proposed to play a role in sequential memory retrieval Deco and Rolls (2005), slow activity propagation Setareh et al. (2018), perceptual bistability Shpiro et al. (2009) and decision making Theodoni et al. (2011). Moreover, SFA has beneficial consequences both for reservoir computing approaches Nicola and Clopath (2017) and for spiking neuron-based machine learning architectures Bellec et al. (2018). Further work could explore the relation between long correlation time induced by adaptation and computational properties.

In the presence of strong oscillatory input, chaos can be suppressed Schwalger et al. (2006); Rajan et al. (2010). In particular, in the presence of adaptation, chaos is more easily suppressed when the driving frequency is close to the resonance frequency. In the presence of weak input, chaos is not fully suppressed. Interestingly, we found that in the chaotic regime the presence of adaptation shapes the SNR in frequency space. In particular, adaptation increases the SNR for low-frequency signals, a possibly important feature since behaviorally relevant stimuli can have information encoded in slow signal components Melamed et al. (2004). Crucially, this effect is not present in the sub-critical regime ($g<g_{c}$), since signal and external noise are shaped together Deger et al. (2014). It is known that the properties of biological neurons, including adaptation parameters, can be dynamically adjusted using neuromodulators McCormick and Williamson (1989); Stiefel et al. (2009). In view of our results, this would allow to dynamically shape the SNR depending on the requirements imposed by the behavioral context.

While our theory is applicable to single units with $D$ interacting variables, the effect of a single adaptation variable ($D=2$) on the dynamics of random recurrent networks was also studied independently and simultaneously by another group Beiran and Ostojic (2018), who reached results consistent with ours Muscinelli et al. (2018). The authors of Beiran and Ostojic (2018) used a slightly different network architecture and did not focus on the relation between single neuron response and spectral properties, but rather on the correlation time of the network activity and on the effect of white noise input. One major difference is the conclusion reached regarding correlation time: by using a different definition, in Beiran and Ostojic (2018) the authors conclude that the correlation time does not scale with the adaptation timescale. Based on our analysis, we infer that the definition of correlation time used in Beiran and Ostojic (2018) captures only the oscillatory contribution to the correlation time, and not its long tail.

Current mean-field theories for spiking neural networks Brunel (2000) are self-consistent only with respect to mean activities (firing rates), whereas second-order statistics such as autocorrelation function or power spectral density of inputs and outputs are inconsistent Wieland et al. (2015). While iterative numerical procedures are available Lerchner et al. (2006); Dummer et al. (2014); Wieland et al. (2015), a self-consistent analytical calculation of the autocorrelation (or power spectrum) via DMFT for networks of spiking neurons is known to be a hard theoretical problem. In the present manuscript, the rate-based modeling framework allowed us to put forward explicit expressions for the map of autocorrelations. For a general nonlinearity $\phi(x)$, this map takes the form of an infinite series (according to Eq. 36 in appendix B). However, for polynomial nonlinearities the series simplifies to a finite sum, e.g. Eq. 71 in appendix F, which permits a closed-form analytical expression for the iterative map. Therefore, our study offers a unique method for the calculation of the autocorrelation in biologically constrained random neural networks, and thus represents a promising step towards a self-consistent mean-field theory beyond first-order rate models Sompolinsky et al. (1988); Kadmon and Sompolinsky (2015).

Extensions and generalizations

We see four extensions to the work presented in the present manuscript. First, our study is limited to rate neurons while it would be interesting to extend the analysis to spiking neuron models. As a first step in this direction, previous work has already investigated the introduction of white noise in random rate networks Kadmon and Sompolinsky (2015); Schuecker et al. (2018); Beiran and Ostojic (2018), which would be straightforward to include in the case of $D$-dimensional rate units. Second, our framework can readily be extended to multiple adaptation variables (see Fig. 8 for two examples). This is a key feature in order to account for realistic SFA, which is known to have multiple timescales and it has been shown to have power-law structure Fairhall et al. (2001); La Camera et al. (2006); Lundstrom et al. (2008); Pozzorini et al. (2013). Interestingly, our framework can be extended to power-law adaptation, since we require only the knowledge of the linear frequency-response function of the single neurons. We expect that in this situation the internal noise generated by the network will also have a power-law profile of the type $f^{\alpha}$, with $\alpha>0$. With such a noise spectrum, the signal that maximizes information transmission should be dominated by low-frequencies in a power-law fashion Rieke et al. (1996); Pozzorini et al. (2013). Third, the introduction of additional structure in the connectivity, such as low-rank perturbations Mastrogiuseppe and Ostojic (2018), attractor structure Pereira and Brunel (2018), or large scale connectivity of the brain Deco et al. (2011), could give rise to interesting dynamics when combined with single units with multiple adaptation variables. Finally, while our study focused on neural networks, random network models are used in other areas of biology and physics Dorogovtsev and Mendes (2002). By extending mean-field theory techniques to more complex node dynamics, our approach also contributes to understanding the interaction between node dynamics and network structure in more general settings.

Acknowledgments

The authors would like to thank Johanni Brea for helpful discussions and comments. This research was supported by the Swiss national science foundation, grant agreement 200020_165538.

Appendix A Stability of the fixed point

General theory.

The system of $N\cdot D$ coupled nonlinear differential equations (Eq. 2) becomes intractable for large $N$. However, because $\phi(0)=0$, the system has a fixed point at the origin $\{x_{i}^{\alpha}=0\}_{i=1,\dotsc,N}^{\alpha=1,\dotsc,D}$, the stability of which can be studied owing to the clustered structure of the system. The Jacobian at the fixed point is given by

[TABLE]

where $\mathrm{J}$ is the random connectivity matrix and $\mathrm{I}_{N}$ is the $N$-dimensional identity matrix. The matrix $\mathrm{B}$ is of size $ND\times ND$ and it therefore has $ND$ eigenvalues. Since all the blocks of $\mathrm{B}$ commute with each other, we can apply the result of Silvester (2000) to find a relation between the eigenvalues of $\mathrm{J}$, $\mathrm{A}$ and $\mathrm{B}$

[TABLE]

where $\mathrm{A}^{-}$ is the matrix obtained by removing the first column and the first row from the matrix $\mathrm{A}$. This expression is valid for all the eigenvalues of $\mathrm{B}$ that are not coincident with those of $\mathrm{A}^{-}$. Eq. 32 can be transformed into a polynomial equation of degree $D$ in $\lambda_{\mathrm{B}}$, so that for every value of $\lambda_{\mathrm{J}}$ we obtain $D$ eigenvalues of $\mathrm{B}$, as expected. From now on we will assume that, without loss of generality, $\phi^{\prime}(0)=1$.

In the $N\rightarrow\infty$ limit, the eigenvalues $\lambda_{\mathrm{J}}$ are known to be uniformly distributed on a disk in the complex plane, centered at zero and of radius $g$ Girko (1985). If one can invert Eq. 32, it becomes computationally fast to compute the eigenvalues of the Jacobian in the $N\rightarrow\infty$ limit without finite-size effects. Whether one can obtain an explicit inverse formula depends on the dimensionality and on the entries of the matrix $\mathrm{A}$.

Network with adaptation.

For the two-dimensional model defined by Eqs. (3, 4), we can invert Eq. 32, and obtain an expression for the eigenvalues of the Jacobian

[TABLE]

Using the mapping between the eigenvalues of the connectivity matrix $\mathrm{J}$ and those of the Jacobian matrix $\mathrm{B}$ ( Eq. 33), we find the critical value of $g$ for which the stability of the fixed point is lost

[TABLE]

where $\beta_{H}(\gamma)=-1-\gamma+\sqrt{2\gamma^{2}+2\gamma+1}$. The critical value $g_{c}$ can also be calculated from dynamical mean-field theory (see appendix E).

The bifurcation that characterizes the loss of stability depends on two parameters, viz. the ratio of timescales $\gamma$ and the strength of the adaptation $\beta$. To further characterize the bifurcation at $g=g_{c}(\gamma,\beta)$, we can study the imaginary part of the critical eigenvalue, i.e. the one with real part equal to zero at $g=g_{c}(\gamma,\beta)$. If the adaptation strength $\beta$ has a value $\beta\leq\beta_{H}(\gamma)$, then the imaginary part of the critical eigenvalue is equal to zero corresponding to a saddle-node bifurcation at $g=g_{c}(\gamma,\beta)$. On the other hand, if $\beta>\beta_{H}(\gamma)$, then the critical eigenvalue is a pair of complex-conjugate, purely imaginary eigenvalues, a signature of a Hopf bifurcation. Therefore, we introduce the curve $\beta=\beta_{H}(\gamma)$, which separates the positive quadrant of the $\gamma-\beta$ plane in two regions: one in which the system becomes unstable at the critical value $g_{c}(\gamma,\beta)$ via a saddle-node bifurcation, and another one in which the instability occurs via a Hopf bifurcation (Fig 7a). In the Hopf-bifurcation region, the imaginary part of the critical eigenvalues can be computed analytically:

[TABLE]

The parameter $f_{m}$ is the frequency of low-amplitude oscillations close to the bifurcation, if $N<\infty$. In the finite-$N$ case, we find numerically that these low-amplitude oscillations are stable. When $N\rightarrow\infty$, however, we find that chaotic dynamics onset right above the bifurcation (see section IV). The frequency $f_{m}$ is monotonic in $\beta$ but non-monotonic in $\gamma$ (Fig. 7b), indicating that a slower adaptation variable (smaller $\gamma$) does not necessarily correspond to slower oscillations. When considering codimension-two bifurcations, we have that for $g=g_{c}$ and $\beta=\beta_{H}(\gamma)$ the system undergoes a Bogdanov-Takens bifurcation.

Appendix B Qualitative study of the iterative map

For a qualitative understanding of the effect of the iterations on the power spectral density, we exploit the fact that $x^{1}$ is a Gaussian process, for which the following formula holds Malakhov (1978)

[TABLE]

where the angular brackets indicate the mean over the statistics of $x^{1}$. Eq. 36 gives the effect of a nonlinearity $\phi$ on a the autocorrelation of a Gaussian process $x^{1}$. By truncating the series after the first term, we get

[TABLE]

Fourier transforming this equation we get an approximation of the power spectral density of $\phi(x^{1})$

[TABLE]

where we we introduced the function $\Psi_{1}\left(\int_{-\infty}^{\infty}S_{x}(f^{\prime})df^{\prime}\right):=\left(\left\langle\phi^{\prime}(x^{1})\right\rangle\right)^{2}$ to highlight the fact that the coefficient that multiplies $S_{x}(f)$ depends on the area under the power spectral density, i.e. on the variance of $x^{1}$, and is therefore nonlocal in frequency space. We stress that retaining only the first term in Eq. 36 is different than considering a linear approximation of $\phi$, since the dependence of the coefficient on the variance would not appear in that case.

Using this approximation, we can express the power spectral density at the $n^{\text{th}}$ iteration of the iterative method, as a function of the initial power spectral density $S_{\phi(x^{1})}^{(0)}(f)$ from which we started to iterate. We obtain

[TABLE]

where $\Psi_{1}^{(n)}:=\Psi_{1}\left(\int_{-\infty}^{\infty}(S_{x})^{(n)}(f^{\prime})df^{\prime}\right)$. If we take $S_{\phi(x^{1})}^{(0)}(f)$ to be constant and we define $a_{n}=\left(\prod_{k=1}^{n-1}\Psi_{1}^{(k)}\right)$, we can rewrite the above expression as

[TABLE]

If $g>g_{c}$, there will be a range of frequencies for which $g^{2}\tilde{G}(f)>1$, which implies that its $n^{\text{th}}$ power diverges when $n$ grows. In a purely linear network, this phenomenon would lead to a blow-up of the power spectral density, in agreement with the fact that activity in a linear network is unbounded for $g>g_{c}$. If $\phi$ is a compressive nonlinearity however, the coefficient $a_{n}$ will tend to zero for growing $n$, counterbalancing the unbounded growth of $\left(g^{2}\tilde{G}(f)\right)^{n}$. Based on Eq. 40, we would predict that all the modes for which $\tilde{G}(f)>1/g^{2}$ will get amplified over multiple iterations, while all the other modes will get suppressed. While this is a highly simplified description, the suppression and the amplification of modes is clearly visible when observing the evolution of the power spectrum over iterations (Fig. 2f) and when comparing the dynamics of the self-consistent solution (Fig. 8c,f) to the corresponding linear response function (Fig. 8b,e). When truncating the series after the first order however, the mean-field network does not admit a self-consistent solution, for which we need to retain also higher order terms. Such terms will balance the progressive sharpening of the power spectrum, allowing for a self-consistent solution.

As an example of higher-order term, consider the next term in the series in Eq. 36, given by

[TABLE]

where $\Psi_{2}$ is defined analogously to $\Psi_{1}$. In general, higher-order terms will contain convolutions of the power spectral density with itself, which are responsible for the creation of higher harmonics. To qualitatively understand this effect, consider the case in which $S_{x}(f)$ is a Dirac $\delta$-function with support in $f_{0}$. In this case, the two-fold convolution of $S_{x}(f)$ with itself is again equal to a Dirac $\delta$-function, but centered in $2f_{0}$. A similar argument can be given for resonant power spectral densities, which implies that a self-consistent solution should exhibit harmonics of the fundamental resonance frequency. Note that in this paper we considered odd functions, for which only odd terms in the series are nonzero.

For higher values of $g$, the relative importance of higher-order terms in the series in Eq. 36 will increase, leading to a broader power spectrum. The self-consistent power spectrum however, seems to be always narrower than the single neuron linear response function. For a possible explanation of this phenomenon, we consider the $g\rightarrow\infty$ limit, which was already studied in Crisanti and Sompolinsky (2018) for the network without adaptation. Using the same technique, we conclude that in this limit the autocorrelation decay tends to be the same as one obtained for a single unit driven by white noise Crisanti and Sompolinsky (2018). In the frequency domain, this is equivalent to say that the power spectral density of the network tends to the one of a single unit driven by white noise.

Appendix C Mean-field theory derivation

In this section, we extend the derivation of dynamical mean-field theory (DMFT) to the case of the network of multi-dimensional rate units. Since there are no additional complication with respect to the standard case, we report here only the main steps. For a review of the path-integral approach to DMFT, see e.g. Schücker et al. (2016); Crisanti and Sompolinsky (2018). The moment-generating functional corresponding to the microscopic system in Eq. 2 is given by

[TABLE]

where

[TABLE]

and we introduced the notation $\tilde{\boldsymbol{x}}^{T}\boldsymbol{x}=\sum_{\alpha}\sum_{i}\int{\tilde{x}_{i}^{\alpha}(t)x_{i}^{\alpha}(t)dt}$. The integral is over paths and bold symbols indicate vectors, over both the network space and the rate model space, so that \mathcal{D}\boldsymbol{x}\mathrel{\hbox to0.0pt{\raisebox{1.29167pt}{\cdot}\hss}\raisebox{-1.29167pt}{\cdot}}=\prod_{\alpha}\prod_{i}\mathcal{D}x_{i}^{\alpha}.

We are interested in properties that are independent of the particular realization of the coupling matrix $\mathrm{J}$. In order to extract those properties, we average over the quenched disorder by defining the averaged generating function

[TABLE]

The average over each $J_{ij}$ can be computed by noticing that the terms corresponding to different $J_{ij}$ factorize and the integral can be solved by completing the square. Since the details of this calculation are analogous to the one-dimensional case, we directly report the result

[TABLE]

We now aim to decouple the interaction term in the last line by introducing the auxiliary field

[TABLE]

We introduce $Q_{1}$ in the generating functional by inserting the following representation of the unity

[TABLE]

where $\delta[\cdot]$ is the delta functional. Using the integral representation of the delta functional leads to the introduction of a second auxiliary field, which we call $Q_{2}$. We obtain

[TABLE]

This expression has the advantage that any interaction between different units is removed and all the contribution coming from different units factorize. It is convenient to rewrite the averaged generating functional as a field theory for two auxiliary fields $Q_{1},Q_{2}$, i.e. we remove the vectorial response terms $\boldsymbol{j}^{T}\boldsymbol{x},\tilde{\boldsymbol{j}}^{T}\tilde{\boldsymbol{x}}$ and we add two scalar response terms for the auxiliary fields. The result is

[TABLE]

where we extended our notation to Q_{1}^{T}Q_{2}\mathrel{\hbox to0.0pt{\raisebox{1.29167pt}{\cdot}\hss}\raisebox{-1.29167pt}{\cdot}}=\int\int Q_{1}(s,t)Q_{2}(s,t)dsdt. The crucial observation to make is that essentially all factors associated to different units factorized yielding the factor $N$. For this reason, the integration is now only over all rate model indices but over only one unit index. The remainder is the problem of one unit, characterized by $D$ variables, interacting with two external fields $Q_{1},Q_{2}$.

The final step is to perform a saddle-point approximation, i.e. replace $Q_{1},Q_{2}$ by their values that make the action stationary. To do this, we need to solve the two saddle-point equations

[TABLE]

These equations are analogous to the ones in the one-dimensional case, and lead to the saddle-point solution

[TABLE]

where $C_{\phi(x^{1})}(s,t)$ is the autocorrelation function of $\phi(x^{1})$ evaluated at the saddle point solution. The averaged generating functional at the leading order in $N$ can be written as

[TABLE]

This is the statistical field theory corresponding to $D$ linearly interacting variables, with $x^{1}$ that receives a Gaussian noise whose autocorrelation is given by $C_{\phi(x^{1})}$. Writing the corresponding differential equations results in our mean-field description (Eq. 8).

Appendix D Mean-field theory with heterogeneous adaptation

In this section, we will extend the derivation of the dynamic mean-field theory (DMFT) for the case of the network with heterogeneous adaptation. We consider the case in which each neuron has different parameters, sampled i.i.d from the same distributions, and different parameters of the same neuron are uncorrelated with each other. More precisely, we sample the elements of the matrix $\mathrm{A}_{i}$ for neuron $i$ as

[TABLE]

where the subscript $i$ runs over the neurons in the network.

In deriving the mean-field theory, most of the steps are identical to those in appendix C, so we will focus on the additional terms due to the new source of disorder. We separate the contribution of mean adaptation parameters $\bar{\mathrm{A}}^{\alpha\beta}$ from the deviations, so that the generating functional reads

[TABLE]

where

[TABLE]

and $\bar{\mathrm{A}}$ is the matrix of the expected values of $\mathrm{A}$.

The action $S_{0}$ is the same as for the network without heterogeneity, and when averaging over the connectivity disorder, we obtain the same result as for homogeneous network. In this case however, we need to also average over the disorder due to heterogeneity, i.e. over all the $\mathrm{A}_{k}^{\alpha\beta}$. The averaged generating functional will then result from the average

[TABLE]

The new terms due to the heterogeneity result in integrations of the type

[TABLE]

that can be solved by completing the square. After averaging over both the connectivity disorder and the heterogeneity disorder, the generating functional reads

[TABLE]

The last term, which is due to the heterogeneity, factorizes into the contributions associated to different units. From this point on, in order to derive the mean-field equations, we follow exactly the same steps as in appendix C, so we do not report those steps here. The mean-field equations read

[TABLE]

where $\eta_{H}^{\alpha\beta}$ are Gaussian processes associated to the heterogeneity, that all have mean zero and autocorrelation

[TABLE]

For the particular case of adaptation with heterogeneity on the parameter $\beta$, as studied in section IV, we have the following mean-field equations

[TABLE]

where $\eta_{H}(t)$ is a Gaussian process with mean zero and autocorrelation

[TABLE]

From Eqs. 60,61, we can find the self-consistent equation for the power spectrum:

[TABLE]

where $\tilde{G}_{H}(f)$ is an effective filter given by

[TABLE]

where $\omega=2\pi f$. The effective filter $\tilde{G}_{H}(f)$ predicts a larger power at low frequencies, similar to what is observed in simulations (cf. Fig. 2d).

Appendix E Fixed point stability in the mean-field network

Here we consider the full matrix of linear response functions (see below), to conclude that the only quantity that matters for the stability at the fixed point is $\tilde{G}(f)$.

Starting from the microscopic network equations (Eq. 2), we derive a set of differential equations, that we write in matrix form

[TABLE]

where $\Delta_{1}=\delta^{\alpha 1}\delta^{\beta 1}$ is a matrix whose only nonzero element is $[\Delta_{1}]^{11}=1$. $\chi_{ik}(\tau)$ is a $D$ by $D$ matrix, whose component are defined as $\chi_{ik}^{\alpha\beta}(\tau)=\frac{\delta x_{i}^{\alpha}(\tau)}{\delta h_{k}^{\beta}(0)}$, where $h_{k}^{\beta}$ is a small perturbation given to the variable $x_{k}^{\beta}$ at time $\tau=0$. Notice that in deriving Eq. 65, we have assumed stationarity and that $\phi^{\prime}(0)=1$. We now Fourier transform Eq. 65 and get

[TABLE]

Inverting the matrix $\left(2\pi if\,\mathrm{I}_{D}-\mathrm{A}\right)$ and recognizing the linear response function of the single unit $\tilde{\chi}_{0}(f)$, we obtain

[TABLE]

where $\tilde{\chi}_{0}(f)$ is a $D$ by $D$ matrix whose elements are $\tilde{\chi}_{0}^{\alpha\beta}(f)$, defined in section III.

Since in the mean-field approximation the mean of the linear response function is zero, we look for the second moments Kadmon and Sompolinsky (2015). We multiply every element of the matrix equation (Eq. 67) by its complex conjugate and average over the quenched disorder. We obtain

[TABLE]

where the absolute value is intended element-wise. Due to the structure of the matrix $\Delta_{1}$, we have that $|\tilde{\chi}_{0}(f)\Delta_{1}\tilde{\chi}(f)|^{2}=\tilde{G}(f)\Delta_{1}|\tilde{\chi}(f)|^{2}$, as it can be verified simply by using the definition of $\Delta_{1}$. Finally, we can solve for $|\tilde{\chi}(f)|^{2}$

[TABLE]

Since the only nonzero eigenvalue of the matrix $\tilde{G}(f)\Delta_{1}$ is $|\tilde{\chi}_{0}^{11}(f)|^{2}$, the stability condition for the fixed point is given by

[TABLE]

Appendix F Effect of nonlinearities on second-order statistics

In this section, we provide some additional details on how to compute the effect of nonlinearities on the second order statistics (autocorrelation or power spectral density) of a Gaussian process. We consider three cases of interest: polynomials, piecewise linear functions and arbitrary nonlinear functions. To simplify our notation, we drop the superscript of and consider a generic Gaussian process $x^{1}$.

The effect of polynomial nonlinearities can be expressed in closed form in time domain. This can be seen by considering again the infinite series expression (Eq. 36), valid for stationary Gaussian processes $x$

[TABLE]

where the angular brackets indicate the average over the statistics of $x$. In the case in which $\phi$ is a polynomial of degree $p$, only the terms in the sum up to $p$ are nonzero. As an example, we can compute the effect of a cubic approximation of the hyperbolic tangent, i.e. $\phi(x)\simeq\phi_{3}(x):=x-\frac{x^{3}}{3}$

[TABLE]

As expected, the effect of the nonlinearity depends on $C_{x}(0)$ i.e. on the variance of $x$ itself. Notice that the coefficient of the first term is compressive (i.e. smaller than one) only if $C_{x}(0)$ is smaller than one itself. This type of behavior is expected since $\phi_{3}$ is unbounded.

Another interesting case are piecewise linear nonlinearities. In this case, we use Price’s theorem twice to get

[TABLE]

For a piecewise linear $\phi$, the second derivative $\phi^{\prime\prime}$ is a sum of Dirac’s delta functions with variable coefficients. More precisely, we consider

[TABLE]

where $x_{p}$ are the points in which the first derivative is discontinuous, $c_{p}$ are some arbitrary coefficients and $\Theta(\cdot)$ is the Heaviside function. The second derivative of $\phi_{PL}$ is given by

[TABLE]

The delta functions allow us to compute the correlation function $C_{\phi_{PL}^{\prime\prime}}(t)$ explicitly

[TABLE]

where we defined $\rho(t):=\frac{C_{x}(t)}{C_{x}(0)}$. Inserting Eq. 76 in Eq. 73 and integrating twice with respect to $C_{x}(t)$ we get

[TABLE]

In the case in which $\phi$ is an odd function, the term $f_{\phi}\left(0;C_{x}(0)\right)$ is equal to zero. For the specific case of the piecewise linear approximation of the hyperbolic tangent considered in this paper, i.e.

[TABLE]

the expression in Eq. 77 reduces to

[TABLE]

For the piecewise linear function, an alternative approach is based on the infinite series in Eq. 36, which yields Stratonovich (1967); Kruscha and Lindner (2016):

[TABLE]

with input variance $\sigma^{2}=C_{x}(0)$ and cumulative Gaussian distribution function $F(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-y^{2}/2}\,dy$. For the figures in this paper, we used the map in Eq. F.

For an arbitrary nonlinear function, we can use two methods. The first method is a semi-analytical approach that relies on the integral form of the autocorrelation of the rate $C_{\phi(x)}(\tau)$ as a functional of the autocorrelation $C_{x}(\tau)$ of $x$ Schücker et al. (2016)

[TABLE]

where $Dx=e^{-x^{2}/2}dx$. Notice that a slightly different version of this formula was already proposed in Sompolinsky et al. (1988). Therefore, to obtain the effect of $\phi$ on the power spectral density, one should 1) inverse Fourier transform $S_{x}(f)$ to get $C_{x}(\tau)$ 2) apply Eq.(81), by computing the two integrals numerically 3) Fourier transform $C_{\phi(x)}(\tau)$ to get $S_{\phi(x)}(f)$. Practically, this procedure requires the application of the fast Fourier transform algorithm and the numerical evaluation of two integrals.

The second method is fully numerical and it can be useful in cases in which the integrals in the first method are expensive to evaluate numerically. This method consists in approximating the power spectral density $S_{\phi(x)}$ via Monte Carlo sampling. More precisely, we sample multiple realizations in frequency domain of the Gaussian process with zero mean and power spectral density $S_{x}(f)$. We then transform each sample to time domain and apply the nonlinearity $\phi(x)$ to each sample $x(t)$ individually. Finally, we transform back to Fourier domain and get $S_{\phi(x)}$ by averaging. Despite being computationally more expensive than the closed form expressions, this sampling method provides a solution of the mean-field theory for an arbitrary nonlinearity and it is computationally much cheaper than running the full microscopic simulation. Moreover, this method can easily be extended to be used in the presence of a non-Gaussian sinusoidal input (cf. section VI and Rajan et al. (2010)).