<figure><img src="image_1.png" /></figure>

HAL
open science

Searching for long timescales without fine tuning

Xiaowen Chen, William Bialek

► To cite this version:

Xiaowen Chen, William Bialek. Searching for long timescales without fine tuning. Physical Review E, 2024, 110 (3), pp.034407. 10.1103/PhysRevE.110.034407. hal-04885441

HAL Id: hal-04885441

https://hal.science/hal-04885441v1

Submitted on 14 Jan 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Searching for long time scales without fine tuning

Xiaowen Chen¹*† and William Bialek¹,²‡

¹Joseph Henry Laboratories of Physics, and Lewis–Sigler Institute
for Integrative Genomics, Princeton University, Princeton NJ 08544

²Initiative for the Theoretical Sciences, The Graduate Center,
City University of New York, 365 Fifth Ave, New York NY 10016

(Dated: August 22, 2024)

Animal behavior occurs on time scales much longer than the response times of individual neurons. In many cases, it is plausible that these long time scales emerge from the recurrent dynamics of electrical activity in networks of neurons. In linear models, time scales are set by the eigenvalues of a dynamical matrix whose elements measure the strengths of synaptic connections between neurons. It is not clear to what extent these matrix elements need to be tuned in order to generate long time scales; in some cases, one needs not just a single long time scale but a whole range. Starting from the simplest case of random symmetric connections, we combine maximum entropy and random matrix theory methods to construct ensembles of networks, exploring the constraints required for long time scales to become generic. We argue that a single long time scale can emerge generically from realistic constraints, but a full spectrum of slow modes requires more tuning. Langevin dynamics that generates patterns of synaptic connections drawn from these ensembles involves a combination of Hebbian learning and activity-dependent synaptic scaling.

# I. INTRODUCTION

Living systems face various challenges over their lifetimes, and responding to these challenges often involves behaviors that occur over multiple time scales [1]. As an example, a migratory bird needs both to react to instantaneous gusts and to navigate its course over months [2]. Thanks to recent development of behavior tracking and image analysis techniques, long time scales have been identified and quantified in the behaviors of many animals, including the scale invariant correlations spanning three decades in time in the fruit fly [1], the near-marginal modes in locally linear approximations to neural and behavioral dynamics of the nematode *C. elegans* [3], and the non-Markovian behavior in the larval zebrafish [4].

How do those time scales emerge from the generator of the behavior, the nervous system? While transient responses of individual neurons decay on the time scale of tens of milliseconds, autonomous behaviors can last orders of magnitude longer. We see this when we hold a string of numbers in our heads for tens of seconds before dialing a phone, and when a musician plays a piece from memory that lasts many minutes. While this emergence of long time scales can arise from the coupling of living systems to a slowly fluctuating environment [5–7] such as the day-night cycle, experimental evidences support that long time scales can be generated by interactions of the neurons. On one hand, single long time

scales in behavior have been associated with persistent neural activities, where after a pulse stimulation, some neurons are found to hold their firing rate at specific values that encode the transient stimuli [8–13]. On the other hand, a wide range of persistent time scales has been observed, e.g. in various cortical areas of monkeys [14, 15] and macaque prefrontal cortex [16], as well as within the same local circuit such as the larval zebrafish oculomotor system [17]. Moreover, in the mouse hippocampus, coarse-graining neuronal activities has uncovered reproducible scaling behaviors of networks of 1000+ neurons [18]. As verified both biologically and computationally, these multidimensional dynamics allow the network to integrate information at a broad range of timescales, which are beneficial for decision making [19], memory storage [20], task switching [16], and better prediction of multiscale dynamics [21]. All of these observations motivate a search for mechanisms allowing the emergence of long time scales.

It is plausible that long time scales emerge from the recurrent dynamics of electrical activity in the network of neurons. In the simplest linear model, the relaxation times of the system depend on the eigenvalues of a matrix representing the synaptic connection strengths among neurons, and we can imagine this being tuned so that time scales become arbitrarily long [22]. This simple model has successfully explained the long time scale in the oculomotor system of goldfish, where the nervous system tunes its dynamics to be slow and stable using constant feedback from the environments [10, 11]. In general, long time scales in linear dynamical systems require fine tuning, as the modes need to be slow, but not unstable.

To circumvent this biologically-implausible fine tuning requirement, current approaches of modeling interacting networks which generate long time scales include either changing the setup completely by adding in non-

<footnote>*xiaowenc@princeton.edu</footnote>

<footnote>†Present address: Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, PSL Université, CNRS, Sorbonne Université, Université Paris Cité, F-75005 Paris, France.</footnote>

<footnote>xiaowen.chen@phys.ens.fr</footnote>

<footnote>‡wbialek@princeton.edu</footnote>
linearity [23], fine tuning the structure of the networks using biological structures [19] or highly structured networks such as feedforward [24], ring [25], or grid networks [26], imposing clustered connections [27], adding hetereogeneity to the self-interaction [28], or promoting the interaction matrix to dynamical variables regulating the overall neural activity [29–31]. Nonetheless, it remains unclear how much fine-tuning of the interactions is required to generate long time scales. If we allow all possible connections, what constraints on the strengths of these connections are required for the emergence of long time scales? We can ask this question both for the case of a single long time scale and for a continuum spectrum time scales, as observed in animal behaviors [32].

In this manuscript, we address the fine-tuning question by asking whether we can find ensembles of random connection matrices, subject to biologically plausible constraints, such that the resulting time scales of the system grow with increasing system sizes. We also discuss the conditions for systems to exhibit a continuous spectrum of slow modes as opposed to single slow modes. Finally, we present plausible dynamics for the system to tune its connection matrix towards these desired ensembles, with the hope that these dynamics may correspond to actual regulating process in real biological systems.

## II. SETUP

The problem of characterizing time scales in fully nonlinear neural networks (or any high dimensional dynamical system) is very challenging. To make progress, we follow the example of Ref. [22] and consider the case of linear networks, as to approximate the dynamics near fixed points. For linear systems, time scales are related to the eigenvalues of the dynamical matrix that embodies the pattern of synaptic connectivity. The question of whether behavior is generic can be made precise by drawing these matrices at random from some probability distribution, connecting with the large literature on random matrix theory [33, 34]. Importantly, we expect that some behaviors in these ensembles of networks become sharp as the networks become large, a result which has been exploited in thinking about problems ranging from energy levels of quantum systems [35] to ecology [36] and finance [37].

Concretely, we represent the activity each neuron $i = 1, 2, \dots, N$ by a continuous variable $x_i$, which we might think of as a smoothed version of the sequence of action potentials, and assume a linear dynamics

$$ \dot{x}_i = -x_i + \sum_j M_{ij}x_j + \eta_i(t). \quad (1) $$

See Fig. 1(a) for schematics. For simplicity, we assume the interaction matrix (i.e. the synaptic strength) $M$ is symmetric, i.e. $M_{ij} = M_{ji}$ (see Sec. II A and Appendix E for discussions on general $M$ without symmetry constraint). If the neurons were unconnected ($M = 0$), their

activity $x$ would relax exponentially on a time scale which we choose as our unit of time. In what follows it will be important to imagine that the system is driven, at least weakly, and we take these driving terms to be independent in each cell and uncorrelated in time, where $\langle \eta_i(t) \rangle = 0$ and

$$ \langle \eta_i(t) \eta_j(t') \rangle = 2\delta_{ij}\delta(t-t'). \quad (2) $$

The choice of white noise is conventional, but also important because we want to understand how time scales emerge from the network dynamics rather than being imposed upon the network by outside inputs.

In linear systems we can rotate to independent modes, corresponding to weighted combinations of the original variables. If the matrix $M$ is symmetric, then the dynamics are described by the relaxation times of these modes,

$$ \tau_i \equiv \frac{1}{1 - \lambda_i} = \frac{1}{k_i}, \quad (3) $$

where $\{\lambda_i\}$ are the eigenvalues of $M$; the system is stable only if all $\lambda_i < 1$.

If the matrix $M$ is chosen from a distribution $P(M)$ then the eigenvalues are random variables, but their density, for example, becomes smooth and well defined in the limit $N \to \infty$,

$$ \rho(\lambda) \equiv \lim_{N \to \infty} \frac{1}{N} \sum_{i=1}^{N} \delta(\lambda - \lambda_i). \quad (4) $$

The width of the support, i.e. where the density of eigenvalues is nonzero, is denoted by $l$. This density of eigenvalues, $\rho(\lambda)$, is the spectrum of the interaction matrix $M$, and determines the timescales in the system (see Fig. 1(b).) Systems with long time scales are defined such that the time scales grow as the system size $N$ or faster. Specifically, we focus on two types of time scales:

1. Longest time scale. The longest time scale, $\tau_{\max}$, is the time constant given by the slowest mode of the system, which dominates the dynamics after long enough times after a perturbation (Fig. 1(c)). This time scale is determined by the gap, $g_0$, between the largest eigenvalue and the stability threshold, which with our choice of units is $w = 1$. Mathematically, we define

$$ \tau_{\max} \equiv \frac{1}{g_0} = \frac{1}{1 - \lambda_{\max}}. \quad (5) $$

In the thermodynamic limit, the gap is taken to be between the stability threshold and the right edge of the support of the spectral density<fnref n="1" /> (see Fig. 1(b)).

<footnote n="1">¹ We note that this approximation is not ideal, as the spectral distribution of eigenvalues does not converge uniformly. In some cases, the fluctuation of the largest eigenvalue can be more meaningful than the average.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1: Schematics for symmetric linear dynamical systems, and the corresponding timescales. (a) A linear dynamical system with damping and pairwise interaction $M$ has time scales determined by the eigenvalue spectrum of $M$, especially the gap $g_0$ to the stability threshold (panel (b)) We focus on two timescales. (c) The longest time scale is defined as $\tau_{max} = 1/g_0$, and gives the time scale for the norm activity to decay when the system is perturbed at the time indicated by the red arrow. (d) The correlation time $\tau_{corr}$ is the characteristic time scale of de-correlation when the system is unperturbed. In cases where the system has a continuous range of long time scales, the correlations decay as a power law (red curve).</figcaption></figure>

2. Correlation time. The correlation time scale is the time scale it takes for the system to de-correlate, i.e. the decaying time constant of the correlation function (Fig. 1(d)). To get at the correlation time, let's take seriously the idea that the network is driven by noise. Then $x(t)$ become a stochastic process, and from Eqs. (1) and (2) we can calculate the correlation function

$$ 
\begin{aligned}
C_N(t) &\equiv \frac{1}{N} \sum_i \langle x_i(0)x_i(t) \rangle \\
&= \frac{1}{N} \sum_i \frac{1}{1-\lambda_i} e^{-(1-\lambda_i)|t|} \\
&= \frac{1}{N} \sum_i \tau_i e^{-|t|/\tau_i}.
\end{aligned} 
\quad (6) $$

The normalized correlation function

$$ R_N(t) \equiv \frac{C_N(t)}{C_N(0)} = \frac{\sum_i \tau_i e^{-|t|/\tau_i}}{\sum_i \tau_i} \quad (7) $$

has the intuitive behavior of starting at $R_N(0) = 1$ and decaying monotonically. Then there is a natural definition of the correlation time, by analogy with single exponential decays,

$$ \tau_{corr}(\{\lambda\}) \equiv \int_{0}^{\infty} dt R_N(t) = \frac{\sum_{i} \tau_i^2}{\sum_{i} \tau_i}. \quad (8) $$

In the thermodynamic limit, the autocorrelation coefficient $R(t)$ and the correlation time $\tau_{corr}$ becomes the ratio of two integrals over the eigenvalue density $\rho(\lambda)$.

Importantly, $\tau_{max}$ depends only on the largest eigenvalue, while $\tau_{corr}$ depends on the entire spectrum, and hence can be used to differentiate cases where the system is dominated by a single vs. a continuous spectrum of slow modes. The two time scales satisfy $\tau_{corr} \le \tau_{max}$, with the equality assumed only when all eigenvalues are equal, i.e. the spectral density is a delta function at $\lambda = \lambda_{max}$.

With these definitions, we refine our goal as to find biologically plausible ensembles for the connection matrix $M$, such that the resulting stochastic linear dynamics has time scales, $\tau_{max}$ and $\tau_{corr}$, that are “long,” growing as a power of the system size $N$, perhaps even extensively. To avoid fine tuning, we will construct examples of such ensembles by imposing global constraints on measurable observables of the dynamical system. We then compute the spectral density and the corresponding time scales using a combination of mean-field theory and numerical sampling of finite systems.

### A. Symmetric vs. general choice of interaction matrix $M$

The most general choice for the interaction matrix $M$ is without the symmetry constraint. Nonetheless, there are a few problems with this generality. First, the eigenvalues can now be complex, and the dynamics is composed of different modes of damped oscillations, introducing extra measures of time scales. Second, the correlation function for general interaction matrix depends not only on the eigenvalues, but also on the eigenvectors, vastly complicating the calculation (see Appendix E). In fact, in later parts of this manuscript, when we introduce stability constraints to the interaction matrix, random matrices without symmetry constraints becomes analytically unsolvable, while numerical results using brute force Monte Carlo sampling are limited to small system size. As suggested by the numerical results, asymmetric random matrices will exhibit similar fundamentals as symmetric random matrices, such as where possible phase transitions occur, and whether one can obtain slow modes. However, the scaling coefficient is different (Appendix E). For starting somewhere, we propose to leave solving the dynamics for general interaction matrices for the future, and will focus our discussion in this manuscript on symmetric matrices.

### III. TIME SCALES FOR ENSEMBLES WITH DIFFERENT GLOBAL CONSTRAINTS

To construct linear dynamical systems that generate long time scales without fine tuning individual parameters, we want to find probability distributions $P(M)$
<figure><img src="image_4.png" /><figcaption>FIG. 2: Spectral density for the connection matrix $M$ drawn from (a) Model 1, the Gaussian Orthogonal Ensemble (GOE); (b) Model 2, the GOE with a hard threshold enforcing stability; and (c) Model 3, the GOE with an additional global constraint on the norm activity. Three representative parameters are chosen for each ensemble such that the system is subcritical ($c = 0.6$, red), critical ($c = 1/\sqrt{2}$, blue), and supercritical ($c = 0.8$, black). The stability threshold is visualized as the dashed gray line at $\lambda_w = 1$. (d) Schematic for constraining the averaged norm activity in Model 3, generating ensembles with the soft threshold.</figcaption></figure>

for the connection matrix such that time scales are long, on average. We start with the Gaussian distribution, and gradually add constraints. We will see that for the Gaussian distribution itself, there is a critical value for the scaled variance of the elements of $M$, $c_{crit} = 1/\sqrt{2}$. For $c < c_{crit}$ the system is stable but time scales are short, while for $c > c_{crit}$ the system is unstable. Exactly at the critical point $c = c_{crit}$ time scales are long in the sense that we have defined, diverging with system size. The essential challenge is to find the weakest constraints on $P(M)$ that make these long time scales generic. Many of the results that we need along the way are known in the random matrix theory literature, but we will arrive at some new theoretical questions.

### A. Model 1: the Gaussian Orthogonal Ensemble

The simplest ensemble for the interaction matrix $M$ is the Gaussian Orthogonal Ensemble (GOE), which has been studied since the beginning of random matrix theory [33, 35, 38, 39]. GOE has symmetric matrix elements, $M_{ij} = M_{ji}$, and the matrix elements are independent

Gaussian random variables, with variances such that

$$M_{ii} \sim \mathcal{N}(0, c^2/N), \quad (9)$$

$$M_{ij}|_{i \neq j} \sim \mathcal{N}(0, c^2/2N). \quad (10)$$

As the parameter $c$ sets the scale of connections, we will refer to this parameter as the interaction strength. The factor of $N$ in the variance ensures that the density $\rho(\lambda)$ has support over a range of eigenvalues that are $O(1)$ at large $N$. The eigenvalue distribution of the GOE was solved a long time ago [35]. Here, we review these results, and use the eigenvalue distributions to characterize the longest and the correlation time scales. Despite its simplicity, Model 1 (GOE) provides a foundation for understanding the challenge of obtaining long time scales, before we introduce more complex random matrix ensembles.

The probability distribution for the interaction matrix $M$ can be written as

$$P_{\text{GOE}}(M) \propto \exp \left( -\frac{N}{2c^2} \text{Tr} M^{\top} M \right). \quad (11)$$

Because the distribution only depends on matrix traces, it is invariant to rotations of $M$. Equivalently, if we think of decomposing the matrix $M$ into its eigenvectors and eigenvalues, the probability depends only on the eigenvalues. Thus, we can integrate out the eigenvectors, and obtain the joint distribution of eigenvalues (see Appendix A for details),

$$P_{\text{GOE}}(\{\lambda_i\}) \propto \exp \left[ -\frac{N}{2c^2} \sum_i \lambda_i^2 + \frac{1}{2} \sum_{j \neq k} \ln |\lambda_j - \lambda_k| \right], \quad (12)$$

where the logarithmic repulsion term emerges from the Jacobian when we change variables from matrix elements to eigenvalues and eigenvectors. For a bit more intuition, we can identify that this joint density function of eigenvalues is mathematically identical to the Boltzmann distribution of two-dimensional Coulomb gas confined to one dimension (the line of real numbers), where each of the $N$ eigenvalues experiences an external quadratic potential, and pairs of electrons repel logarithmically (hence the two-dimensional Coulomb repulsion).

The spectral density can then be found using a mean field approximation, which becomes exact as $N \to \infty$. The result is Wigner's well-known semicircle distribution [35],

$$\rho_{\text{GOE}}(\lambda) = \frac{1}{\pi c} \sqrt{2 - \frac{\lambda^2}{c^2}}, \quad \lambda \in [-\sqrt{2c}, \sqrt{2c}]. \quad (13)$$

The width of the support is $l = 2\sqrt{2c}$. For completeness, we review the derivation of the spectral density in Appendix A.

Equation (13) for the spectral density, together with Fig. 2(a), shows that there is a phase transition at $c_{crit} = 1/\sqrt{2}$. If the interaction strength is greater than
this critical strength (supercritical), then $\lambda_{max} > 1$ and the system becomes unstable. If the interaction strength is smaller (subcritical), then the gap size between the largest eigenvalue and the stability threshold $\lambda = 1$ is of order 1, and the time scales remain finite as system size increases. The only case when the system has slow modes is at the critical value of the interaction strength, $c = c_{crit} = 1/\sqrt{2}$, where the spectral density becomes tangential to the stability threshold.

At criticality, corresponding to the blue curve in Fig. 2(a), we can estimate the size of the gap by asking that the gap be large enough to contain ~ 1 mode, that is

$$\int_{1-g_0}^{1} d\lambda N \rho_{\text{GOE}}(\lambda) \sim 1. \qquad (14)$$

Near $\lambda \to 1$, the spectral distribution can be approximated by $\rho_{GOE}(\lambda) \sim (1-\lambda)^{1/2}$, and this gives $Ng_0^{3/2} \sim 1$ or $g_0 \sim N^{-2/3}$. Thus the longest time scale grows with system size, $\tau_{max} \sim N^{2/3}$.

In the same way, we can estimate the full correlation function from the correlation time [Eq. (6)],

$$C_N(t) = \frac{1}{N} \sum_i \tau_i e^{-t/\tau_i}.$$

In the thermodynamic limit, $N \to \infty$,

$$C_{N \to \infty}(t) = \int d\lambda \frac{\rho_{\text{GOE}}(\lambda)}{(1-\lambda)} e^{-(1-\lambda)|t|} \\ \sim \int d\lambda \frac{1}{(1-\lambda)^{1/2}} e^{-(1-\lambda)|t|} \qquad (15) \\ \sim |t|^{-1/2}.$$

This has the power-law behavior expected for a critical system, where there is a continuum of slow modes.

Note that in the GOE system, slow modes with time scales growing as system size are only possible at a single value of the interaction strength. Nonetheless, we need to distinguish the fine tuning here as happening at an ensemble level, which is different from the element-wise fine tuning that might have been required if we considered particular interaction matrices.

## B. Model 2: GOE with hard stability threshold

Drawing interaction matrices from the GOE leads to long time scales only at a critical value of interaction strength. Can we modify the ensemble such that long time scales can be achieved without this fine tuning? In particular, in the GOE, if the interaction strength is too large, then the system becomes unstable. What will the spectral distribution look like if we allow $c > c_{crit}$ but edit out of the ensemble any matrix that leads to instability?

Mathematically, a global constraint on the system stability requires all eigenvalues to be less than the stability

threshold, $\lambda_w = 1$. This modifies the distribution of the interaction matrices with a Heaviside step function:

$$P_{\text{hard}}(\mathbf{M}) \propto P_{\text{GOE}}(\mathbf{M}) \prod_{i} \Theta(1 - \lambda_{i}). \qquad (16)$$

Conceptually, what this model does is to pull matrices out of the GOE and discard them if they produce unstable dynamics; the distribution $P_{hard}(\{\lambda\})$ describes the matrices that remain after this editing. Importantly we do not introduce any extra structure, and in this sense $P_{hard}$ is a maximum entropy distribution, as discussed more fully below.

The spectral density $\rho_{hard}(\lambda)$ that follows from $P_{hard}$ was first found by Dean and Majumdar [40–42]. Again, there is a phase transition depending on the interaction strength. For ensembles with interaction strength less than the critical value $c_{crit} = 1/\sqrt{2}$, the stability threshold is away from the bulk spectrum, so the spectral density remains as Wigner’s semicircle. On the other hand, if the interaction strength is greater than the critical value, the spectral density becomes

$$\rho_{\text{hard}}(\lambda) = \frac{1}{c^2} \frac{\sqrt{\lambda + l^* - 1}}{2\pi\sqrt{1-\lambda}} (l^* - 2\lambda), \qquad (17)$$

where the width of the support

$$l^* = \frac{2}{3} \left( 1 + \sqrt{1 + 6c^2} \right). \qquad (18)$$

Building on these existing results for the spectral distribution, we now can ask for an interpretation of the emergent time scales. As shown in Fig. 2(b), the stability threshold acts as a wall pushing the eigenvalues to pile up. More precisely, near the stability threshold $\lambda = 1$ we have $\rho_{hard}(\lambda) \sim (1 - \lambda)^{-1/2}$, which [by the same argument as in Eq. (14)] indicates that the longest time scale increases as system size with $\tau_{max} \sim N^2$.

The autocorrelation function also is dominated by the eigenvalues close to the stability threshold. Following steps detailed in Appendix B, we can write the resulting autocorrelation coefficient as

$$R_{\text{hard}}(t) = \frac{C(t)}{C(0)} = 1 - \sqrt{\pi}(\varepsilon t)^{1/2} + \varepsilon t + \mathcal{O}((\varepsilon t)^2), \qquad (19)$$

where $\varepsilon \sim N^{-2}$ is an infrared cut-off introduced to compute the correlation function in the thermodynamic limit, and the correlation time as

$$\tau_{\text{corr}} \sim \varepsilon^{-1} \sim N^2. \qquad (20)$$

We see that for supercritical systems, both the longest time scale $\tau_{max}$ and the correlation time $\tau_{corr}$ increase as a power of the system size; the rate is even faster than the system at criticality. In fact, there are divergently many slow modes. Meanwhile, the interaction strength can undertake a range of values, as long as they are greater
than a certain threshold. Conceptually, with the two-dimensional electron gas analogy, the hard threshold imposes a semi-infinite well in addition to the quadratic potential and logarithmic repulsion. Although divergently many eigenvalues cluster near the stability threshold, the repulsion among the eigenvalues indeed pushes the smallest eigenvalue further away from the stability threshold, introducing shorter timescales to the system.

Nonetheless, we cannot quite claim that the problem is solved. First, in the supercritical phase, the correlation function does not decay as a power law. Instead, the correlation function stays at 1 for a time period $\tau_{corr} \sim \tau_{max}$, and then decays exponentially. This means the system has a single long time scale, rather than a continuous spectrum of slow modes. Second, in order for a system to impose a hard constraint on its stability, it needs to measure its stability. Naively, checking for stability, especially in the presence of slow modes, requires access to infinitely long measuring times; implementing a sharp threshold may also be challenging.

### C. Model 3: Constraining mean-square activity

While it can be difficult to check for stability, it is much easier to imagine checking the overall level of activity in the network. One can even think about mechanisms that would couple indirectly to activity, such as the metabolic load. If the total activity is larger than some target level, the system might be veering toward instability, and there could be feedback mechanisms to reduce the overall strength of connections. Regulation of this qualitative form is known to occur in the brain, and is termed synaptic scaling [43–45]; this is hypothesized to play an important role in maintaining persistent neural activities [30, 31]. In this section, we construct the least structured distribution $P(M)$ that is consistent with a fixed mean (square) level of activity, which we can think of as a soft threshold on stability, and derive the density of eigenvalues that follow from this distribution. In the following section we discuss possible mechanisms for a system to generate matrices $M$, dynamically, out of this

ensemble.

#### 1. The spectral density

It is useful to remember that the GOE, Eq. (11), can be seen as the maximum entropy distribution of matrices consistent with some fixed variance of the matrix elements $M_{ij}$ [46, 47]. If we want to add a constraint, we can stay within the maximum entropy framework, and in this way we isolate the generic consequences of this constraint: we are constructing the least structured ensemble of networks that satisfies the added condition.

Recall that if we want to constrain the mean values of several functions $f_\mu(M)$, the maximum entropy distribution has the form

$$P(M) \propto \exp \left[ - \sum_{\mu} g_{\mu} f_{\mu}(M) \right], \quad (21)$$

where the $g_\mu$'s are the Lagrange multipliers used to enforce the constraints [46, 47]. In our case, we are interested in the mean-square value of the individual matrix elements, and the mean-square value of the activity variables $x_i$. Because our basic model Eqs. (1) and (2) predicts that

$$\mu = \frac{1}{N} \sum_{i} \langle x_i^2 \rangle = \frac{1}{N} \sum_{i} \frac{1}{1 - \lambda_i}, \quad (22)$$

the relevant maximum entropy model becomes

$$P_{\text{soft}}(M) \propto \exp \left[ -\frac{N}{2c^2} \text{Tr} M^\top M - N\xi \sum_i \frac{1}{1-\lambda_i} \right], \quad (23)$$

where $\xi$ is the Lagrange multiplier constraining the norm activity $\mu$. The scaling of $N\xi$ ensures a well-defined thermodynamic limit. See Fig. 2(d) for schematics.

Again, this distribution is invariant to orthogonal transformation. After the integration over the rotation matrices, we have

$$P_{\text{soft}}(\{\lambda_i\}) \propto \exp \left[ -\frac{N}{2c^2} \sum_i \lambda_i^2 + \frac{1}{2} \sum_{j \neq k} \ln |\lambda_j - \lambda_k| - N\xi \sum_i \frac{1}{1-\lambda_i} \right]. \quad (24)$$

Mathematically, this is the Boltzmann distribution for two-dimensional electron gas, with a soft wall local potential. Luckily, the same arguments that yield the exact density of eigenvalues in the Gaussian Orthogonal Ensemble also work here (see Appendix A), and we find

$$\rho_{\text{soft}}(\lambda) = \frac{1}{\pi \sqrt{(\lambda - 1 + g_0 + l)(1 - g_0 - \lambda)}} \left[ 1 + \frac{l^2}{8c^2} + \left(1 - g_0 - \frac{1}{2}l\right) \frac{\lambda}{c^2} - \frac{\lambda^2}{c^2} + \frac{\xi (2g_0 - 2g_0^2 + l - 2g_0l) - (2g_0 + l)\lambda}{\sqrt{g_0(g_0 + l)(\lambda - 1)^2}} \right], \quad (25)$$

where the gap size $g_0$ and the width of the support $l$ are

solved numerically by setting the spectral density at the
two ends of the support to zero.

To our surprise, unlike for GOE with hard stability
constraints as in Model 2, here we find a finite gap for all
$\xi > 0$ (see Appendix C for calculation). Intuitively, this
is due to the soft threshold acting as a repelling force
on the eigenvalues to push them away from the stabil-
ity threshold. This finite gap means that there exists a
maximum time scale even when the system is infinitely
large. This upper limit of longest time scale depends on
the Lagrange multiplier $\xi$ and the interaction strength $c$.
Because the Lagrange multiplier $\xi$ is used to constrain the
averaged norm activity $\mu$, the maximum time scale is set
by the allowed norm of the activity $\mu$, measured in units
of expected norm for independent neurons (Fig. 3(e)). As
we explain below, the greater dynamic range the system
leads to longer time scales.

The dependence of the gap on the Lagrange multiplier $\xi$ is shown in Fig. 2(c) at each of several fixed values of the interaction strength $c$. In the limit of $\xi \to 0$, we recover the hard wall case, where there is a phase transition at $c_{crit} = 1/\sqrt{2}$. This phase transition remains present for $0 < \xi \ll 1$, with the eigenvalues close to the stability threshold pushed into the bulk spectrum (see details in Section III C 2 and Appendix D). For large $\xi > 1$, the entire spectrum is pushed away from the wall. A closer look at the longest time scale $\tau_{max}$ vs. Lagrange multiplier in Fig. 3(a) confirms that amplification of time scales occurs only when $\xi < 1$, corresponding to an amplification of mean-square activity $\mu$.

## 2. The scaling of time scales in three phases

We now discuss the dependence of time scales on the
interaction strength $c$. In contrast to the ensemble with a
hard stability threshold, we find a finite gap for all values
of interaction strength $c$, but the scaling of time scales
vs. the Lagrange multiplier (and hence the mean-square
activity) is different in the different phases. Details of
the calculation can be found in Appendix D.

For the subcritical and critical phases, the results are
as expected from the hard wall case. In the subcrit-
ical phase, as $\xi \to 0$, the spectral distribution con-
verges smoothly to the familiar semicircle. The time
scales and the mean-square activity both approach con-
stants (Fig. 3(a,b,e), red curves). On the other hand,
when $c = c_{crit}$, we find that the longest time scale
grows as $\tau_{max} \sim \xi^{-2/5}$, and the correlation time scale
$\tau_{corr} \sim \xi^{-1/5}$, i.e. both time scales can be large if $\xi$
is small enough (Fig. 3(a, b), blue curves). Meanwhile,
the norm activity $\mu$ asymptotically approaches a con-
stant value of $\mu = 2$ (Fig. 3(c, d), blue curves). This
suggests that if a system is poised at criticality, then the
system can exhibit long time scales, even when the dy-
namic range of individual components is well controlled.
The autocorrelation function exhibits a power law-like
decay, as expected; see the blue curve in Fig. 3(g).

The most interesting case is the supercritical phase,

<figure><img src="image_5.png" /><figcaption>FIG. 3: Mean field results for Model 3 constraining the mean-square activity. (a) The longest time scale $\tau_{max}$ vs. the Lagrange multiplier $\xi$ which constrains the mean-square activity, for interaction strength $c = 0.6$ (subcritical, red), $c = c_{crit} = 1/\sqrt{2}$ (critical, blue), $c = 0.8$ (supercritical, black). The dashed lines bound the scaling behavior at small $\xi$ (see Appendix D). (b) Same as (a) but for the correlation time scale $\tau_{corr}$. (c, d) Same as (a, b), but plotting against the average norm activity $\mu$. At $c = c_{crit}$, the norm activity asymptotically approaches $\mu = 2$. (e) Average norm activity $\mu$ as a function of the Lagrange multiplier $\xi$. (f) The exact scaling and the amplification of time constant depending on whether the system is subcritical (with interaction strength $c < 1/\sqrt{2}$), critical ($c = 1/\sqrt{2}$), or supercritical. The different shade indicates different values of the Lagrange multiplier $\xi$: from lightest to darkest, $\xi = 10^{-1}, 10^{-4}, 10^{-7}, 10^{-10}$, respectively (see Fig. 4(e) legend). (g) Normalized correlation function $R(t)$ computed at $\xi = 10^{-10}$. For each curve from left to right, the interaction strength is $c = 0.2, c = 0.6$ (red), $c = c_{crit}$ (blue), $c = 1, c = 2$ and $c = 3$ (black), respectively.</figcaption></figure>

where the interaction strength $c > c_{crit} = 1/\sqrt{2}$. As
$\xi \to 0$, the spectrum does not converge to the spec-
trum with the hard constraint. We find that both time scales and the norm activity increase as power laws of the Lagrange multiplier $\xi$ (Fig. 3(a-d), black curves), with $\tau_{max} \approx 3\tau_{corr} \sim \xi^{-2/3}$, and $\mu \sim \xi^{-1/3}$. This implies that the time scales grow as a power of the allowed dynamic range of the system, although not with the size of the system. The question of whether the resulting time scales are “long” then becomes more subtle. Quantitatively, we see from Fig. 3(c), that if the system has an allowed dynamic range just $10 \times$ that of independent neurons, the system can generate time scales $\tau_{max}$ almost $10^4 \times$ longer than the time scale of isolated neurons. Similar amplification of the correlation time can been seen in Fig. 3(d).

Interestingly, once the system is in the supercritical phase, the ratio of amplification has only a small dependence on the interaction strength $c$ (Fig. 3(f)). Intuitively, an increasing interaction strength $c$ implies that without constraints more modes will be unstable, while with constraints more modes concentrate near the stability threshold, but the entire support of the spectrum also expands (see Fig. 2(d) and Appendix A Eq. (A7)), so the density of slow modes and their distance to the stability threshold remain similar. This is perhaps another indication for long time scales without fine tuning when the system uses its norm activity to regulate its connection matrix.

We note that, although both the critical phase and the supercritical phase can reach time scales that are as long as the dynamic range allows, there are significant differences between the two phases. One difference is that in the critical phase, locally the dynamic range for each neuron can remain finite, while for the supercritical phase, the variance of activity for individual neurons can be much greater. Moreover, as shown by Fig. 3(g), systems in the supercritical phase are dominated by a single slow mode, rather than by a continuous spectrum of slow modes. While the autocorrelation function decays as a power law in the critical phase, in the supercritical phase, it holds at $R(t) = 1$ for a much longer time compared to the subcritical case, but then decays exponentially. While a single long time scale can be achieved without fine tuning, it seems that a continuous spectrum of long time scales is much more challenging.

### 3. Finite Size Effects

If we want these ideas to be relevant to real biological systems, we need to understand what happens at finite $N$. We investigate this numerically using direct Monte-Carlo sampling of the (joint) eigenvalue distribution. As the system size grows, the time scales $\tau_{max}$ and $\tau_{corr}$ also increase, up to the upper limit given by the mean field results; see Fig. 4(a, b). The relation between the timescales and the average norm activity does not depend on the system size (Fig. 4(c, d)). Interestingly, the scaling exponent $\alpha$ for the gap difference,

$$ \Delta g_0(N) \equiv g_0(N) - g_0 \sim N^{\alpha}, \quad (26) $$

depends on the Lagrange multiplier $\xi$ (Fig. 4(e,f)). This is supported by the analytical argument that the scaling interpolates between two limiting cases: for small $\xi$, the gap scales as $\Delta g_0 \sim N^{-2}$, with the same scaling factor as in the case with the hard threshold; for large $\xi$, the gap scales as $\Delta g_0 \sim N^{-2/3}$, which has the same scaling factor as in the Gaussian ensemble without any additional constraint. In any case, thousands of neurons will be well described by the $N \to \infty$ limit.

<figure><img src="image_6.png" /><figcaption>FIG. 4: Finite size effects on the time scales $\tau_{max}$ and $\tau_{corr}$ for Model 3, super-critical phase, with interaction strength $c = 1$. (a) The longest time scale $\tau_{max}$ as a function of the Lagrange multiplier $\xi$. Results from direct Monte Carlo sampling of the eigenvalue distribution are plotted together with the mean field results. (b) Same as panel (a) but for the correlation time scale $\tau_{corr}$. (c) The longest timescale $\tau_{max}$ as a function of the norm activity $\mu$. The difference among different-sized systems is small. (d) Same as (c) but for the correlation time scale $\tau_{corr}$. (e) The gap difference $\Delta g_0(N) \equiv g_0(N) - g_0$ vs. system size $N$. (f) The exponent $\alpha$ that explains the convergence of $\tau_{max}$ when system size increases depends on the Lagrange multiplier $\xi$, interpolating between $\alpha = 2$ and $\alpha = 2/3$. The red curve shows a sigmoidal fit.</figcaption></figure>
### 4. Distribution of matrix elements

Now that we have examples of network ensembles that exhibit long time scales, we need to go back and check what these ensembles predict for the distribution of individual matrix elements. In particular, because we did not constrain the self interaction $M_{ii}$ to be 0, we want to check whether the long time scales emerge as a collective behavior of the network, or trivially from an effective increase of the intrinsic time scales for individual neurons. A similar question arises in real networks, where there have been debates about the importance of feedback within single neurons vs. the network dynamics in maintaining long time scales; see Ref. [12] for a review. We confirm that, at least in our models, this is not an issue: the constraint on the norm activity, in fact, pushes the average eigenvalues to be negative, and hence the effective self interaction for individual neurons actually leads to a shorter intrinsic time scale. We can impose additional constraints on the distribution of matrix elements, such that $\langle M_{ii} \rangle = 0$; for this ensemble, we can again solve for the spectral distribution, and we find that scaling behaviors described above do not change (figures not shown).

### D. Dynamic tuning

So far, we have established that a distribution constraining the norm activity can lead generically to long time scales, but we haven't really found a mechanism for implementing this idea. But if we can write the distribution of the interaction matrices $M$ in the form a Boltzmann distribution, we know that we can sample this distribution by allowing the matrix elements be dynamical variables undergoing Brownian motion in the effective potential. We will see that this sort of dynamics is closely related to previous work on self-tuning to criticality [29], and we can interpret the dynamics as implementing familiar ideas about synaptic dynamics, such as Hebbian learning and metaplasticity.

We can rewrite our model in Eq. (23) as

$$ P_{\text{soft}}(M) = \frac{1}{Z} \exp \left[ -\frac{N}{2c^2} \text{Tr} M^{\text{T}} M - N\xi \sum_i \frac{1}{1-\lambda_i} \right] \\ = \frac{1}{Z} \exp [-V(M)/T], \quad (27) $$

$$ V(M) = \frac{1}{2} \text{Tr} M^{\text{T}} M + c^2 \xi \text{Tr}(1 - M)^{-1}, \quad (28) $$

with a temperature $T = c^2/N$. The matrix $M$ will be drawn from the distribution $P(M)$, as $M$ itself performs Brownian motion or Langevin dynamics in the potential

$V(M)$:

$$ \begin{aligned} \tau_M \dot{M} &= -\frac{\partial V(M)}{\partial M} + \zeta(t) \\ &= -M - c^2 \xi (\mathbf{1} - M)^{-2} + \zeta(t), \end{aligned} \quad (29) $$

where the noise has zero mean, $\langle \zeta \rangle = \mathbf{0}$, and is independent for each matrix element,

$$ \langle \zeta_{ij}(t) \zeta_{kl}(t') \rangle = 2T \tau_M \delta_{ik} \delta_{jl} \delta(t-t'). \quad (30) $$

It is useful to remember that, in steady state, our dynamical model for the $\{x_i\}$, Eqs. (1) and (2), predicts that

$$ \langle x_i x_j \rangle = [(\mathbf{1} - M)^{-1}]_{ij}. \quad (31) $$

This means that we can rewrite the Langevin dynamics of $M$, element by element, as

$$ \begin{aligned} \tau_M \dot{M}_{ij} &= -M_{ij} - c^2 \xi [(\mathbf{1} - M)^{-2}]_{ij} + \zeta_{ij}(t) \\ &= -M_{ij} - c^2 \xi \langle x_i x_k \rangle \langle x_k x_j \rangle + \zeta_{ij}(t). \end{aligned} $$

Because the $x_i$ are Gaussian, we have

$$ \begin{aligned} \sum_k \langle x_i x_k \rangle \langle x_k x_j \rangle &= \frac{1}{2} \sum_k (\langle x_i x_k x_k x_j \rangle - \langle x_i x_j \rangle \langle x_k x_k \rangle) \\ &= \frac{1}{2} \left\langle x_i x_j \left( \sum_k x_k^2 - \sum_k \langle x_k^2 \rangle \right) \right\rangle. \end{aligned} \quad (32) $$

We now imagine that the dynamics of $M$ is sufficiently slow that we can replace averages by instantaneous values, and let the dynamics of $M$ do the averaging for us. In this approximation we have

$$ \tau_M \dot{M}_{ij} = -M_{ij} - \frac{1}{2} c^2 \xi x_i x_j \left( \sum_k x_k^2 - \theta \right) + \zeta_{ij}(t), \quad (33) $$

where the threshold $\theta = \sum_k \langle x_k^2 \rangle$.

The terms in this Langevin dynamics have a natural neuroscience interpretation. First, the connection strength decays with an overall time constant $\tau_M$. Second, the synaptic connection $M_{ij}$ is driven by the correlation between pre- and post-synaptic activity, $\sim x_i x_j$, as in Hebbian learning [48, 49]. The response to correlations is modulated depending on whether the global neural activity is greater or less than a threshold value. If the network is highly active, then the connection between neurons with correlated activity will decrease, i.e. the dynamics are anti-Hebbian; if the overall network is quiet, the dynamics are Hebbian.

We still have the problem of setting the threshold $\theta$. The first option is to set the threshold fixed using a constant, $\tilde{\theta}$. Ideally, for the dynamics to generate samples out of the correct $P(M)$, we need

$$ \tilde{\theta} = \theta^* \equiv \langle x^\text{T} x \rangle_{\text{s.s.}} = N \mu(c, \xi), \quad (34) $$
<figure><img src="image_7.png" /><figcaption>FIG. 5: (a) The spectral distribution for connection matrix drawn from the dynamics with fixed threshold approaches the static distribution as the time constant $\tau_M$ increases. (b) Dynamics of representative eigenvalues with the optimal fixed threshold $\tilde{\theta} = \theta^*$, with $N = 64$, and time constant $\tau_M = 1$ and 100. If the connection matrix updates too fast (top), i.e. $\tau_M$ is too small, the system exhibits quasi-periodic oscillations, and does not reach a steady state distribution. (c) The expected eigenvalues for the dynamics with fixed threshold $\tilde{\theta}$, and for the ones with sliding threshold, with $N = 32$, $c = 1$, $\xi = 2^{-5}$ and $\tau_M = 1000$. (d) Example traces for the eigenvalues vs. time, for fixed thresholds at different threshold values (left), and for sliding thresholds at different updating timescales $\tau_\theta$. (e) Heatmap for expected $\lambda_{max}$ as a function of dynamical model parameters, i.e. the time constant $\tau_M$, the fixed threshold $\tilde{\theta}$ (left), and the sliding threshold time scale $\tau_\theta$ (right). At equal $\tau_M$, the dynamics with a sliding threshold has a much larger range of parameter $\tau_\theta$ for the system to exhibit slow modes ($\lambda_{max} \lesssim 1$), while the dynamics with fixed threshold requires fine-tuning of its parameter $\tilde{\theta}$.</figcaption></figure>

where as above $\mu$ is the mean activity whose value is enforced by the Lagrange multiplier $\xi$. This means that the fixed threshold $\tilde{\theta}$ needs to be tuned in relation to $\xi$, and it is challenging to have a mechanism that does this directly, and just pushes the fine tuning problem back one step. Importantly, if $\tilde{\theta} = \theta^*$ then the steady state spectral density of the connection matrix approaches the desired equilibrium distribution as the update time constant increases (Fig. 5(a)). As expected, if the connection matrix updates too fast, i.e. $\tau_M$ is too small, the system exhibits quasi-periodic oscillations, and does not reach a steady state distribution. In contrast, long $\tau_M$ leads to the adiabatic approximation for the steady state distribution (Fig. 5(b)). Nonetheless, if $\tilde{\theta}$ deviates from $\theta^*$, then the steady distribution does not have slow modes. As shown in Fig. 5(c), if the threshold is too small, then the entire spectrum is shifted away from the stability threshold, and the system no longer exhibits long time scales; if the threshold is too large, then the largest eigenvalue oscillates around the stability threshold $\lambda = 1$, and a typical connection matrix drawn from the steady state distribution is unstable (Fig. 5(d)).

But we can once again relieve the fine tuning problem

by promoting $\theta$ to a dynamical variable,

$$ \tau_{\theta} \dot{\theta} = \sum_{k} x_{k}^{2} - \theta, \qquad (35) $$

which we can think of as a sliding threshold in the spirit of models for metaplasticity [50]. We pay the price of introducing yet another new time scale, $\tau_\theta$, but in Fig. 5(d,e) we see that this can vary over at least three orders of magnitude without significantly changing the spectral density of the eigenvalues.

To see whether the sliding threshold really works, we can compare results where $\theta$ is fixed to those where it changes dynamically; we follow the mean value of $\lambda_{max}$ as an indicator of system performance. We choose parameters in the supercritical phase, specifically $c = 1$ and $\xi = 2^{-5}$, and study a system with $N = 32$. Figure 5(e) shows that with fixed threshold, even in the adiabatic limit where $\tau_M \gg 1$, there is only a measure-zero range for the fixed threshold $\tilde{\theta}$ such that $\lambda_{max}$ is very close to, but smaller than 1, and hence requires fine-tuning of its parameters. In contrast, for the dynamics with sliding threshold, at $\tau_M \gg 1$ there is a large range of values for the time constant $\tau_\theta$ such that the system hovers just below instability, generating a long time scale.
IV. DISCUSSION

Living systems can generate behaviors on time scales that are much longer than the typical time scales of their component parts, and in some cases correlations in these behaviors decay as an approximate power-law, suggesting a continuous spectrum of slow modes. In order to understand how surprised we should be by these behaviors, it is essential to ask whether there exist biologically plausible dynamical systems that can generate these long time scales. We have constructed a mechanism for a broad class of models for interacting living systems–symmetric linear dynamical models–to reach long time scales with the least stringent fine-tuning condition: when the interaction strength of the connection matrix is large enough, imposing a global constraint on the stability of the system leads to divergently many slow modes. To impose a biologically plausible mechanism for living systems, we constrain the averaged norm activity as a proxy for global stability; in this case, the time scales for the slow modes are set by the allowed dynamic range of the system.

Further, we showed that ensembles of interaction ma-
trix *M* with long time scales can be achieved dynamically
by updating *M* with a sliding threshold for the norm
activity, a mechanism that resembles homeostatic synap-
tic plasticity in neuronal networks, particularly through
mechanisms such as synaptic scaling and switching be-
tween Hebbian and anti-Hebbian learning by threshold-
ing the average neural firing rate [43, 44]. Biologically,
it has been suggested that global modulation of synap-
tic strengths across a spatially extended network can
be effected by activity-dependent secretion of molecules
such as brain-derived neurotrophic factor, tumor necrosis
factor-α, and retinoic acid [51, 52]. From the modeling
perspective, we can compare the Langevin dynamics of
our system to the BCM theory of metaplasticity in neu-
ral networks, in that both models involve a combination
of Hebbian learning and a threshold on the neural activ-
ity [50]. Both dynamics have more controlled dynamics
when the activity threshold is chosen to be sliding com-
pared to fixed, but there are two key differences. First,
the BCM theory imposes a threshold on the activity of
locally connected neurons, while here the threshold is on
the overall neural activity. Second, the BCM dynamics is
Hebbian when the postsynaptic activities are larger than
the threshold, and anti-Hebbian otherwise, which is the
opposite of the dynamics for our system. It is interest-
ing that in some other models for homeostasis, plastic-
ity requires the activity detection mechanism to be fast
(τθ/τM ≪ 1) for the system to be stable [53, 54], which
we do not observe for our system.

Importantly, the slow modes achieved through con-
straining norm activity typically lead to exponentially
decaying correlations; only when the interaction strength
of the matrix is at a critical value do we find power-law
decays. This suggests that a continuous range of slow
modes is more difficult to achieve than a single long time
scale. A natural follow-up question is whether there exist

mechanisms which can tune the system to criticality in a
self-organized way, for example by coupling the interac-
tion strength to the averaged norm activity.

To achieve long time scales a high dimensional dynam-
ical system requires some degree of fine tuning: from the
most to the least stringent, examples include setting in-
dividual elements of the connection matrix by coupling
the system to the environment [10], or choosing a par-
ticular network architecture such as a ring or grid net-
work [24–26]. All of those network designs require tuning
of either the structure (i.e. whether certain neurons are
connected) or the interaction strength (i.e. how strong
is the connection). Comparing to other models with a
focus on tuning the network structure, our work focuses
on the distribution of interaction strength. Starting from
all-to-all connected networks allow us to be as general as
possible, which is especially useful when one considers
a dynamical process that leads to such functional net-
works. Possible future work includes pruning the network
where the interaction strength is small, which will allow a
comparison between the emergent network structure and
other works which start with a fixed structure.

Both for simplicity and to understand the most ba-
sic picture, we have been focusing on linear networks
with symmetric connections. Realistically, many biolog-
ical networks are asymmetric, which gives rise to more
complex and even chaotic dynamics [55, 56]. In the asym-
metric case, the eigenvalue spectrum for the Gaussian
ensemble is well known (uniform distribution inside a
unit circle) [57, 58], but a similar global constraint on
the norm activity leads to a dependence on the overlap
among the left and right eigenvectors (see Appendix E
and Ref. [59, 60]). In particular, the matrix distribution
can no longer be separated into the product of eigenvalues
and eigenvectors, and it is difficult to solve for the spec-
tral distribution analytically. Two new features emerge
in the asymmetric case. First, the time scales given by
real eigenvalues vs. complex eigenvalues may be differ-
ent, leading to more (or less) dominant oscillatory slow
modes in large systems [61]. Second, asymmetric con-
nection matrices can lead to complicated transient activ-
ity when the system is perturbed, with the time scales
mostly dominated by the eigenvector overlaps, and can
be very different from the time scales given by the eigen-
values [62]. In the limit of strong asymmetry, the net-
work is organized into a feed-forward line, which with
fine tuning of interactions, information can be held for
a time that is extensive in system size; see examples in
Refs. [24] and [63]. It will be interesting to check whether
systems can store information in these transients without
fine-tuning the structure of the network.

While we compare this dynamical updating mechanism
to synaptic updating mechanisms in neuroscience, note
that the methods we propose is conceptual and can be
applied generally to models of interacting systems, and
is not limited to neurons. The system we study can be
extended to consider more specific ensembles for partic-
ular biological systems. For example, real neural net-
works have inhibitory and excitatory neurons, so that elements belonging to the same column need to share the same sign, and the resulting spectral distribution has been shown to differ from the unit sphere [64]; more generally recent work explores how structured connectivity can lead to new but still universal dynamics [65]. Another limitation of our work is that it only considers linear dynamics, or only dynamical systems where all fixed points to be equally likely. In contrast, non-linear dynamics such as the Lotka-Volterra model in ecology [66], discrete spiking dynamics in neuronal nets [67], and the gating neural network in machine learning [68, 69] have been shown to drive systems to non-generic, marginally stable fixed points, around which there exists an extensive number of slow directions for the dynamics. In summary, we believe the issue of whether a continuum of slow modes can arise generically in neural networks remains open, but we hope that our study of very simple models has helped to clarify this question.

## Code availability

Custom codes were developed in C to implement standard Monte Carlo sampling methods and stochastic dynamics. The codes are available at https://doi.org/10.5281/zenodo.12103680.

## Acknowledgments

We thank Hanrong Chen and Kamesh Krishnamurthy for fruitful discussions. This work was supported in part by the National Science Foundation, through the Center for the Physics of Biological Function (PHY-1734030) and Grant PHY-1607612, and by the National Institutes of Health through Grants NS104889 and R01EB026943-01.

[1] W. Bialek and J.W. Shaevitz. Long timescales, individual differences, and scale invariance in animal behavior. *Phys Rev Lett*, 132(4):048401, 2024.

[2] A. Flack, M. Nagy, W. Fiedler, I.D. Couzin, and M. Wikelski. From local collective behavior to global migratory patterns in white storks. *Science*, 360(6391):911-914, 2018.

[3] A.C. Costa, T. Ahamed, and G.J. Stephens. Adaptive, locally linear models of complex dynamics. *Proc Natl Acad Sci (USA)*, 116(5):1501–1510, 2019.

[4] G. Sridhar, M. Vergassola, J. Marques, M.B. Orger, A.C. Costa, and C. Wyart. Uncovering multiscale structure in the variability of larval zebrafish navigation. *arXiv preprint arXiv:2405.17143*, 2024.

[5] A.C. Costa, G. Sridhar, C. Wyart, and M. Vergassola. Fluctuating landscapes and heavy tails in animal behavior. *PRX Life*, 2(2):023001, 2024.

[6] M.C. Morrell, A.J. Sederberg, and I. Nemenman. Latent dynamical variables produce signatures of spatiotemporal criticality in large biological systems. *Phys Rev Lett*, 126(11):118302, 2021.

[7] D.J. Schwab, I. Nemenman, and P. Mehta. Zipfs law and criticality in multivariate data without fine-tuning. *Phys Rev Lett*, 113(6):068102, 2014.

[8] E. Aksay, G. Gamkrelidze, H.S. Seung, R. Baker, and D.W. Tank. In vivo intracellular recording and perturbation of persistent activity in a neural integrator. *Nat Neurosci*, 4(2):184–193, 2001.

[9] C.D. Brody, A. Hernández, A. Zainos, and R. Romo. Timing and neural encoding of somatosensory parametric working memory in macaque prefrontal cortex. *Cereb Cortex*, 13(11):1196–1207, 2003.

[10] G. Major, R. Baker, E. Aksay, B. Mensh, H.S. Seung, and D.W. Tank. Plasticity and tuning by visual feedback of the stability of a neural integrator. *Proc Natl Acad Sci (USA)*, 101(20):7739–7744, 2004.

[11] G. Major, R. Baker, E. Aksay, H.S. Seung, and D.W. Tank. Plasticity and tuning of the time course of analog persistent firing in a neural integrator. *Proc Natl Acad Sci (USA)*, 101(20):7745–7750, 2004.

[12] G. Major and D.W. Tank. Persistent neural activity: prevalence and mechanisms. *Curr Opin in Neurobiol*, 14(6):675–684, 2004.

[13] R. Srimal and C.E. Curtis. Persistent neural activity during the maintenance of spatial position in working memory. *Neuroimage*, 39(1):455–468, 2008.

[14] J.D. Murray, A. Bernacchia, D.J. Freedman, R. Romo, J.D. Wallis, X. Cai, C. Padoa-Schioppa, T. Pasternak, H. Seo, D. Lee, et al. A hierarchy of intrinsic timescales across primate cortex. *Nat Neurosci*, 17(12):1661–1663, 2014.

[15] A. Bernacchia, H. Seo, D. Lee, and X.-J. Wang. A reservoir of time constants for memory traces in cortical neurons. *Nat Neurosci*, 14(3):366–372, 2011.

[16] M. Rigotti, D.B.D. Rubin, X.-J. Wang, and S. Fusi. Internal representation of task rules by recurrent dynamics: the importance of the diversity of neural responses. *Front Comput Neurosci*, 4:24, 2010.

[17] A. Miri, K. Daie, A.B. Arrenberg, H. Baier, E. Aksay, and D.W. Tank. Spatial gradients and multidimensional dynamics in a neural integrator circuit. *Nat Neurosci*, 14(9):1150–1159, 2011.

[18] L. Meshulam, J.L. Gauthier, C.D. Brody, D.W. Tank, and W. Bialek. Coarse graining, fixed points, and scaling in a large population of neurons. *Phys Rev Lett*, 123(17):178103, 2019.

[19] R. Chaudhuri, K. Knoblauch, M.-A. Gariel, H. Kennedy, and X.-J. Wang. A large-scale circuit mechanism for hierarchical dynamical processing in the primate cortex. *Neuron*, 88(2):419–431, 2015.

[20] S. Fusi, P.J. Drew, and L.F. Abbott. Cascade models of synaptically stored memories. *Neuron*, 45(4):599–611, 2005.

[21] G. Tanaka, T. Matsumori, H. Yoshida, and K. Aihara. Reservoir computing with diverse timescales for prediction of multiscale dynamics. *Phys Rev Res*, 4(3):L032014, 2022.

[22] H.S. Seung. How the brain keeps the eyes still. *Proc Natl Acad Sci (USA)*, 101(20):7745–7750, 2004.
Acad Sci (USA), 93(23):13339-13344, 1996.

[23] C.D. Brody, R. Romo, and A. Kepecs. Basic mechanisms for graded persistent activity: discrete attractors, continuous attractors, and dynamic representations. *Curr Opin Neurobiol*, 13(2):204-211, 2003.

[24] M.S. Goldman. Memory without feedback in a neural network. *Neuron*, 61(4):621-634, 2009.

[25] Y. Burak and I.R. Fiete. Fundamental limits on persistent activity in networks of noisy neurons. *Proc Natl Acad Sci (USA)*, 109(43):17645-17650, 2012.

[26] Y.-L. Shi, R. Zeraati, A. Levina, and T.A. Engel. Spatial and temporal correlations in neural networks with structured connectivity. *Phys Rev Res*, 5(1):013005, 2023.

[27] A. Litwin-Kumar and B. Doiron. Slow dynamics and high variability in balanced cortical networks with clustered connections. *Nat Neurosci*, 15(11):1498-1505, 2012.

[28] M. Stern, N. Istrate, and L. Mazzucato. A reservoir of timescales emerges in recurrent circuits with heterogeneous neural assemblies. *Elife*, 12:e86552, 2023.

[29] M.O. Magnasco, O. Piro, and G.A. Cecchi. Self-tuned critical anti-Hebbian networks. *Phys Rev Lett*, 102(25):258102, 2009.

[30] A. Renart, P. Song, and X.-J. Wang. Robust spatial working memory through homeostatic synaptic scaling in heterogeneous cortical networks. *Neuron*, 38(3):473-485, 2003.

[31] C. Tetzlaff, C. Kolodziejski, M. Timme, M. Tsodyks, and F. Wörgötter. Synaptic scaling enables dynamically distinct short- and long-term memory formation. *PLoS Comput Biol*, 9(10), 2013.

[32] W. Bialek. On the dimensionality of behavior. *Proc Natl Acad Sci (USA)*, 119(18):e2021860119, 2022.

[33] G. Livan, M. Novaes, and P. Vivo. *Introduction to Random Matrices: Theory and Practice*, volume 26. Springer International Publishing, 2018.

[34] R. Marino. Number statistics in random matrices and applications to quantum systems. PhD thesis, Université Paris-Saclay, 2016.

[35] E.P. Wigner. On the statistical distribution of the widths and spacings of nuclear resonance levels. *Math Proc Cambridge*, 47(4):790798, 1951.

[36] R.M. May. Will a large complex system be stable? Nature, 238(5364):413-414, 1972.

[37] J.-P. Bouchaud and M. Potters. *Theory of Financial Risk and Derivative Pricing: From Statistical Physics to Risk Management*. Cambridge University Press, 2003.

[38] F.J. Dyson. The threefold way: Algebraic structure of symmetry groups and ensembles in quantum mechanics. *J Math Phys*, 3(6):1199-1215, 1962.

[39] F.J. Dyson. A Brownian motion model for the eigenvalues of a random matrix. *J Math Phys*, 3(6):1191-1198, 1962.

[40] D.S. Dean and S.N. Majumdar. Large deviations of extreme eigenvalues of random matrices. *Phys Rev Lett*, 97(16):160201, 2006.

[41] D.S. Dean and S.N. Majumdar. Extreme value statistics of eigenvalues of Gaussian random matrices. *Phys Rev E*, 77(4):041108, 2008.

[42] S.N. Majumdar and G. Schehr. Top eigenvalue of a random matrix: large deviations and third order phase transition. *J Stat Mech*, 2014(1):P01012, 2014.

[43] G.G. Turrigiano, K.R. Leslie, N.S. Desai, L.C. Rutherford, and S.B. Nelson. Activity-dependent scaling of quantal amplitude in neocortical neurons. *Nature*,

391(6670):892-896, 1998.

[44] G.G. Turrigiano and S.B. Nelson. Homeostatic plasticity in the developing nervous system. *Nat Rev Neurosci*, 5(2):97-107, 2004.

[45] L.F. Abbott and S.B. Nelson. Synaptic plasticity: taming the beast. *Nat Neurosci*, 3(S11):1178-1183, 2000.

[46] E.T. Jaynes. Information theory and statistical mechanics. *Phys Rev*, 106(4):620, 1957.

[47] S. Pressé, K. Ghosh, J. Lee, and K.A. Dill. Principles of maximum entropy and maximum caliber in statistical physics. *Rev Mod Phys*, 85(3):1115-1141, 2013.

[48] D.O. Hebb. *The Organization of Behavior: A Neuropsychological Theory*. J. Wiley; Chapman & Hall, 1949.

[49] J.C. Magee and C. Grienberger. Synaptic plasticity forms and functions. *Annu Rev Neurosci*, 43:95-117, 2020.

[50] E.L. Bienenstock, L.N. Cooper, and P.W. Munro. Theory for the development of neuron selectivity: orientation specificity and binocular interaction in visual cortex. *J Neurosci*, 2(1):32-48, 1982.

[51] K. Pozo and Y. Goda. Unraveling mechanisms of homeostatic synaptic plasticity. *Neuron*, 66(3):337-351, 2010.

[52] G. Turrigiano. Homeostatic synaptic plasticity: local and global mechanisms for stabilizing neuronal function. *Cold Spring Harb Perspect Biol*, 4(1):a005736, 2012.

[53] F. Zenke, G. Hennequin, and W. Gerstner. Synaptic plasticity in neural networks needs homeostasis with a fast rate detector. *PLoS Comput Biol*, 9(11):e1003330, 2013.

[54] F. Zenke, W. Gerstner, and S. Ganguli. The temporal paradox of Hebbian learning and homeostatic plasticity. *Curr Opin in Neurobiol*, 43:166-176, 2017.

[55] H. Sompolinsky, A. Crisanti, and H.J. Sommers. Chaos in Random Neural Networks. *Phys Rev Lett*, 61(3):259-262, 1988.

[56] C. van Vreeswijk and H. Sompolinsky. Chaos in Neural Networks with Balanced Excitatory and Inhibitory Activity. *Science*, 274(5293):1724-1726, 1996.

[57] J. Ginibre. Statistical ensembles of complex, quaternion, and real matrices. *J Math Phys*, 6(3):440-449, 1965.

[58] P.J. Forrester and T. Nagao. Eigenvalue statistics of the real Ginibre ensemble. *Phys Rev Lett*, 99(5):050603, 2007.

[59] J.T. Chalker and B. Mehlig. Eigenvector statistics in non-Hermitian random matrix ensembles. *Phys Rev Lett*, 81(16):3367-3370, 1998.

[60] B. Mehlig and J.T. Chalker. Statistical properties of eigenvectors in non-Hermitian Gaussian random matrix ensembles. *J Math Phys*, 41(5):3233-3256, 2000.

[61] G. Akemann and E. Kanzieper. Integrable structure of Ginibres ensemble of real random matrices and a Pfaffian integration theorem. *J Stat Phys*, 129(5):1159-1231, 2007.

[62] J. Grela. What drives transient behavior in complex systems? *Phys Rev E*, 96(2):022316, 2017.

[63] S. Ganguli, D. Huh, and H. Sompolinsky. Memory traces in dynamical systems. *Proc Natl Acad Sci (USA)*, 105(48):18970-18975, 2008.

[64] K. Rajan and L.F. Abbott. Eigenvalue spectra of random matrices for neural networks. *Phys Rev Lett*, 97(18):188104, 2006.

[65] W. Tarnowski, I. Neri, and P. Vivo. Universal transient behavior in large dynamical systems on networks. *Phys Rev Res*, 2(2):023333, 2020.

[66] G. Biroli, G. Bunin, and C. Cammarota. Marginally
stable equilibria in critical ecosystems. New J Phys, 20(8):083051, 2018.

[67] O. Harish and D. Hansel. Asynchronous rate chaos in spiking neuronal circuits. PLoS Comput Biol, 11(7):e1004266, 2015.

[68] T. Can, K. Krishnamurthy, and D.J. Schwab. Gating creates slow modes and controls phase-space complexity in GRUs and LSTMs. In Mathematical and Scientific Machine Learning, pages 476–511. PMLR, 2020.

[69] K. Krishnamurthy, T. Can, and D.J. Schwab. Theory of gating in recurrent neural networks. Phys Rev X, 12(1):011011, 2022.

[70] F.G. Tricomi. Integral equations. Interscience Publishers, Inc., London & New York, 1957.

## Appendix A: Spectral distribution for the Gaussian Orthogonal Ensemble with eigenvalue constraints

For completeness, we sketch here the derivation of the spectral density for random matrices drawn from the Gaussian Orthogonal Ensemble [35, 38, 39], subject to constraints on the eigenvalues. Excellent pedagogical discussions can be found in Refs [33, 34]. These methods allow us to derive the spectral densities in all cases that we consider in the main text.

Let $M$ be a matrix with size $N \times N$. Assume $M$ is real and symmetric ($M_{ij} = M_{ji}$), and that the individual elements of the $M$ matrix are independent gaussian random numbers,

$$ M_{ii} \sim \mathcal{N}(0, c^2/N), \\ M_{ij} | i \neq j \sim \mathcal{N}(0, c^2/2N). $$

Here, $c$ is the interaction strength. This is the Gaussian Orthogonal Ensemble (GOE). Together with its complex and quaternion counterparts, the Gaussian Ensembles

are the only random matrix ensembles that both have independent entries and are invariant under orthogonal (unitary, symplectic) transformations, which is more obvious when we write the probability distribution of $M$ in terms of its trace:

$$ P_{\text{GOE}}(M) \propto \exp \left( -\frac{N}{2c^2} \text{Tr} M^\top M \right). \quad (\text{A1}) $$

More generally, we consider the ensemble of random matrices with additional constraints on the eigenvalues, i.e. where the probably of a given matrix $M$ can be written as

$$ P(M) \propto \exp \left( -N \sum_{i} f(\lambda_{i}) - \frac{N}{2c^2} \text{Tr} M^{\top} M \right). \quad (\text{A2}) $$

Symmetric matrices can be diagonalized by orthogonal transformations,

$$ M = O^{\top} \Lambda O, $$

where the matrix $O$ is constructed out of the eigenvectors of $M$ and the matrix $\Lambda$ is diagonal with elements given by the eigenvalues $\{\lambda_i\}$. Because $P(M)$ is invariant to orthogonal transformations of $M$, it is natural to integrate over these transformations and obtain the joint distribution of eigenvalues. To do this we need the Jacobian, which is given by the Vandermonde determinant,

$$ dM = \prod_{i \le j} dM_{ij} = \prod_{i < j} |\lambda_i - \lambda_j| d\mu(O) \prod_{i=1}^{N} d\lambda_i, $$

where $d\mu(O)$ is the Haar measure of the orthogonal group under its own action. Now we can integrate over the matrices $O$, or equivalently over the eigenvectors, to obtain

$$ P(\{\lambda_i\}) \prod_{i=1}^{N} d\lambda_i = \int d\mu(O) \prod_{i<j} |\lambda_i - \lambda_j| P(M) \prod_{i=1}^{N} d\lambda_i, \\ P(\{\lambda_i\}) \propto \exp \left( -N \sum_i f(\lambda_i) - \frac{N}{2c^2} \sum_i \lambda_i^2 + \frac{1}{2} \sum_{j \neq k} \ln |\lambda_j - \lambda_k| \right) \\ \propto \exp \left( -N \sum_i u(\lambda_i) + \frac{1}{2} \sum_{j \neq k} \ln |\lambda_j - \lambda_k| \right), $$

where we have defined

$$ u(\lambda) = f(\lambda) + \frac{1}{2c^2} \lambda^2. $$

Mathematically, the distribution $P(\{\lambda_i\})$ is equivalent to

the Boltzmann distribution of a two-dimensional electron gas confined to one dimension, where the electrons are also experiencing a local potential with strength $u(\lambda_i)$. In addition, each pair of the eigenvalues repel each other with logarithmic strength; this term comes from the Van-
dermonde determinant, and gives all the universal features for Gaussian ensembles.

Now we would like to solve for the spectral density. In mean-field theory, which for these problems becomes exact in the thermodynamic limit $N \to \infty$, we can replace sums over eigenvalues by integrals over the spectral density,

$$ \rho(\lambda) = \frac{1}{N} \sum_{i} \delta(\lambda - \lambda_{i}). $$

Then the eigenvalue distribution can be approximated by²

$$ P(\rho(\lambda)) \propto \exp[-N^2 S[\rho(\lambda)]], $$

where

$$ S[\rho(\lambda)] = \int d\lambda \rho(\lambda) u(\lambda) - \frac{1}{2} \int d\lambda d\lambda' \rho(\lambda) \rho(\lambda') \ln|\lambda - \lambda'|. $$

Because $N$ is large, the probability distribution is dominated by the saddle point, $\rho^*$, such that

$$ \left. \frac{\delta \tilde{S}}{\delta \rho} \right|_{\rho=\rho^*} = 0, $$

where $a$ and $b$ are the edges of the support, and

$$ C = \int_{a}^{b} \rho(\lambda) d\lambda. \qquad (A5) $$

If the distribution has a single region of support, then $C = 1$. If the distribution has more than one region of support, then we need to solve with Tricomi's solution separately for each support, and the normalization changes accordingly. In general, solving the equation reduces to finding the edges of the support.

We now apply the general solution (Eq. (A4)) to the three models discussed in the manuscript to obtain the corresponding spectral distributions.

**Model 1, GOE:** For the Gaussian Orthogonal Ensemble without any additional constraints on the eigenvalues, we have $f(\lambda) = 0$. This means the local potential $u(\lambda) = \lambda^2/2c^2$ is symmetric, and hence the resulting eigenvalue distribution is invariant for $\lambda \to -\lambda$, and hence the edges of the support satisfies $a = -b$. Substituting $g(\lambda) = \lambda/c^2$ into Tricomi's solution, we have the integral

$$ \frac{1}{\pi} \mathrm{Pr} \int_{-b}^{b} d\lambda' \frac{\sqrt{\lambda' - b}\sqrt{b - \lambda'}}{\lambda - \lambda'} \lambda' = \lambda^2 - \frac{b^2}{2}. $$

where

$$ \tilde{S} = S + \kappa \int d\lambda \rho(\lambda) $$

has a term with the Lagrange multiplier $\kappa$ to enforce the normalization of the density. Then, the spectral distribution satisfies

$$ u(\lambda) - \int d\lambda' \rho^*(\lambda') \ln|\lambda - \lambda'| = -\kappa. $$

To eliminate $\kappa$ we can take a derivative with respect to $\lambda$, which gives us

$$ g(\lambda) = \frac{du(\lambda)}{d\lambda} = \mathrm{Pr} \int \frac{d\lambda' \rho(\lambda')}{\lambda - \lambda'}, \quad (A3) $$

with the integral defined by its Cauchy principal value. Two methods are common in solving equations of this form. One is the resolvent method, which we will not discuss in detail; see Ref [33]. The other is the Tricomi solution [70], which states that for smooth enough $g(\lambda)$, the solution of Eq (A3) for the density $\rho(\lambda)$ is

$$ \rho(\lambda) = \frac{1}{\pi \sqrt{\lambda - a\sqrt{b - \lambda}}} \left[ C - \frac{1}{\pi} \mathrm{Pr} \int_{a}^{b} d\lambda' \frac{\sqrt{\lambda' - a}\sqrt{b - \lambda'}}{\lambda - \lambda'} g(\lambda') \right], \quad (A4) $$

We expect the density to fall to zero at the edges of the support, rather than having a jump. Thus, we impose $\rho(a) = \rho(b) = 0$, which sets $b = \sqrt{2c}$, and the spectral density becomes

$$ \rho_{\mathrm{GOE}}(\lambda) = \frac{1}{\pi c} \sqrt{2 - \frac{\lambda^2}{c^2}}. $$

This is Wigner's semicircle law.

**Model 2, GOE with hard threshold:** What is the spectral distribution for random matrices drawn from the GOE, but with all eigenvalues constrained to be less than a given threshold $\lambda_{\mathrm{thres}} = 1$? This problem was first solved by Dean and Majumbda [40–42], by recognizing that the eigenvalues remain to satisfy Eq. (A4) with the same $g(\lambda) = \lambda$ as in the GOE case, while the conditions for the edges of the support become different: namely, the two edges of the support are no longer symmetric with respect to $\lambda = 0$. By setting the right edge of the support $b = 1$ and using $\rho(a) = 0$, one can solve for the spectral distribution and find if $c > c_{\mathrm{crit}} = 1/\sqrt{2}$,

$$ \rho_{\mathrm{hard}}(\lambda) = \frac{1}{c^2} \frac{\sqrt{\lambda + l^* - 1}}{2\pi\sqrt{1-\lambda}} (l^* - 2\lambda), $$
where the width of the support is

$$l^* = \frac{2}{3} \left( 1 + \sqrt{1 + 6c^2} \right).$$

Intuitively, one can think about the eigenvalues with hard threshold using the analogy of the two dimensional electrons confined by a wall, or a semi-well potential.

**Model 3, GOE with soft threshold:** In this model, the soft threshold we impose on the eigenvalues

can be incorporated in the eigenvalue distribution (Eq. (23)) by recognizing

$$u(\lambda) = \frac{\lambda^2}{2c^2} + \frac{\xi}{1-\lambda}, \quad g(\lambda) = \frac{\partial u}{\partial \lambda} = \frac{\lambda}{c^2} + \frac{\xi}{(1-\lambda)^2}.$$

Plugging in Tricomi's solution (Eq. (A4)), and redefining the variable $x \equiv \lambda - a$, where $a$ is the left edge of the support, we obtain the spectral distribution,

$$f(x) \equiv \rho(x+a) = \frac{1}{\pi \sqrt{x(l-x)}} \left[ 1 - \frac{1}{\pi} \text{Pr} \int_{0}^{l} dx' \frac{\sqrt{x'(l-x')}}{x-x'} \left( \frac{x'+a}{c^2} + \frac{\xi}{(1-a-x')^2} \right) \right]. \quad (A6)$$

where $l$ is the width of the support.

The spectral distribution can be solved analytically. The first term of the Cauchy principle value in Eq. (A6) can be evaluated exactly, using

$$\text{Pr} \int_{0}^{l} dx' \frac{\sqrt{x'(l-x')}}{x-x'} x' = -\frac{\pi}{8} (l^2 + 4lx - 8x^2)$$

and

$$\text{Pr} \int_{0}^{l} dx' \frac{\sqrt{x'(l-x')}}{x-x'} a = a\pi(x - \frac{l}{2}).$$

The second term of the integral can only be evaluated

when $l > 1 - a$, which always holds when $N$ is finite. The Cauchy principle value is

$$\text{Pr} \int_{0}^{l} dx' \frac{\sqrt{x'(l-x')}}{x-x'} \frac{\xi}{(1-a-x')^2} = -\frac{\pi\xi}{2} \frac{(l-2x)(1-a) + lx}{\sqrt{(1-a)(1-a-l)(1-a-x)^2}}.$$

Denoting the gap between the right edge of the support and the critical eigenvalue $\lambda_{\text{thres}} = 1$ as $g_0$, and recognizing that $a = 1 - g_0 - l$, we obtain the spectral distribution (Eq. (A6)) as given in the main text:

$$\rho_{\text{soft}}(\lambda) = \frac{1}{\pi \sqrt{(\lambda - 1 + g_0 + l)(1 - g_0 - \lambda)}} \left[ 1 + \frac{l^2}{8c^2} + \left(1 - g_0 - \frac{1}{2}l\right) \frac{\lambda}{c^2} - \frac{\lambda^2}{c^2} + \frac{\xi}{2} \frac{(2g_0 - 2g_0^2 + l - 2g_0l) - (2g_0 + l)\lambda}{\sqrt{g_0(g_0 + l)(\lambda - 1)^2}} \right].$$

The only unknowns are the value of the gap $g_0$ and the width $l$ of the support, which we can solve for by setting $f(0) = f(l) = 0$, and obtaining set of coupled equations

$$\left\{ \begin{array}{l} 1 + \frac{l^2}{8c^2} + \frac{(1 - g_0 - l)l}{2c^2} + \frac{\xi l}{2} (g_0 + l)^{-3/2} g_0^{-1/2} = 0, \\ 1 - \frac{3l^2}{8c^2} - \frac{(1 - g_0 - l)l}{2c^2} - \frac{\xi l}{2} (g_0 + l)^{-1/2} g_0^{-3/2} = 0. \end{array} \right. \quad (A7)$$

This set of coupled equation gives us the scaling behavior of the gap $g_0$ and the correlation time $\tau_{\text{corr}}$ for different interaction strength $c$. Since the support $l$ shows up in many terms in the form of $l^2/c^2$, we can deduce that $l$ expands with the interaction strength $c$, as in both the GOE case and the GOE with hard threshold case. For plotting in the main text, the exact values of $g_0$ and $l$ are

solved numerically.

## Appendix B: Compute autocorrelation function in Model 2

The computation of autocorrelation function in Model 2 is subtle, since

$$C_{\text{hard}}(t) = \int d\lambda \frac{\rho(\lambda)}{(1-\lambda)} e^{-(1-\lambda)|t|} \sim \int dk \, k^{-1/2} k^{-1} e^{-k|t|}$$

is not integrable. After introducing an infrared (IR) cut-off at $\epsilon \sim g_0 \sim N^{-2}$, we can write the autocorrelation function as

$$C_{\text{hard}}(t) \approx \int_{\epsilon} dk k^{-3/2} e^{-kt} = -2t^{1/2} \Gamma\left(\frac{1}{2}, \epsilon t\right) + 2\epsilon^{-1/2} e^{-\epsilon t},$$
where $\Gamma(a, z)$ is the incomplete gamma function

$$ \Gamma(a, z) = \int_{z}^{\infty} u^{a-1} e^{-u} du. $$

The resulting autocorrelation coefficient is

$$ R_{\text{hard}}(t) = \frac{C_{\text{hard}}(t)}{C_{\text{hard}}(0)} \approx e^{-\epsilon t} - (\epsilon t)^{1/2} \Gamma\left(\frac{1}{2}, \epsilon t\right) \\ = 1 - \sqrt{\pi}(\epsilon t)^{1/2} + \epsilon t + \mathcal{O}((\epsilon t)^2), \qquad (\text{B1}) $$

and the correlation time

$$ \tau_{\text{corr}} = \frac{\int_{\epsilon} dk \rho(k) k^{-2}}{\int_{\epsilon} dk \rho(k) k^{-1}} \sim \frac{\int_{\epsilon} dk k^{-5/2}}{\int_{\epsilon} dk k^{-3/2}} \sim \epsilon^{-1} \sim N^2. \quad (\text{B2}) $$

## Appendix C: Scaling argument for finite gap in symmetric systems with soft threshold

In this Appendix we will use scaling argument to show that in Model 3 (GOE with soft threshold constraining the norm activity), there is always a finite gap between the right edge of the support of the spectral density $\rho(\lambda)$ and the stability threshold $\lambda^* = 1$.

Let the interaction matrix $M$ be symmetric. Recall Equation (24), the joint probability distribution for the eigenvalues is

$$ P(\{\lambda_i\}) \propto \exp \left[ -\frac{N}{2c^2} \sum_i \lambda_i^2 - N\xi \sum_i \frac{1}{(1-\lambda_i)^\beta} + \frac{1}{2} \sum_{j \neq k} \ln |\lambda_j - \lambda_k| \right]. \quad (\text{C1}) $$

Here, we use $\beta$ as the exponent for the purpose of generality (soft threshold has $\beta > 0$). For Model 3, $\beta = 1$. After coarse-graining and saddle point approximation, we obtain that the spectral density $\rho(\lambda)$ satisfies

$$ \frac{\lambda}{c^2} + \frac{\beta\xi}{(1-\lambda)^{\beta+1}} = \text{Pr} \int \frac{d\lambda' \rho(\lambda')}{\lambda - \lambda'}. \quad (\text{C2}) $$

Let us denote $x \equiv 1 - \lambda$. If there is no gap between largest eigenvalue and $\lambda^* = 1$, we can write the spectral density near $\lambda^*$ as $\tilde{\rho}(x) \equiv \rho(1 - x) = x^\alpha$. Our maximum entropy constraint states that

$$ \int d\lambda \rho(\lambda) \frac{1}{(1-\lambda)^\beta} = \text{const.} $$

$$ \int_0 dx x^{\alpha - \beta} = \text{const.}, $$

which implies $\alpha > \beta - 1$.

For Equation (C2) to make sense, we need

$$ \frac{1-x}{c^2} + \frac{\beta\xi}{x^{\beta+1}} = \text{Pr} \int_0 \frac{dx' \tilde{\rho}(x')}{x' - x} = \text{Pr} \int_0 \frac{dx' x'^{\alpha}}{x' - x} = x^{\alpha}. $$

As $x \to 0$, this equation requires $\beta + 1 = -\alpha$, for both sides to have the same speed of divergence. This requires $\beta < 0$, which is contradictory to the soft threshold assumption.

This contradiction means there has to be a non-vanishing gap $g_0$ between the largest eigenvalue and the stability threshold $\lambda^*$. We can re-write the density function using $y = 1 - g_0 - \lambda$, then Eq. (C2) becomes

$$ \frac{1 - g_0 - y}{c^2} + \frac{\beta\xi}{(y + g_0)^{\beta+1}} = \text{Pr} \int_0 \frac{dy' y'^{\alpha}}{y' - y}. \quad (\text{C3}) $$

The right hand side does not explicitly depend on $g_0$, so we have

$$ g_0 \sim \xi^{1/(\beta+2)}. \quad (\text{C4}) $$

For Model 3, since $\beta = 1$, the gap $g_0 \sim \xi^{1/3}$.

## Appendix D: Derivation for the scaling of time constants in Model 3

For the ensemble of connection matrices with a maximum entropy constraint on the norm activity (Model 3), we are interested in how the time constants, $\tau_{\text{max}}$ and $\tau_{\text{corr}}$, depend on the parameters of the system: the interaction strength $c$, and the constrained norm activity $\mu$. In particular, how the time constants scale with the norm activity when the system is set in different phases by the interaction strength $c$. Here, we analyze the spectral distribution $\rho(\lambda)$, and investigate how the gap size $g_0$, the length of the support $l$, and the constrained norm activity $\mu$ scales with the Lagrange multiplier that we used to constrain the norm, $\xi$. The results are summarized in Table I, and can be visualized in Fig. 3.

We have shown that the spectral distribution is (also see Eq. (25))
$$ 
\begin{aligned}
\rho(\lambda) &= \frac{1}{\pi\sqrt{(\lambda - 1 + g_0 + l)(1 - g_0 - \lambda)}} B(\lambda), \\
B(\lambda) &= \left[ 1 + \frac{l^2}{8c^2} + \left(1 - g_0 - \frac{1}{2}l\right) \frac{\lambda}{c^2} - \frac{\lambda^2}{c^2} + \frac{\xi(2g_0 - 2g_0^2 + l - 2g_0l) - (2g_0 + l)\lambda}{\sqrt{g_0(g_0 + l)(\lambda - 1)^2}} \right].
\end{aligned}
\quad (D1) $$

By setting the spectral density at the edge of the support to zero, we can solve for the scaling dependence of gap size $g_0$ and the length of the support $l$ on the Lagrange multiplier $\xi$. Mathematically, we have

$$ 
\begin{aligned}
B(1 - g_0 - l) &= 0, \\
B(1 - g_0) &= 0.
\end{aligned}
\quad (D2) $$

After simple algebraic manipulation, this set of constraints becomes

$$ 
\begin{aligned}
(2g_0 + l)(8c^2 - 3l^2) + 4l^2 &= 0, & (D3) \\
8 - \frac{l^2}{c^2} - 2\xi l^2 (l + g_0)^{-3/2} g_0^{-3/2} &= 0, & (D4)
\end{aligned}
$$

which sets the scaling relation between the longest time scale $\tau_{max} = 1/g_0$ and the Lagrange multiplier $\xi$.

For the correlation time $\tau_{corr}$, we can express it as the ratio

$$ \tau_{corr} = \frac{\nu}{\mu}, \quad (D5) $$

where the denominator is the expectation value for averaged norm activity,

$$ 
\begin{aligned}
\mu = \langle x_i^2 \rangle &= \int d\lambda \rho(\lambda) \frac{1}{1-\lambda} & (D6) \\
&= \frac{1}{c^2} + \frac{1}{c^2 \sqrt{g(g+l)}} \left( c^2 - g - \frac{l}{2} + \frac{l^2}{8} \right) - \frac{\xi}{8} \frac{l^2}{g^2(g+l)^2},
\end{aligned}
$$

and the numerator can be written as

$$ 
\begin{aligned}
\nu &\equiv \int d\lambda \rho(\lambda) \frac{1}{(1-\lambda)^2} & (D7) \\
&= -\frac{1}{c^2} + \frac{1}{c^2} g^{-3/2} (g+l)^{-3/2} \left( c^2 g + g^3 + \frac{1}{2} c^2 l + \frac{3}{2} g^2 l - \frac{1}{4} l^2 + \frac{5}{8} g l^2 + \frac{1}{16} l^3 \right) - \frac{\xi}{8} \frac{l^2(2g+l)}{g^3(g+l)^3}.
\end{aligned}
$$

Because in the limit of $\xi \gg 1$, $g_0$ is large and there is no long time scales, we focus our discussions on cases where $\xi \ll 1$.

**Case 1:** $c < 1/\sqrt{2}$ - As $\xi$ approaches 0, we recover the semicircle spectral density with the wall far away from the spectrum. This suggests we can write $g_0 = 1 - \sqrt{2c} + A_{-\xi}^{\gamma}$ and $l = 2\sqrt{2c} - B_{-\xi}^{\gamma}$. The expected value for the norm is

$$ 
\begin{aligned}
\lim_{\xi \to 0} \langle x_i^2 \rangle |_{c < 1/\sqrt{2}} &= \int_{-\sqrt{2c}}^{\sqrt{2c}} d\lambda \frac{1}{c\pi} \frac{\sqrt{2 - \lambda^2/c^2}}{1 - \lambda} & (D8) \\
&= \frac{1}{c^2} - \frac{\sqrt{1 - 2c^2}}{c^2}.
\end{aligned}
$$

Taylor expansion leads to $\gamma = 1$. As $\xi$ decreases, both $g_0$ and $\mu$ reaches an upper limit. There is no long time scales in this case.

**Case 2** $c = 1/\sqrt{2}$ - This is the critical case in the hard wall limit, when the gap $g_0$ is 0, and the spectral density remains a semicircle. To solve for the scaling behavior when $\xi \ll 1$, we follow similar procedure as in the previous section and assume $g_0 = A_c \xi^{\gamma}$, $l = l_0 + B_c \xi^{\gamma}$. And again, $l_0 = 2\sqrt{2c} = 2$. Now, the 0th order terms of Eq. (D4) has no terms that scales with $\xi$, requiring us to go to higher orders to solve for $\gamma$, which gives

$$ \gamma = \frac{2}{5}. \quad (D9) $$

The norm is

$$ \langle x_i^2 \rangle = 2 - D\xi^{1/5} + O(\xi^{3/5}). \quad (D10) $$

This is interesting: as $\xi$ decreases, the norm activity $\mu$ approaches the limit $\mu = 2$, and the corresponding $g_0$
<table><caption>TABLE I: Scaling of inverse slowest time scale (gap) $g_0$, width of the support of spectral density $l$, and averaged norm per neuron $\langle x_i^2 \rangle$ versus the Lagrange multiplier $\xi$ (to leading order) in different regimes.</caption><thead><tr><th></th><th></th><th>g<sub>0</sub></th><th>l</th><th>&lang;x<sub>i</sub><sup>2</sup>&rang;</th></tr></thead><tbody><tr><td>&xi; &ll; 1</td><td>0 &lt; c &lt; 1/&radic;2</td><td>1 &minus; &radic;2c + A<sub>&minus;</sub>&xi;</td><td>2&radic;2c &minus; B<sub>&minus;</sub>&xi;</td><td><sup>1</sup>&frasl;<sub>c<sup>2</sub></sub> &minus; <sup>&radic;(1 &minus; 2c<sup>2</sup>)</sup>&frasl;<sub>c<sup>2</sup></sub> &minus; D<sub>&minus;</sub>&xi;</td></tr><tr><td>&xi; &ll; 1</td><td>c = 1/&radic;2</td><td>A<sub>c</sub>&xi;<sup>2/5</sup></td><td>2 &minus; B<sub>c</sub>&xi;<sup>2/5</sup></td><td>2 &minus; D<sub>c</sub>&xi;<sup>1/5</sup></td></tr><tr><td>&xi; &ll; 1</td><td>c &gt; &radic;2</td><td>A<sub>+</sub>&xi;<sup>2/3</sup></td><td>l<sub>0</sub> &minus; B<sub>+</sub>&xi;<sup>2/3</sup></td><td>D<sub>+</sub>&xi;<sup>&minus;1/3</sup> + <sup>1</sup>&frasl;<sub>c<sup>2</sup></sub></td></tr></tbody></table>

continuous to decrease. The system can reach an infinitely long timescale with a bounded dynamic range for individual neurons.

**Case 3** $c > 1/\sqrt{2}$ - In this regime, the spectral density at $\xi = 0$ is divergent near $\lambda = 1$. But for any $\xi > 0$, there is a finite gap between the wall and right edge of the spectrum, so a limit for the spectral density when $\xi \to 0$ is not well-defined.

We can assume the gap $g_0$ takes the form $g_0 = A_+\xi^\gamma$. Because $l$ is of order 1, we can assume $l = l_0 + B_+\xi^\gamma$. Plugging these to Eq. (D3), we can solve the 0th order equation and get

$$l_0 = \frac{2}{3} \left( 1 + \sqrt{1 + 6c^2} \right). \qquad (\text{D11})$$

This $l_0$ is equal to the length of the spectrum at the hard wall limit, i.e. when $\xi = 0$. Plugging $g_0$ and $l$ into Eq. (D4), we get from 0-th order solution

$$\mu = \frac{2}{3} \qquad (\text{D12})$$

and

$$A_+ = 2^{-2/3} l_0^{1/3} \left(2 - \frac{l_0^2}{4c^2}\right)^{-2/3}. \qquad (\text{D13})$$

We notice that the prefactor $A$ is not a monotonic function of $c$. Rather, $A$ takes a minimum $A_{\min} = 1$ at $c = 2$. This is interesting, as we don't want the interaction strength to be too large.

The norm becomes

$$\mu = \int d\lambda \rho(\lambda) \frac{1}{1 - \lambda} = D_+ \xi^{-1/3} + \frac{1}{c^2} + \mathcal{O}(\xi)^{1/3}, \quad (\text{D14})$$

where

$$D_+ = \frac{1}{c^2} A^{-1/2} l_0^{-1/2} \left( c^2 + \frac{l^2}{8} - \frac{l}{2} \right) - \frac{A^{-2}}{8}. \qquad (\text{D15})$$

Similarly, the prefactor $D_+$ has a maximum $D_{\max} = 3/8$ at $c = 2$.

For the correlation time, the denominator is the averaged norm, and the numerator is

$$\begin{aligned} \nu &\equiv \int d\lambda \rho(\lambda) \frac{1}{(1 - \lambda)^2} \\ &= E\xi^{-1} - \frac{1}{c^2} + \mathcal{O}(\xi). \end{aligned} \qquad (\text{D16})$$

The prefactor $E$ has maximum $E_{\max} = 1/8$ at $c = 2$. Then, the correlation time is

$$\tau_{\text{corr}}^{\text{sw}} = \frac{E}{D} \xi^{-2/3} + \mathcal{O}(\xi). \qquad (\text{D17})$$

Interestingly, when we compute the ratio of the two time scales, $\tau_{\max}$ and $\tau_{\text{corr}}$, we found that the ratio is a constant, as

$$\frac{\tau_{\text{corr}}^{\text{sw}}}{\tau_{\max}^{\text{sw}}} = \frac{1}{3} - \frac{AE}{c^2 D^2} \xi^{1/3} + \mathcal{O}(\xi). \qquad (\text{D18})$$

In this case, as $\xi$ approaches 0, the gap $g_0$ decreases, but the norm also increases to infinite. If we only look at the relation between the parametrized variable $g_0$ and $\mu$, we find that $g_0 \sim \mu^{-2}$, which matches what we get from numerically solving the exact Eqs. (D3) and (D4). To get an infinitely long time scale, an unlimited growth of dynamic range for individual neuron is required, which is impossible to realize in biological systems.

## Appendix E: Asymmetric interaction

Recall that the system we consider is a linear dynamical system

$$\frac{d\dot{x}_i}{dt} = -x_i + \sum_j M_{ij}x_j + \eta_i(t), \qquad (\text{E1})$$

with the entries of the interaction matrix $M$ real, and white noise $\eta$ where

$$\langle \eta_i(t) \rangle = 0, \qquad \langle \eta_i(t) \eta_j(t') \rangle = 2\delta_{ij}\delta(t - t').$$

If $M$ can be asymmetric, the solution of this dynamical system gives the steady state autocorrelation at equal time

$$\mu = \frac{1}{N} \langle \sum_i x_i^2 \rangle = \frac{1}{N} \sum_{i,j} O_{ij} \frac{2}{2 - \lambda_i - \lambda_j^\dagger}. \qquad (\text{E2})$$

Here, $O_{ij}$ is the overlap between the left and right eigenvectors, defined as

$$O_{ij} = \langle l_i | l_j \rangle \langle r_j | r_i \rangle, \qquad (\text{E3})$$

where $\langle l_i |$ is the left eigenvector (and $|r_i\rangle$ is the right eigenvector) corresponding to the $i$-th eigenvalue. The
normalization is chosen such that the left and right eigen-vectors form a bi-orthogonal basis, i.e. $\langle l_i|r_j \rangle = \delta_{ij}$.

In Model 3, we consider systems which is stable by constraining the autocorrelation when the system is in steady state. The corresponding maximum entropy distribution for the matrix element is

$$P(M) \propto \exp \left[ -\frac{N}{2c^2} \mathrm{Tr} M^\mathrm{T} M - N\xi \sum_{i,j} O_{ij} \frac{2}{2 - \lambda_i - \lambda_j^\dagger} \right], \quad (\mathrm{E4})$$

<figure><img src="image_8.png" /><figcaption>FIG. 6: Relation between the gap $g_0$ and the Lagrange multiplier $\xi$ for real asymmetric matrix $M$ subject to constraints on the norm. The connection strength $c = 2$. Slow modes are in general more likely to be found in eigenvalues which form complex conjugate pairs (right), than real eigenvalues (left), suggesting slow modes are more likely oscillatory than exponential decay.</figcaption></figure>

This probability measure of the matrix can no longer be separate into the product measure of eigenvalue and eigenvectors due to the term with eigenvector overlaps. Another source of complexity comes from additional constraints for eigenvalues of real asymmetric matrices, where all eigenvalues have to be real or form complex conjugate pairs.

To get an idea about whether the gap in the symmetric case remains in asymmetric systems, we directly sample matrices from the distribution [Eq. (E4)] using Monte Carlo Markov Chains. As shown in Fig. 6, there is still a finite gap $g_0$ between the right edge of the spectrum and the stability threshold, the same as in the symmetric case. It is interesting to notice that the real eigenvalues and the complex eigenvalues behave differently: the gap between the complex eigenvalue and the wall can be fitted to a $g_0 \sim \xi^{1/3}$ curve, while the gap between the largest real eigenvalue and the wall seems to be independent of $\xi$, suggesting that slower modes are more likely to be found oscillatory than exponential decay. Other interesting thing to notice is that this 1/3 scaling is slower than what we got from symmetric systems.