A reservoir of timescales emerges in recurrent circuits with heterogeneous neural assemblies

  1. Merav Stern
  2. Nicolae Istrate
  3. Luca Mazzucato  Is a corresponding author
  1. Institute of Neuroscience, University of Oregon, United States
  2. Faculty of Medicine, The Hebrew University of Jerusalem, Israel
  3. Departments of Physics, University of Oregon, United States
  4. Mathematics and Biology, University of Oregon, United States
13 figures, 1 table and 1 additional file


Summary of the main results.

(a) Left: Microscopic model based on a recurrent network of spiking neurons with excitatory and inhibitory cell types, arranged in neural assemblies of heterogeneous sizes. Right: Phenomenological model based on a recurrent network of rate units. Each unit corresponds to an E/I neural assembly, whose size is represented by the unit’s self-couplings si. (b) A lognormal distribution of self-couplings (representing assemblies of different sizes) generates time-varying activity whose heterogeneous distribution of timescale fits population activity recorded from awake monkey orbitofrontal cortex (data from Cavanagh et al., 2016). (c) When driving our heterogeneous network with broadband time-varying input, comprising a superposition of sine waves of different frequencies, large and small assemblies preferentially entrain with low and high spectral components of the input, respectively, thus demixing frequencies into responses of different populations.

Dynamical and fixed point properties of networks with two self-couplings.

(a) Ratio of autocorrelation timescales τ2/τ1 of units with self-couplings s2 and s1, respectively (τi is estimated as the half width at half max of a unit’s autocorrelation function, see panels iii, iv), in a network with n1=n2=0.5 and g=2 and varying s1,s2. A central chaotic phase separates four different stable fixed point regions with or without transient activity. Black curves represent the transition from chaotic to stable fixed point regimes, which can be found by solving consistently Equation 15, Equation 16, and Equation 18 (using equal to 1 in the latter), see Methods ('Fixed points and transition to chaos' and 'Stability conditions') for details. (i, ii) Activity across time during the initial transient epoch (left) and distributions of unit values at their stable fixed points (right), for networks with N=1000 and (i) s1=3.2,s2=1.5, (ii) s1=3.2,s2=1.2. (iii, iv) Activity across time (left) and normalized autocorrelation functions C(τ)/C(0), (right) of units with (iii) s1=0.8,s2=1.5, (iv) s1=0.8,s2=3.2. (b) Timescales τ2,τ1 (left) and their ratio τ2/τ1 (right) for fixed s1=1 and varying s2, as a function of the relative size of the two populations n1=N1/N,n2=N2/N (at g=2, N=2000; average over 20 network realizations). All points in b. were verified to be within the chaotic region using Equation 18.

Continuous distributions of self-couplings.

(a–i) In a network with a Gaussian distribution of self-couplings (mean μ=1 and variance σ2=9), and g=2.5, the stable fixed point regime exhibits a distribution of fixed point values interpolating between around the zero fixed point (for units with si1) and the multi-modal case (for units with si>1). The purple curve represents solutions to x=stanh(x). (a,b–ii) A network with a lognormal distribution of self-couplings (parameters for (a,b) μ=0.2,0.5 and σ2=1,0.62, and g=2.5 ;autocorrelation timescales τi in units of ms) in the chaotic phase, span several orders of magnitude as functions of the units’ self-couplings si. (a-ii) Mean-field predictions for the autocorrelation functions and their timescales (purple curve) were generated from Equation 13 and Equation 14 via an iterative procedure, see Methods: 'Dynamic mean-field theory with multiple self-couplings' , 'An iterative solution'. (b) Populations of neurons recorded from orbitofrontal cortex of awake monkeys exhibit a lognormal distribution of intrinsic timescales (data from Cavanagh et al., 2016) (panel b-i, red), consistent with neural activity generated by a rate network with a lognormal distribution of self-couplings (panel b-i, blue; panel b-ii). ; We note that Cavanagh et al., 2016 use fitted exponential decay time constants of the autocorrelation functions as neurons’ timescales, while we use the half widths at half max of the autocorrelation functions as the timescales. To bridge these two definitions, we multiplied (Cavanagh et al., 2016) data by a factor of ln(2) before comparing it with our model and presenting it in this figure. The model membrane time constant was assumed to be 3 ms in this example.

Figure 4 with 1 supplement
Separation of timescales and metastable regime.

(a) Examples of bistable activity. (i, iv,i) - time courses; (ii, v) - histograms of unit’s value across time; (iii, vi) - histograms of dwell times. (a–i, ii, iii) An example of a probe unit x2 with s2=5, embedded in a neural network with N=1000 units, N1=N1 units with s1=1 and g=1.5. (a–iv, v, vi) An example of a probe unit driven by white noise. Note the differences in the x-axis scalings of the timecourses (a–i vs. a–iv and a–iii vs. a–vi).(b) The unified colored noise approximation stationary probability distribution (dark blue curve, Equation 19, its support excludes the shaded gray area) from the effective potential Ueff (dashed blue curve) captures well the activity histogram (light blue area; same as (a–ii)); whereas the white noise distribution (dark green curve, obtained from the naive potential U, dashed green curve) captures the probe unit’s activity (light green area; 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 of the full distribution of the dwell times is given in (a-iii)). Mean-field predictions (purple curve) were generated by calculating the average dwell times from a trace of x2, which was produced by solving the mean-field equations; Equation 2 simultaneously and consistently with Equation 3 with n1=1 and n2=0. The mean first passage time from the unified colored noise approximation (Equation 22, black curve), and for a simplified estimate thereof (Equation 4, gray dashed line) capture well T. When driven by white noise (green curve and light green curve are simulation results and simplified estimate using Equation 4, respectively), the probe’s average dwell times are orders of magnitude shorter than with colored noise, exhibiting substantial support of the probe distribution in the region where the crossing between wells happens (allowing frequent crossing,(a-iv) green line at x=0) and, equivalently, the low value of the potential around its maxima ((b) green dashed line at x=0). Comparison of white and colored noise demonstrates the central role of the self-consistent colored noise to achieve long dwell times.

Figure 4—figure supplement 1
Validation of universal colored noise approximation (UCNA) approach to estimate escape times.

’Escape points’ from one metastable state to the next were estimated from simulations as follows: for a transition from the x=s2 well towards the x+=x2 well (Left panel; analogous calculation for transitions x+x, Center panel), the escape point xesc was defined as the point where the trajectory starts accelerating towards the target well (positive second derivative). The distribution of escape points (Right panel) predominantly lied outside of the forbidden region (93.8% of escape points had |xesc|>xc). Parameters: same as in Figure 4a–i, s2=5. See Methods ('Universal colored-noise approximation to the Fokker-Planck theory') for further details.

Heterogeneity of timescales in E-I spiking networks.

(a) Top: Schematic of a spiking network with excitatory (black) and inhibitory populations (red) arranged in assemblies with heterogeneous distribution of sizes. Bottom: In a representative trial, neural assemblies activate and deactivate at random times generating metastable activity (one representative E neuron per assembly; larger assemblies on top; representative network of N=10,000 neurons), where larger assemblies tend to activate for longer intervals. (b) The average activation times <T> of individual assemblies (blue dots; the average was calculated across 100s simulation and across all neurons within the same assembly for all assemblies in 20 different network realizations; self-coupling units are in [mV], see Methods section). Fit of log(T)=a2sE2+a1sE+a0 with a2=0.14,a1=1.97,a0=5.51 (pink curve). Inset: cross-validated model selection for polynomial fit. As the assembly strength (i.e. the product of its size and average recurrent coupling) increases, <T> increases, leading to a large distribution of timescales ranging from 20 ms to 100s. (c) Eigenvalue distribution of the full weight matrix J (brown) and the mean-field-reduced weight matrix JMF (pink). (d) The Schur eigenvectors of the weight matrix JMF show that the slow (gapped) Schur eigenvalues (top) are associated with eigenvectors corresponding to E/I cluster pairs (bottom). See Appendix (e) Spiking network model for more details and for the scaling to larger networks.

Network response to broadband input.

(a) Power spectrum density of a network driven by time-dependent input comprising a superposition of 11 sinusoidal frequencies (see main text for details). Maroon and navy curves represent average power spectrum density in s1 and s2 populations, respectively; circles indicate the peak in the power spectrum density amplitudes at each frequency; amplitude A = 0.5, g=3.0, s1=1, and s2=4. (b) Modulation index of the power spectrum density amplitudes as a function of frequency in networks with s1=1 and various s2. The blue circles mark the cutoff frequency fc where the modulation index changes sign. (c) Cutoff period, 2πωc1, as a function of self-coupling s2 for different input amplitudes. An inversely proportional relation between the cutoff period and the amplitude of the broadband signal is present.

Appendix 1—figure 1
Transition to chaos with multiple self-couplings: Fixed point solutions and stability.

(a–i) The fixed point curve xαsαtanhxα, from Equation 15, for sα>1. Stable solutions are allowed within the dark green region. (b–i) The shape of a unit’s contribution to stability q1=(sαcoshxα2)2, from Equation 18. Stable solutions of xαsαtanhxα=η, 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 q1) is stable, light green region in (a–i).(a–ii) The curve xαsαtanhxα 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).

Appendix 1—figure 2
Network dynamics with identical self-couplings, adopted from Stern et al., 2014.

(a,b) Graphical solutions to Equation 29. (a) For s>1 there are two stable non-zero solutions (full black circles) and an unstable solution at zero (open black circle). The green background over the x axis denotes the regions of allowed activity values at the stable fixed point on the transition line to chaos (solid red curve in (c)).(b) For s<1 there is a single stable solution (full black circle) at zero. (c) Regions of the network dynamics over a range of s and g values. Below the long dashed blue line, any initial activity in the network decays to zero. Above the solid red curve, the network exhibits transient irregular activity that eventually settles into one out of a number of possible nonzero stable fixed points. 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; corrected version of Figure 4a in Stern et al., 2014). (e) Widths at half peak (values of τ’s in the main text notation) of the autocorrelation functions in (d).

Appendix 2—figure 1
Synaptic weight matrices of clustered spiking networks.

Synaptic weight matrices of a clustered spiking networks. (a) Left: Synaptic weight matrix J of a clustered spiking network with N=2000 neurons, exhibiting the E/I four submatrices, where diagonal blocks within each submatrix reveal the E/I clustered structure (larger to smaller clusters, top to bottom in each submatrix). Right: 2(p+1)-dimensional mean-field-reduced synaptic weight matrix JMF corresponding to the full matrix on the left. Populations are ordered as follows: p+1 excitatory clusters (larger to smaller, top to bottom), background excitatory population, p+1 inhibitory clusters, background inhibitory population. For each pair of populations, the value in JMF is obtained from the corresponding block in J by summing its column values (presynaptic inputs) and averaging over the rows (postsynaptic neurons belonging to the population). (b) Right: Schur decomposition of JMF=VTVT into an upper triangular matrix T. Left: Schur eigenvectors (columns of V, bottom) sorted from larger to smaller Schur eigenvalue (top). Larger Schur eigenvalues are associated with eigenvectors whose loadings are on larger E/I cluster pairs.

Appendix 2—figure 2
Metastable attractor dynamics in clustered spiking networks.

Metastable attractor dynamics in clustered spiking networks. (a) Metastable activity from all clustered neurons in a representative trial of a network with N=2000 neurons (larger to smaller clusters, top to bottom; action potentials from excitatory and inhibitory neurons, black and red, respectively). (b) Firing rate distributions of four representative neurons from the network in (a) exhibit bimodal distributions (colors represent clusters with different self-couplings; spike counts estimated in 100ms bins over 20 trials of 200 seconds duration). (c) Left: Metastable attractor dynamics unfolds through different network configurations with a range of simultaneously active clusters from one to four (average occurrence of each configuration across 20 networks as fraction of total simulation time). Inset: Metastable dynamics are uncorrelated across clusters (distribution of pairwise Pearson correlations between clusters’ firing rates, 0.01±0.12, mean±SD; 20 networks). Right: The firing rate of neurons in active clusters depend on how many clusters are simultaneously active, with higher firing rates in configurations with less co-active clusters (mean and SD across 20 networks).

Appendix 2—figure 3
Metastable activity and timescale distribution in networks of N=2000,5000 neurons (top and bottom rows, respectively).

Panels (a) and (b) have the same notations as Figure 5b, and Figure 5d, respectively. The fits of average activation time T vs. self-coupling log(T)=a2sE2+a1sE+a0 yielded a2=0.44,a1=2.06,a0=1.04 for N=2000; and a2=0.093,a1=0.88,a0=0.45 for N=5000.

Appendix 2—figure 4
Chaotic rate network of N=200 units (self-couplings drawn from a lognormal distribution with parameters μ=0.2,σ=2).

(a) Eigenvalue distribution (brown) and Schur eigenvalues (pink) of the synaptic weight matrix J. (b) The Schur eigenvectors of J corresponding to large positive Schur eigenvalues have loadings localized on units with larger self-couplings, with a similar structure to the Schur eigenvectors in the spiking networks in Appendix 2-figure 3 .

Appendix 3—figure 1
Timescale analysis for an RNN with two-time constants θi, Equation 31, governing equal populations of neurons (N1=N2=1000) and gain g=2.5.

(a) Average autocorrelation function for each population. The insert shows the dynamics of individual neurons from each population: blue for neurons with timeconstant θ1 and green for neurons with timeconstant θ2. In the networks considered here, θ1=0.1 ms is kept constant while: θ2=0.1 ms (i), θ2=1.0 ms (ii), θ2=10.0 ms (iii), θ2=100.0 ms (iv), θ2=1000.0 ms (v). (b) Population timescale ratio τ2/τ1 for fixed timeconstant θ1=0.1 ms and varying θ2.


Table 1
Parameters for the clustered network used in the simulations.
Model parameters for clustered network simulations
JEEE-to-E synaptic weights0.9/N [mV]
JIEE-to-I synaptic weights0.9/N [mV]
JEII-to-E synaptic weights2.7/N [mV]
JIII-to-I synaptic weights5.4/N [mV]
JE0E-to-E synaptic weights3.7/N [mV]
JI0I-to-I synaptic weights3.3/N [mV]
JEE+Potentiated intra-assembly E-to-E weight factor14N/2000
JII+Potentiated intra-assembly I-to-I weight factor5N/2000
gEIPotentiation parameter for intra-assembly I-to-E weights
gIEPotentiation parameter for intra-assembly E-to-I weights8N/2000
rextAverage baseline afferent rate to E and I neurons5 [spks/s]
VEthrE-neuron threshold potential1.43 [mV]
VIthrI-neuron threshold potential0.74 [mV]
VresetE- and I-neuron reset potential0 [mV]
τmE- and I-neuron membrane time constant20 [ms]
τrefrE- and I-neuron absolute refractory period5 [ms]
τsE- and I-neuron synaptic time constant5 [ms]

Additional files

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Merav Stern
  2. Nicolae Istrate
  3. Luca Mazzucato
A reservoir of timescales emerges in recurrent circuits with heterogeneous neural assemblies
eLife 12:e86552.