{'fraction_non_alphanumeric': 0.04084224135940711, 'fraction_numerical': 0.013898829450742271, 'mean_word_length': 4.528976257339801, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A reservoir of timescales in random neural networks

The temporal activity of many biological systems, including neural circuits, exhibits fluctuations simultaneously varying over a large range of timescales. The mechanisms leading to this temporal heterogeneity are yet unknown. Here we show that random neural networks endowed with a distribution of self-couplings, representing functional neural clusters of different sizes, generate multiple timescales activity spanning several orders of magnitude. When driven by a time-dependent broadband input, slow and fast neural clusters preferentially entrain slow and fast spectral components of the input, respectively, suggesting a potential mechanism for spectral demixing in cortical circuits.

The temporal activity of many biological systems, including neural circuits, exhibits fluctuations simultaneously varying over a large range of timescales. The mechanisms leading to this temporal heterogeneity are yet unknown. Here we show that random neural networks endowed with a distribution of self-couplings, representing functional neural clusters of different sizes, generate multiple timescales of activity spanning several orders of magnitude. When driven by a time-dependent broadband input, slow and fast neural clusters preferentially entrain slow and fast spectral components of the input, respectively, suggesting a potential mechanism for spectral demixing in cortical circuits.

I. INTRODUCTION
Experimental evidence shows that the temporal activity of many physical and biological systems exhibits fluctuations simultaneously varying over a large range of timescales. In condensed matter physics for example, spin glasses typically exhibit aging and relaxation effects whose timescales span several orders of magnitude [1]. In biological systems, metabolic networks of E. coli generate fluxes with a power-law distribution of rates [2, 3]. Gas release in yeast cultures exhibit frequency distributions spanning many orders of magnitude [4], endowing them with robust and flexible responses to the environment [5].
In the mammalian brain, a hierarchy of timescales in the activity of single neurons is observed across different cortical areas from occipital to frontal regions [6][7][8]. Moreover, neurons within the same local circuit exhibit a large range of timescales from milliseconds to minutes [9][10][11]. This heterogeneity of neuronal timescales was observed in awake animals during periods of ongoing activity, in the absence of external stimuli or behavioral tasks, suggesting that multiple timescales of neural activity may be an intrinsic property of recurrent cortical circuits. Recent studies highlighted the benefits of leveraging computations on multiple timescales when performing complex tasks in primates [12] as well as in artificial neural networks [13]. However, the neural mechanisms underlying the emergence of multiple timescales are not yet understood. We suggest here such mechanism.
We focus on random neuronal networks whose units are recurrently connected, with couplings that are chosen randomly. In our model, each network unit represents a functional cluster of cortical neurons with similar response properties. We interpret the unit's selfcoupling as the neural cluster strength, reflecting the product of the cluster size and the average value of the recurrent synaptic coupling between its neurons. In the case where the self-couplings are zero or weak (order 1/ √ N ), random networks are known to undergo a phase transition from silence to chaos when the variance of the random couplings exceeds a critical value [14]. When the self-couplings are strong (order 1) and are all equal, a third phase appears featuring multiple stable fixed points accompanied by long transient activity [15]. In all these cases, all network units exhibit the same intrinsic timescale, estimated from their autocorrelation function. Here, we demonstrate a novel class of recurrent networks, capable of generating temporally heterogeneous activity whose multiple timescales span several orders of magnitude. We show that when the self-couplings are heterogeneous, a reservoir of multiple timescales emerges, where each unit's intrinsic timescale depends both on its own self-coupling and the network's self-coupling distribution. In particular, we find an exponential relationship between (a power of the) unit's self-coupling and its timescales. We analytically study the dynamics of a single unit in the limit of large selfcoupling, revealing a new metastable regime described by colored noise-driven transitions between potential wells. We show that these results generalize to biologically plausible models of cortical circuits based on spiking networks with cell-type specific clustered architectures, where a reservoir of timescales emerges in presence of a heterogeneous distribution of cluster strengths. We then study the stimulus-response properties of our networks with heterogeneous self-couplings. In networks with zero or weak self-couplings chaotic activity is suppressed best at a single resonant frequency [16]. However, when we drive our networks with a time-dependent broadband input featuring a superposition of multiple frequencies, we find that the chaotic activity is suppressed across multiple frequencies which depend on the units' respective selfcouplings. We see that units with large and small selfcouplings are preferentially entrained by the low and high frequency components of the input, respectively. This spectral specificity suggests that a reservoir of timescales may be a natural mechanism for cortical circuits to flexibly demix different spectral features of complex time-  4), for sα > 1. Stable solutions are allowed within the dark green region. b-i) The shape of a unit's contribution to stability q −1 = (sα − cosh xα 2 ) −2 , from Eq. (7). Stable solutions of xα − sα tanh xα = η, filled blue circles in (a-i), with different |x| values contribute differently to stability. At the edge of chaos only a fixed point configuration with all units contributing most to stability (minimal q −1 ) is stable, light green region in (a-i). a-ii) The curve xα − sα tanh xα for sα < 1. a-iii) A possible distribution of the Gaussian mean-field η. A representative fixed point solution is illustrated by the dashed blue line: for sα < 1 a single solution exists for all values of η, (filled blue circle in a-ii);For sα > 1 multiple solutions exist (a-i) for some values of η; some of them lead to instability (empty blue circle in a-i). The other two solutions may lead to stability (filled blue circles in a-ii), although only one of them will remain stable at the edge of chaos (encircled with green line in a-i).

II. RANDOM NETWORKS WITH HETEROGENEOUS SELF-COUPLINGS
We consider a recurrent network of N rate units obeying the dynamical equations where the random couplings J ij from unit j to unit i are drawn independently from a Gaussian distribution with mean 0 and variance 1/N ; g represents the network gain and we chose a transfer function φ(x) ≡ tanh(x). We measure time in units of 1 ms. The self-couplings s i are drawn from a distribution P (s). The special case of equal self-couplings (s i = s) was studied in [15] and a summary of the results are brought in the Appendix A for convenience. Here we study the network properties in relation to both discrete and continuous distributions P (s).
Using standard methods of statistical field theory [17,18], in the limit of large N we can average over realizations of the disordered couplings J ij to derive a set of self-consistent dynamic mean field equations for each population of units x α with self-coupling strengths (2) In our notation, S denotes the set of different values of self-couplings s α , indexed by α ∈ A, and we denote by N α the number of units with the same self-coupling s α , and accordingly by n α = N α /N their fraction. The mean field η(t) is the same for all units and has zero mean η(t) = 0 and autocorrelation where · denotes an average over the mean field.

III. STABLE FIXED-POINTS AND TRANSITION TO CHAOS
Networks with heterogeneous self-couplings exhibit a complex landscape of fixed points x * α , obtained as the self-consistent solutions to the static version of Eq. (2) and Eq. (3), namely where the mean field η has zero mean and its variance is given by The solution for each unit depends on its respective s α (Fig. 1). If s α < 1 a single interval around zero is available. For s α > 1, for a range of values of η, x * α can take values in one of three possible intervals. However, the available solutions in the latter case are further restricted by stability conditions. We can derive the stability condition by expanding the dynamical equations (1) around the fixed point and requiring that all eigenvalues of the corresponding stability matrix are negative. To determine the onset of instability we look for conditions such that at least one eigenvalue develops a positive real part. An eigenvalue of the stability matrix exists at a point z in the complex plane if [15,19] g 2 α∈A n α 1 − tanh 2 (x α ) 2 z + 1 − s α 1 − tanh 2 (x α ) 2 > 1. (6) Since the denominator of the expression above is z plus the slope of the curve in Fig. 1a-i, a solution whose value x * α gives a negative slope (available when s α > 1) leads to a vanishing value of the denominator at some positive z and hence to a positive eigenvalue and instability. Hence, the n α fraction of units with s α > 1 at a stable fixed point are restricted to have support on two disjoint . We refer to this regime as multi-modal, a direct generalization of the stable fixed points regime found in [15] for a single self-coupling s > 1, characterized by transient dynamics leading to an exponentially large number of stable fixed points. For the n α portion of units with s α < 1, the stable fixed point is supported by a single interval around zero.
A fixed point solution becomes unstable as soon as an eigenvalue occurs at z = 0, obtaining from Eq. (6) the stability condition where q α = s α − cosh 2 (x α ) 2 . For s α > 1 the two possible consistent solutions to (4) that may result in a stable fixed point (from the two disjoint intervals in Fig. 1a-i), contribute differently to q α . Larger |x * α | decreases q −1 α ( Fig. 1b-i), thus improving stability. Choices for distributions of x * α along the two intervals become more restricted as g increases or s α decreases, since both render higher values for the stability condition, Eq. 7, forcing more solutions of x i to decrease q −1 α . This restricts a larger fraction of x * α at the fixed points to the one solution with higher absolute value. At the transition to chaos, a single last and most stable solution exists with all x i values chosen with their higher absolute value x * α ( Fig.  1a-i, light green segments). For those with s α < 1 only one solution is available, obtained by the distribution of η through consistency (4) at the fixed point. In this configuration, the most stable solution is exactly transitioning from stability to instability where (7) reaches unity. Hence the transition from stable fixed points to chaos occurs for a choice of g and P (s) such that solving consistently (4) and (5) leads to saturate the stability condition (7) at one.
We illustrate the discussion above in the case of a network with two sub-populations with n 1 and n 2 = 1 − n 1 portions of the units with self-couplings s 1 and s 2 , respectively. In the (s 1 , s 2 ) plane, this model gives rise to a phase portrait with a single chaotic region separating four disconnected stable fixed-point regions (Fig. 2a). A unit's activity is determined by its own self-coupling, the network's distribution of self-couplings and g. We will first discuss the stable fixed points, which present qualitatively different structures depending on the values of the self-couplings. When both self-couplings s 1 , s 2 < 1, the only possibility for a stable fixed point is the trivial solution, with all x i = 0 (Fig. 2a), where the network activity quickly decays to zero. When at least one selfcoupling in greater than one, there are three stable fixed point regions ( Fig. 2a-i, a-ii, a-iii); in these three regions, the network activity starting from random initial conditions unfolds via a long-lived transient periods, then it eventually settles into a stable fixed point. This transient activity with late fixed points is a generalization of the network phase found in [15]. When both self-couplings are greater than one (s 1 , s 2 > 1) the fixed point distri-bution in each sub-population is bi-modal ( Fig. 2a-ii,iii). When s 1 > 1 and s 2 < 1, the solutions for the respective sub-populations are localized around bi-modal fixed points and around zero, respectively ( Fig. 2a-i). In the case of a Gaussian distribution of self-couplings in the stable fixed point regime, a complex landscape of stable fixed points emerges. The unit values at the stable fixed points continuously interpolates between the zero (for units with s i < 1) and the multiple values bi-modal cases (for units with s i > 1) within the same network (Fig. 3a).

IV. A RESERVOIR OF TIMESCALES
In the chaotic phase we can estimate the intrinsic timescale τ i of a unit x i from its autocorrelation function as the half width at its autocorrelation half maximum ( Fig. 2a-vi, τ 1 and τ 2 ). The chaotic phase in the network, Eq. (1), is characterized by a large range of timescales that can be simultaneously realized across the units with different self-couplings. In a network with two self-couplings s 1 and s 2 in the chaotic regime, we found that the ratio of the timescales τ 2 /τ 1 increases as we increase the self-couplings ratio s 2 /s 1 (Fig.  2b). The separation of timescales depends on the relative fraction n 2 /n 1 of the slow and fast populations. When this fraction approaches zero, (with n 1 → ∞), the log of the timescale ratio exhibits a supralinear dependence on the self-couplings ratio, as described analytically in Section V, leading to a vast separation of timescales. Other self-couplings ratios s 2 /s 1 approach the timescale supralinear separation as the self-couplings ratio increases, up to a saturation and decay back to one (not shown).
In a case of a lognormal distribution of self-couplings, in the chaotic regime the network generates a reservoir of multiple timescale τ i of chaotic activity across network units, spanning across several orders of magnitude ( Fig.  3b). For long tailed distributions such as the lognormal, mean field theory can generate predictions for rare units with large self-couplings from the tail end of the distribution by solving (2) and the continuous version of (3), see Appendix B, highlighting the exponential relation between a unit's self-coupling and its autocorrelation decay time.

V. SEPARATION OF TIMESCALES IN THE BISTABLE CHAOTIC REGIME
To gain an analytical understanding of the parametric separation of timescales in networks with heterogeneous self-couplings, we consider the special case of a network with two self-couplings where a large subpopulation (N 1 = N − 1) with s 1 = 1 comprises all but one slow probe unit, x 2 , with large self-coupling s 2 s 1 . The probe unit obeys the dynamical equa- Continuous distributions of self-couplings. a) In a network with a Gaussian distribution of self-couplings (mean µ = 1 and variance σ 2 = 9), and g = 2, the stable fixed point regime exhibits a distribution of fixed point values interpolating between around the zero fixed point (for units with si ≤ 1) and the multi-modal case (for units with si > 1). The purple curve represents solutions to x = s tanh(x). b) In a network with a lognormal distribution of self-couplings (parameters µ = 0.2 and σ 2 = 1), and g = 2, autocorrelation timescales τi in the chaotic phase span several orders of magnitude as functions of the units' self-couplings si (purple curve shows the dynamic mean-field predictions for τi).
the large N limit, we can neglect the backreaction of the probe unit on the mean field and approximate the latter as an external Gaussian colored noise η(t) with autocor- The noise η(t) can be parameterized by its strength, defined as D = ∞ 0 dτ C(τ ) and its timescale (color) τ 1 . For large s 2 , the dynamics of the probe unit x 2 can be captured by a bi-stable chaotic phase whereby its activity is localized around the critical points x 2 = x ± ±s 2 ( Fig. 4a-i) and switches between them at random times. In the regime of strong colored noise (as we have here, with τ 1 7.9 1), the stationary probability distribution p(x) (for x ≡ x 2 , Fig. 4a-ii,b) satisfies the unified colored noise approximation to the Fokker Planck equation [20,21]: given by The distribution p(x) has support in the region h(x) > 0 comprising two disjoint intervals |x| > x c where tanh(x c ) 2 = 1 − 1+τ1 τ1s2 (Fig. 4b). p(x) is concentrated around the two minima x ± ±s 2 of U ef f . The main effect of the strong color τ 1 1 is to sharply decrease the variance of the distribution around the minima x ± . This is evident from comparing the colored noise with a white noise, when the latter is driving the same bi-stable probe  8), its support excludes the shaded gray area) from the effective potential U ef f (dashed blue curve) captures well the activity histogram (same as (a-ii)); whereas the white noise distribution p(x2) (dark green curve, obtained from the naive potential U , dashed green curve) captures the probe unit's activity (same as (a-v)) when driven by white noise, and deviates significantly from the activity distribution when the probe is embedded in our network. (c) Average dwell times, T , in the bistable states. Simulation results, mean and 95% CI (blue curve and light blue background, respectively; An example for the full distribution of T is given in (a-iii)). Mean-field prediction (purple curve). The mean first passage time from the unified colored noise approximation (Eq. (10), black curve) captures well the simulated dwell times. An approximation for the unified colored noise (Eq. (11), gray dashed line) achieves good results as well. the white noise average dwell times are significantly different (green curve).
noise with an equivalent strength to the colored noise, Fig. 4a-iv,v,vi. The naive potential for the white noise case U = x 2 /2 − s 2 log cosh(x) is obtained from (8) by sending τ 1 → 0 in the prefactor h and in potential U ef f . It results in wider activity distribution compared to our network generated colored noise, in agreement with the simulations, Fig. 4a,b.
In our network generated colored noise the probe unit's temporal dynamics is captured by the mean first passage time T for the escape out of the potential well: We evaluated the integrals by steepest descent around The agreement between (10) and simulation results improves with increasing s 2 , as expected on theoretical ground [20,21] (Fig. 4c). The asymptotic scaling for large s 2 is In this slow probe regime, we thus achieved a parametric separation of timescales between the population x 1 , with its intrinsic timescale τ 1 , and the probe unit x 2 whose activity fluctuations exhibit two separate timescales: the slow timescale T of the bistable switching and the fast timescale τ 1 of the fluctuations around the metastable states (obtained by expanding the dynamical equation around the meta-stable values x ± = ±s 2 ). One can generalize this metastable regime to a network with N − p units which belong to a group with s 1 = 1 and p N slow probe units x α , for α = 2, . . . , p + 1, with large self-couplings s α . The slow dynamics of each probe unit x α is captured by its own bistable switching time T α in (10) and all slow units are driven by a shared external colored noise η(t) with timescale τ 1 . In summary, in our model multiple timescales can be robustly generated with specific values, varying over several orders of magnitude.
Is the relationship between the unit's self-coupling and its timescale relying on single-unit properties, or does it rely on network effects? To answer this question, we compare the dynamics of a unit when driven by a white noise input vs. the self-consistent input generated by the rest of the recurrent network (i.e., the mean field). If the neural mechanism underlying the timescale separation was a property of the single-cell itself, we would observe the same effect regardless of the details of the input noise. We found that the increase in the unit's timescale as a function of s 2 is absent when driving the unit with white noise, and it only emerges when the unit is driven by the self-consistent mean field. We thus concluded that this neural mechanism is not an intrinsic property of a single unit but requires the unit to be part of a recurrently connected network.

VI. A RESERVOIR OF TIMESCALES IN E-I SPIKING NETWORKS
We next investigated whether the neural mechanism endowing the random neural network (1) with a reservoir of timescales could be implemented in a biologically plausible model exhibiting spiking activity and excitatory/inhibitory cell-type specific connectivity. We modeled the local cortical circuit as a recurrent network of excitatory (E) and inhibitory (I) current-based leakyintegrated-and-fire neurons (see Appendix C for details), where both E and I populations were arranged in neural clusters (Fig. 5A) [22][23][24]. Synaptic couplings between neurons in the same cluster were potentiated compared to those between neurons in different clusters. Using mean field theory, we found that the recurrent interactions of cell-type specific neurons belonging to the same cluster can be interpreted as a self-coupling, which can be expressed in terms of the underlying network parameters as EE is the average synaptic couplings between E neurons within the cluster. The spiking network time-varying activity unfolds through sequences of metastable attractors [23,24], characterized by the activation of different subsets of neural clusters (Fig. 5A). This metastable regime is similar to the one observed in the random neural network for large self-couplings (Fig. 4) and the characteristic timescale T of a cluster's metastable dynamics can be estimated from its average activation time. We tested whether the heterogeneity in the cluster self-coupling dis-tribution could lead to a heterogeneous distribution of timescales. We endowed the clustered network with a heterogeneous distribution of cluster sizes; other sources of heterogeneity include the variability in average coupling strengthsJ (in) and the Erdos-Renyi connectivity, yielding altogether a heterogeneous distribution of selfcouplings (Fig. 5C). We found that the cluster activation timescales T varied across clusters spanning a large range from 20ms to 100s (Fig. 5C). The cluster timescale T was proportional to the value of a cluster's self-coupling. In particular, the functional dependence of log(T ) vs. selfcoupling s E was best fit by a quadratic polynomial (Fig.  5C, see Appendix C for details), in agreement with the analytical calculation in the random network model (11). We thus concluded that a reservoir of timescales can naturally emerge in biologically plausible models of cortical circuits from a heterogeneous distribution of cluster sizes. Both the range of timescales (20ms-100s) [10] and the distribution of cluster sizes (50-100 neurons) [25,26] are consistent with experimental observations.

VII. DEMIXING OF TIME-DEPENDENT BROADBAND INPUT
Previous work in random networks with no selfcouplings (s i = 0 in (1)) showed that stimulus-driven suppression of chaos is enhanced at a particular input frequency, related to the network's intrinsic timescale [16]. We investigated whether in our network with two different self-couplings s 1 < s 2 , in the chaotic regime, the stimulus-dependent suppression of chaos exhibited different features in the two sub-populations, depending on their different intrinsic timescale. We drove each network unit x i with an external broadband stimulus I i (t) = A L l=1 sin(2πf l t + θ i ) consisting of the superposition of L sinusoidal inputs of different frequencies f l in the range 1 − 200 Hz, with an equal amplitude A = 0.5 and random phases θ i . We found that the sub-population with a slow, or fast, intrinsic timescale preferentially entrained its activity with slower, or faster, spectral components of the broadband stimulus respectively (Fig. 6a). We quantified this effect using a spectral modulation index m(f ) = [(P 2 (f ) − P 1 (f ))/(P 2 (f ) + P 1 (f ))] (12) where P α (f ) is the power-spectrum peak of subpopulation α at the frequency f (Fig. 6b). A positive, or negative, value of m(f ) reveals a stronger, or weaker, respectively, entrainment at frequency f in the sub-population s 2 compared to s 1 . m(f ) exhibited a crossover behavior whereby the low frequency component of the input predominantly entrained the slow population s 2 , while the fast component of the input predominantly entrained the fast population s 1 . When fixing s 1 = 1 and varying s 2 , we found that the dependence of the crossover frequency f c on s 2 was strong at low input amplitudes and was progressively tamed at larger input amplitudes c , as a function of self coupling s2 for different input amplitudes. An inversely proportional relation between the cut off period and the amplitude of the broadband signal is present. (Fig. 6c). This is consistent with the fact that the input amplitude completely suppresses chaos beyond a certain critical value, as previously reported in network's with no self-couplings [16].

VIII. DISCUSSION
We demonstrated a new robust and biologically plausible network mechanism whereby multiple timescales emerge across units with heterogeneous self-couplings. In our model, units are interpreted as neural clusters, or functional assemblies, consistent with experimental evidence from cortical circuits [25][26][27][28][29] and theoretical modeling [23,24]. We found that the neural mechanism underlying the large range of timescales is the heterogeneity in the distribution of self-couplings (representing neural cluster strengths). We showed that this mechanism can be naturally implemented in a biologically plausible model of a neural circuit based on spiking neurons with excitatory/inhibitory cell-type specific connectivity. This spiking network represents a microscopic realization of our mechanism where neurons are arranged in clusters and a cluster's self-coupling represents the strength of the recurrent interactions between neurons belonging to that cluster. A heterogeneous distribution of cluster sizes, in turn, generates a reservoir of timescales.
Several experimental studies uncovered heterogeneity of timescales of neural activity across brain areas and species. Comparison of the population-averaged autocorrelations across cortical areas revealed a hierarchical structure, varying from 50ms to 350ms along the occipital-to-frontal axis [6]. Neurons within the same area exhibit a wide distribution of timescales as well. A heterogeneous distribution of timescales (from 0.5s to 50s) was found across neurons in the oculomotor system of the fish [11] and primate brainstem [30], suggesting that timescale heterogeneity is conserved across phylogeny. During periods of ongoing activity, the distribution of single-cell autocorrelation timescales in primates was found to be right-skewed and approximately lognormal, ranging from 10ms to 10s [10]. Single neuron activity was found to encode long reward memory traces in primate frontal areas over a wide range of timescales up to 10 consecutive trials [9]. In these studies, autocorrelation timescales where estimated using parametric fits, which may be affected by statistical biases, although a new Bayesian generative approach might overcome this issue [31]. In this study, we estimated timescales nonparametrically as the half-width at half-maximum of the autocorrelation function. In our biologically plausible model based on a spiking network with cell-type specific connectivity, the distribution of timescales was in the range between 20ms and 100s, similar to the range of timescales observed in experiments [10,11,30]. Moreover, the distribution of cluster sizes in our model is within the 50-100 neurons range, consistent with the size of functional assemblies experimentally observed in cortical circuits [25,26]. A fundamental new prediction of our model, to be tested in future experiments, is the direct relationship between cluster strength and its timescale.
Previous neural mechanisms for generating multiple timescales of neural activity relied on single cell biophysical properties, such as membrane or synaptic time constants [32]. In feedforward networks, developmental changes in single-cell conductance can modulate the timescale of information transmission, explaining the transition from slow waves to rapid fluctuations observed in the developing cortex [33]. However, the extent to which this single-cell mechanism might persist in pres-ence of strong recurrent dynamics was not assessed. To elucidate this issue, we examined whether a heterogeneous distribution of single-unit integration time constants could lead to a separation of timescales in a random neural network (see Appendix D for details). In this model, half of the units were endowed with a fast time constant which we held fixed, and the other half with a slow time constant, whose value we varied across networks. We found that, although the average network timescale increased proportionally to the value of the slower time constants, the difference in autocorrelation time between the two populations remained negligible. These results suggest that, although the heterogeneity in single-cell time constants may affect the dynamics of single neurons in isolation or within feedforward circuits [33], the presence of strong recurrent dynamics fundamentally alter these single-cell properties in a counterintuitive way. Our results suggest that a heterogeneity in single cell time constants may not lead to a diversity of timescales in presence of recurrent dynamics.
Our results further clarified that the relationship between a cluster's self-coupling and its timescale relies on the strong recurrent dynamics. This relationship is absent when driving an isolated cluster with white noise external input (Fig. 4). Indeed, the mechanism linking the self-coupling to the timescale only emerged when driving the unit with a mean field whose color was selfconsistently obtained from an underlying recurrent network of self-coupled units.
Previous models showed that a heterogeneity of timescales may emerge from circuit dynamics through a combination of structural heterogeneities and heterogeneous long-range connections arranged along a spatial feedforward gradient [34,35]. These networks can reproduce the population-averaged hierarchy of timescales observed across cortex in the range of 50-350ms [6,35]. A similar network architecture can also reproduce the heterogeneous relaxation time after a saccade, found in the brainstem oculomotor circuit [11,30], in a range between 10-50s [36,37]. This class of models can explain a timescale separation within a factor of 10, but it is not known whether they can be extended to several orders of magnitude, as observed between neurons in the same cortical area [10]. Moreover, while the feedforward spatial structure underlying these two models is a known feature of the cortical hierarchy and of the brainstem circuit, respectively, it is not known whether such a feedforward structure is present within a local cortical circuit. Our model, on the other hand, relies on strong recurrent connectivity and local functional assemblies, two hallmarks of the architecture of local cortical circuits [25][26][27][28][29]. Other network models generating multiple timescales of activity fluctuations were proposed based on self-tuned criticality with anti-hebbian plasticity [38], or multiple blockstructured connectivity [39].
In our model, the dynamics of units with large selfcouplings, exhibiting slow switching between bistable states, can be captured analytically using the universal colored noise approximation to the Fokker-Planck equation [20,21]. This is a classic tool from the theory of stochastic processes, which we successfully applied to investigate neural network dynamics for the first time. This slow switching regime may underlie the emergence of metastable activity, ubiquitously observed in the population spiking activity of behaving mammals [37,[40][41][42][43][44][45].
What is the functional relevance of neural circuits exhibiting a reservoir of multiple timescales? The presence of long timescales deeply in the chaotic regime is a new feature of our model which may be beneficial for memory capacity away from the edge of chaos [46]. Moreover, we found that, in our model, time-dependent broadband inputs suppress chaos in a population-specific way, whereby slow (fast) subpopulations preferentially entrain slow (fast) spectral components of the input. This mechanism may thus endow recurrent networks with a natural and robust tool to spatially demix complex temporal inputs [13] as observed in visual cortex [47]. Third, the presence of multiple timescales may be beneficial for performing flexible computations involving simultaneously fast and slow timescales, such as in role-switching tasks [12]; or as observed in time cells in the hippocampus [48,49]. A promising direction for future investigation is the exploration of the computational properties of our model in the context of reservoir computing [50] or recurrent networks trained to perform complex cognitive tasks [51].

ACKNOWLEDGMENTS
We would like to thank G. Mongillo and G. La Camera for discussions. LM was supported by National Institute of Neurological Disorders and Stroke grant R01-NS118461 and by National Institute on Drug Abuse grant R01-DA055439 (CRCNS). MS was supported by The Hebrew University of Jerusalem "Emergency response to covid19" grant.

A. DYNAMICAL REGIONS OF NETWORKS WITH IDENTICAL SELF-COUPLINGS, A SUMMARY
It is constructive to quickly survey the results of Stern et al. [15] who studied the special case of including a single value self-coupling s for all clusters in the network, P (s i ) = δ s,si . In this case the dynamics of all units in the network follow: Two variables determine the network dynamics, the network gain g and the self-coupling value s. The network gain g defines the strength of the network impact on its units. It brings the network into chaotic activity, without self-coupling (s = 0), for values g > 1 [14]. The self-coupling s generates bi-stability. Without network impact (g = 0) the dynamical equation 13 for each unit becomes which has two stable solutions for s > 1 (Appendix Fig.  1a), both at x = 0. For s < 1 (Appendix Fig. 1b) a single stable solution exists at x = 0.
When small values of network gain g are introduced to the network dynamics, Eq. (13), with identical bi-stable units (s > 1), each unit solution jitters around one of its two possible fixed points. After an irregular activity the network settle into a stable fixed point. This generates a region of transient irregular activity with stable fixed points (Appendix Fig. 1c). As g increases and s decreases, different possible fixed point configurations lose their stability (as a result, the typical time spent in the transient activity increases). When the critical line s c ≈ 1 + 0.157 ln (0.443g + 1) is crossed, no fixed point remains stable and the network activity becomes chaotic [15]. The "last" stable fixed point at the transition line has a unique configuration with all unit values located farthest from x = 0 (Appendix Fig. 1a, light green lines). Additional decrease of s and g leads to a region where any initial activity of the network decays and the trivial solution (x i = 0 for all i) is stable (Appendix Fig. 1c).

B. MEAN FIELD THEORY WITH MULTIPLE SELF-COUPLINGS
We derive the dynamic mean-field theory in the limit N → ∞ by using the moment generating functional [52,53]. For the derivation we follow the Martin-Siggia-Rose-De Dominicis-Janssen path integral approach formalism [54] as appears extensively in [18], we borrow their notations as well. For our model, Eq. (1), the moment generating functional is given by: In the region between these two curves, the network activity is chaotic. Colored circles denote, according to their locations on the phase diagram and with respect to their colors, the values of s (ranging from 1.6 and decreasing with steps of 0.2) and g = 1.5, used for the autocorrelation functions C(τ ) in (d) a . e) Widths at half peak (values of τ 's in the main text notation) of the autocorrelation functions in (d).
a Corrected version where Dx = i Dx i and Dx = i Dx i /2πi. To start, we calculate Z(J) J . We take advantage of the selfaveraging nature of our model, particularly by averaging over the quenched disorder, J. The couplings, J ij , are i.i.d. variables extracted from a normal distribution and appear only in the last term in (15). We hence focus our current calculation step on that term, and we derive the result to the leading term in N , yielding: The result above suggests that all the units in our network are coupled to one another equivalently (by being coupled only to sums that depend on all units' activity).
To further decouple the network, we define the quantity We enforce this definition by multiplying the disordered averaged moment generating functional with the appropriate Dirac delta function, δ, in its integral form: where dQ 2 is an integral over the imaginary axis (including its 1/(2πi) factor). We can now rewrite the disordered averaged moment generating functional, using (16) to replace its last term, the definition of Q 1 , and with multiplying the functional by the δ function above. All together we get: − sαφ(xα(t)) with n α = N α /N the fraction of units with self-couplings s α across the population, for α ∈ A. In the expression above we made use of the fact that Q 1 and Q 2 , now in a role of auxiliary fields, couple to sums of the fields x 2 i and φ 2 i and hence the generating functional for x i and x i can be factorized with identical multiplications of Z α . Note that in our network, due to the dependency on s i , x i -s are equivalent as long as s i -s are equivalent. Hence, the factorization is for Z α for all x i with s i = s α . Now each factor Z α includes the functional integrals Dx α and Dx α for a single unit with self-coupling s α .
In the large N limit we evaluate the auxiliary fields in (17) by the saddle point approach (we note variable valued at the saddle point by * ), obtaining: and yielding the saddle point values (Q * 1 , Q * 2 ): where C(τ ), with τ = f (t, t ), represents the average autocorrelation function of the network (as was defined in the main text). The second saddle point Q * 2 = 0 vanishes due to x α (t)x α (t ) = 0 as can be immediately extended from [18,52]. The action at the saddle point reduces to the sum of actions for individual, non-interacting units with self-coupling s α . All units are coupled to a common external field Q * 1 . Inserting the saddle point values back into Eq. (17), we obtain Thus in the large N limit the network dynamics are reduced to those of a number of A units x α (t), each represents the sub-population with self-couplings s α and follows dynamics governed by for all α ∈ A and where η(t) is a Gaussian mean field with autocorrelation The results above can be immediately extended for the continuous case of self-coupling distribution P (s) yielding: η(t)η(t ) = g 2 p(s)φ(x(s, t))φ(x(s, t ))ds (23) with p(s) the density function of the self-couplings distribution in the network and the units dynamics dependent on their respective self-couplings: C. SPIKING NETWORK MODEL Clustered network architecture. We simulated a recurrent network of N = 2000 excitatory (E) and inhibitory (I) spiking neurons with relative fractions n E = 80% and n I = 20% and connection probabilities p EE = 0.2 and p EI = p IE = p II = 0.5 (Fig. 5). Non-zero synaptic weights from pre-synaptic neuron j to post-synaptic neuron i were J ij = j ij / √ N , with j ij sampled from a normal distribution with mean j αβ , for α, β = E, I, and standard deviation δ 2 . Neurons were arranged in p cell-type specific clusters. E clusters had heterogeneous sizes drawn from a uniform distribution with mean of N clust E = 80 E-neurons and 30% standard deviation. The number of clusters was determined as p = round(n E N (1 − n bgr )/N clust E ), where n bgr = 0.1 is the fraction of background neurons in each population, i.e., not belonging to any cluster. I clusters were paired to E clusters and the size of each I cluster was matched to the corresponding E cluster with a proportionality factor n I /n E = 1/4. Neurons belonging to the same cluster had potentiated intra-cluster weights by a factor J + αβ , while those belonging to different clusters had depressed inter-cluster weights by a factor J − αβ , where: J + EI = p/(1 + (p − 1)/g EI ), J + IE = p/(1 + (p − 1)/g IE ), J − EI = J + EI /g EI , J − IE = J + IE /g IE and J − αα = 1−γ(J + αα −1) for α = E, I, with γ = f (2−f (p+1)) −1 . f = (1−n bgr )/p is the fraction of E neurons in each cluster. Parameter values are in Table I.
Single neuron dynamics. We simulated current-based leaky-integrate-and-fire (LIF) neurons, with membrane potential V and dynamical equation where τ m is the membrane time constant. Input currents included a contribution I rec from the other recurrently connected neurons and a constant external current I ext = N ext J α0 r ext (units of mV s −1 ), for α = E, I, representing afferent inputs from other brain areas and N ext = n E N p EE . When the membrane potential V hits the threshold V thr α (for α = E, I), a spike is emitted and V is held at the reset value V reset for a refractory period τ ref r . We chose the thresholds so that the homogeneous network (i.e.,where all J ± αβ = 1) was in a balanced state with average spiking activity at rates (r E , r I ) = (2, 5) spks/s. The post-synaptic currents evolved according to where τ s is the synaptic time constant, J ij are the recurrent couplings and t k is the time of the k-th spike from the j-th presynaptic neuron. Parameter values are in Table I. Self-couplings. We can estimate the E-cluster selfcouplings in this model using mean field methods [22,24].
Model parameters for clustered network simulations Parameter Description Value jEE mean E-to-E synaptic weights × √ N 0.6 mV jIE mean E-to-I synaptic weights × √ N 0.6 mV jEI mean I-to-E synaptic weights × √ N 1.9 mV jII mean I-to-I synaptic weights × √ N 3.8 mV jE0 mean E-to-E synaptic weights × √ N 2.6 mV jI0 mean I-to-I synaptic weights × √ N 2.3 mV δ standard deviation of the synaptic weight distribution 20% J + EE Potentiated intra-cluster E-to-E weight factor 14 J + II Potentiated intra-cluster I-to-I weight factor 5 gEI Potentiation parameter for intra-cluster I-to-E weights 10 gIE Potentiation parameter for intra-cluster E-to-I weights 8 rext Average baseline afferent rate to E and I neurons 5 spks/s V thr