Dynamically selected steady states and criticality in non-reciprocal networks

doi:10.1016/j.chaos.2024.114809

Chaos, Solitons & Fractals

Volume 182, May 2024, 114809

https://doi.org/10.1016/j.chaos.2024.114809 Get rights and content

Highlights

•
Two ordered phases emerge in unbalanced neural networks for arbitrary link asymmetry.
•
In spin-glass (SG) phases dynamics is chaotic; in ferromagnetic (F) phases it reaches a fixed point.
•
As link reciprocity increases, chaotic effects weaken, until marginal stability is eventually attained.
•
Both ordered phases shrink as reciprocity decreases, with SG eventually disappearing.
•
Intriguing analogies between chaotic and SG phases are put forward.

Abstract

We consider a simple neural network model, evolving via non-linear coupled stochastic differential equations, where neural couplings are random Gaussian variables with non-zero mean and arbitrary degree of reciprocity. Using a path-integral approach, we analyse the dynamics, averaged over the network ensemble, in the thermodynamic limit. Our results show that for any degree of reciprocity in the couplings, two types of criticality emerge, corresponding to ferromagnetic and spin-glass order, respectively. The critical lines separating the disordered from the ordered phases is consistent with spectral properties of the coupling matrix, as derived from random matrix theory. As the non-reciprocity (or asymmetry) in the couplings increases, both ordered phases diminish in size, ultimately resulting in the disappearance of the spin-glass phase when the couplings become anti-symmetric. We investigate non-fixed point steady-state solutions for uncorrelated interactions. For such solutions the time-lagged correlation function evolves according to a gradient-descent dynamics on a potential, which depends on the stationary variance. Our analysis shows that in the spin-glass region, the variance dynamically selected by the system leads the correlation function to evolve on the separatrix curve, limiting different realizable steady states, whereas in the ferromagnetic region, a fixed point solution is selected as the only realizable steady state. In the spin-glass region, stationary solutions are unstable against perturbations that break time-translation invariance, indicating chaotic behaviour in large single network instances. Numerical analysis of Lyapunov exponents confirms that chaotic behaviour emerges throughout the spin-glass region, for any value of the coupling correlations. While negative correlations increase the strength of chaos, positive ones reduce it, with chaos disappearing for reciprocal (i.e. symmetric) couplings, where marginal stability is attained. On the other hand, in finite size non-reciprocal networks, fixed points and limit cycles can arise in the spin-glass region, especially close to the critical line. Finally, we show that chaos is suppressed when the strength of external noise exceeds a certain threshold. Intriguing analogies between chaotic phases in non-equilibrium systems and spin-glass phases in equilibrium are put forward.

Significance statement: Diverse equilibrium systems with heterogeneous interactions lie at the edge of stability. Such marginally stable states are dynamically selected as the most abundant ones or as those with the largest basins of attraction. On the other hand, systems with non-reciprocal (or asymmetric) interactions are inherently out of equilibrium, and exhibit a rich variety of steady states, including fixed points, limit cycles and chaotic trajectories. How are steady states dynamically selected away from equilibrium? We address this question in a simple neural network model, with a tunable level of non-reciprocity. Our study reveals different types of ordered phases and it shows how non-equilibrium steady states are selected in each phase. In the spin-glass region, the system exhibits marginally stable behaviour for reciprocal (or symmetric) interactions and it smoothly transitions to chaotic dynamics, as the non-reciprocity (or asymmetry) in the couplings increases. Such region, on the other hand, shrinks and eventually disappears when couplings become anti-symmetric. Our results are relevant to advance the knowledge of disordered systems beyond the paradigm of reciprocal couplings, and to develop an interface between statistical physics of equilibrium spin-glasses and dynamical systems theory.

1. Introduction

Countless systems, both natural and human-made — such as species in ecosystems, neurons in neural networks, agents in financial markets and atomic magnetic moments in disordered materials — can be visualized as heterogeneous networks of interacting elements. In the case of fully connected systems with symmetric or reciprocal interactions, the general phenomenology is that upon increasing the heterogeneity in pairwise interactions, there is a transition from a phase with a single equilibrium to a “spin-glass” phase with multiple equilibria [1], [2], [3], [4]. The latter is characterized by ergodicity breaking, slow dynamics, memory, and the breakdown of time-translation invariance (TTI) among other non-trivial features. The spin-glass (SG) phase exhibits a complex energy landscape with many minima, separated by large barriers, hence different copies or “replicas” of the system typically relax in distinct regions of the phase space [1]. This leads to the phenomenon of replica-symmetry breaking (RSB), namely the breakdown of symmetry among different replicas, when averaging over the distribution of interactions [1]. It is well-known that there are two main classes of RSB, the one-step RSB (1RSB) and the full RSB (fRSB), and that these are characterized by a different structure of the metastable states: marginal (and unstable under external perturbations) in fRSB and well-shaped and surrounded by large barriers in 1RSB [5], [6], [7].

Many systems in different situations (physics, biology, ecology, neuroscience and economics) are intriguingly observed to lie just at the edge of stability or operate near criticality [8], [9], [10], [11], [12], [13], [14], [15], [16], [17]. The recent application of RSB frameworks to some of such systems has been highly revealing, suggesting that they may dynamically select marginally stable states, as they could be the most numerous and/or those with the largest basins of attraction, underpinning a fRSB scenario [18], [19], [20], [21], [22], [23]. Such a powerful framework, however, is limited to equilibrium systems with reciprocal (or symmetric) interactions, which preserve detailed balance.

For systems with non-reciprocal interactions, a general framework to quantify the rich diversity of emerging steady states, typically including limit cycles and chaotic trajectories, in addition to fixed points, is yet to be established — although recent years have seen some progress in this direction [24], [25], [26], [27], [28] — and there is, as yet, no method to predict which steady state is dynamically selected by the system.

It is known that systems with asymmetric interactions can exhibit a transition to chaotic behaviour as the randomness in their interactions is increased [29], [30] as well as exotic so-called non-reciprocal phase transitions [31], that constitute a topic of current research interest.

In addition, it has been shown that, for asymmetric neural networks close to the transition into a chaotic phase, the mean number of fixed points diverges exponentially with network size, with the rate of divergence, called topological complexity, following the trend of the largest Lyapunov exponent [32]. However, how topological complexity affects dynamical complexity remains unclear; indeed, a formal connection between spin-glass-like phases with exponentially many equilibria and chaotic dynamics has not been established yet.

The aim of this manuscript is to advance the development of an interface between the statistical physics of classical disordered systems and dynamical systems theory, by showing that chaotic phases in non-equilibrium systems exhibit certain analogies with spin glass phases in equilibrium (fRSB) systems.

We revisit a classical model for neural networks dynamics — introduced by Sompolinsky, Crisanti and Sommers in [29] — and modify it to allow for interactions to be unbalanced (i.e. with a non-vanishing mean) and randomly asymmetric, as well as to incorporate “thermal noise” in the dynamics.

We bypass the (hard) problem of calculating the topological complexity of different attractors in the configuration space and focus on dynamical aspects (e.g. time-dependent moments) of the trajectories at stationarity, as originally proposed by SCS in [29].

By computing two order parameters, namely the mean and variance of neural activity, we reveal the emergence of two types of criticality, one corresponding to ferromagnetic order and the other corresponding to spin-glass order. The transition from the disordered to the ordered phase can be rationalized in terms of spectral properties of the interactions matrix.

Numerical analysis of the largest Lyapunov exponents (LLE) shows that chaotic dynamics arises in the spin-glass region of the phase diagram, for any value of the asymmetry of the couplings.

For uncorrelated couplings, we analyse the non-fixed point steady states and we determine the value of the variance that is dynamically selected by the system. We find that whenever the parameters allow for the existence of a multitude of physical (i.e. realizable) steady-state solutions, the system selects the steady state corresponding to the separatrix curve, which delimitates the basins of attractions of different steady states. This is reminiscent of the behaviour of fRSB systems at equilibrium, which select marginally stable states, i.e. saddles in the free-energy landscape, that are linked together by flat directions. Our analysis suggests that chaotic motion is the manifestation, at the level of single network instances, of such an ensemble averaged dynamics lying on the separatrix curve, delimiting different realizable steady states. Such a motion is not time-translation invariant, similarly to ageing dynamics in spin-glass phases. Finally, we show that upon increasing the signal (i.e. the imbalance of interactions), or the thermal noise, the system is brought to a phase where only one bounded solution exists and chaos is suppressed. Such a behaviour can be seen again, as mirroring equilibrium phase transitions from spin-glass phases (with many equilibria) to ferromagnetic or paramagnetic phases (with a single equilibrium), when the signal or the noise, respectively, are raised above a critical value.

The manuscript is organized as follows: In Section 2, we define the model. Section 3 delves into the linear stability of the quiescent state, revealing the existence of two types of criticality. Section 4 explores the dynamics of the model via a path integral formalism. In particular, we first consider the dynamics in the absence of noise for uncorrelated interactions (Section 4.1). Then, we shift our focus to the noisy dynamics for uncorrelated interactions (Section 4.2). Finally, we analyse the case of correlated interactions (Section 4.3). We discuss the main results and conclusions in Section 5. For further technical details and complementary analyses, readers can refer to the Appendices and the Supplementary Material (SM).

2. Rate model

We consider a variant of the classical model introduced by Sompolinsky, Crisanti and Sommers (SCS) in [29]. The model consists of a fully-connected network with

N

neurons (as shown Fig. 1, left panel), each one described by its time-dependent firing rate,

x_{i} (t)

, with

i = 1, 2, \dots, N

. Their dynamics obey the following set of coupled stochastic differential equations (1)

{\dot{x}}_{i} = - x_{i} + f (g \sum_{j} J_{i j} x_{j} + θ_{i}) + ξ_{i}, i = 1, \dots, N

where the first term on the right hand side describes a spontaneous decay of activity; the second one is the response to inputs, mediated by the gain function

f (x)

which is a monotonically increasing function which saturates for large (in absolute value) arguments, taken here, without loss of generality, to be

f (x) = tanh (x)

;

g

is an overall coupling strength and the “synaptic weights”,

J_{i j}

, are (quenched) Gaussian random variables with mean

J_{0} / N

and variance

J^{2} / N

(2)

\bar{J_{i j}} = \frac{J_{0}}{N}; \bar{J_{i j}^{2}} - \frac{J_{0}^{2}}{N^{2}} = \frac{J^{2}}{N}

where overbars stand for network average. We allow for correlations between

J_{i j}

and its reciprocal,

J_{j i}

, (3)

\bar{J_{i j} J_{j i}} - \frac{J_{0}^{2}}{N^{2}} = γ \frac{J^{2}}{N}

quantified by the Pearson’s correlation coefficient

γ \in [- 1, 1]

, so that for

γ = 1

, interactions are symmetric/reciprocal i.e.

J_{i j} = J_{j i}

, for

γ = 0

they are uncorrelated (or so-called “fully asymmetric”), and for

γ = - 1

they are anti-symmetric. The external field

θ_{i}

is added to generate response functions and will be set to zero afterwards. Finally, the third term,

ξ_{i} (t)

, is a zero-mean Gaussian white noise with

〈 ξ_{i} (t) ξ_{j} (t^{'}) 〉 = 2 σ^{2} δ_{i j} δ (t - t^{'})

, so that

σ

is the noise amplitude, and

〈 . . . 〉

stands for noise average.

Let us remark that this model is very similar to the classical one of SCS but differs in a few aspects: (i) the non-linear function applies to the sum of inputs rather than to each of them; this brings the model closer to rate models in neuroscience [33]; (ii) the connectivity matrix is allowed here to have a non-vanishing mean, so that one can study “excitation” (

J_{0} > 0

) and “inhibition” (

J_{0} < 0

) dominated regimes (we refer to these as “unbalanced” cases); (iii) interactions have an arbitrary degree of symmetry/reciprocity (encoded in

γ

); (iv) noise

σ^{2} \neq 0

is included in the dynamics.

In what follows, we discuss two alternative approaches to tackle analytically different aspects of the model: (i) a linear stability analysis of the fixed points of the noiseless dynamics, allowing us to compute phase boundaries and (ii) a more complete path-integral formalism, allowing for a full (dynamic mean-field) solution (exact in the infinite network-size limit).

3. Linear-stability analysis

To scrutinize the possible steady states of the noiseless dynamics we start by considering a linear stability analysis of the “quiescent” solution, where

x_{i} (t) = 0, \forall i = 1, \dots, N

and the neural network remains inactive. For this we fix

θ_{i} = 0

σ_{i} = 0

in Eq. (1), and expand

f (x)

to the first order

{\dot{x}}_{i} = - x_{i} + g \sum_{j} J_{i j} x_{j}, i = 1, \dots, N .

Thus, the stability of the quiescent solution is controlled by the largest eigenvalue,

λ_{M}

, of the connectivity matrix

J

, which depends on the model parameters

(J_{0}, J)

as well as the degree of correlation,

γ

. In particular, the quiescent solution is stable if

g λ_{M} < 1

, and becomes unstable for values of

g

above

g_{c} = 1 / λ_{M}

The maximum eigenvalue can be easily computed (in the infinite network-size limit) using well-known results from random matrix theory (RMT) for the case of fully asymmetric (

γ = 0

) and symmetric (

γ = 1

) networks. In particular, for fully asymmetric networks with

| J_{0} | < J

the eigenvalues obey the Girko’s circle law, meaning they lay uniformly inside a circle of (spectral) radius

J

in the complex plane centred at the origin [34], [35], [36], i.e. the maximum eigenvalue coincides with

J

. If, however,

| J_{0} | > J

, then one of the eigenvalues leaves the circle becoming a real-valued “outlier” with value

J_{0}

[36], [37]. On the other hand, for symmetric networks, with

| J_{0} | < J

, the eigenvalues are all real and obey the Wigner’s semicircle law, so that they are uniformly distributed on a semicircle of radius

2 J

[38], i.e. the maximum eigenvalue coincides with

2 J

. For symmetric networks with

| J_{0} | > J

, the distribution of the bulk keeps following Wigner’s law, but a single outlier eigenvalue appears at

J_{0} + J^{2} / J_{0}

[39].

Hence, the limit of stability,

g_{c}

, is given by the largest element between the outlier and the largest eigenvalue in the bulk, i.e. one has, for

γ = 0

(4)

g_{c} max (J, J_{0}) = 1 \to \frac{1}{g_{c} J} = max (1, \frac{J_{0}}{J})

while, for

γ = 1

(5)

g_{c} max (2 J, J_{0} + \frac{J^{2}}{J_{0}}) = 1 \to \frac{1}{g_{c} J} = max (2, \frac{J_{0}}{J} + \frac{J}{J_{0}}) .

Furthermore, for general

γ

’s, the distribution of eigenvalues has been derived for balanced matrices in [40], and the position of the outlier has been calculated for unbalanced matrices in [41].

Thus, depending on whether the largest eigenvalue is an outlier or is at the edge of the bulk, two different paths to instability emerge. If the outlier is the largest eigenvalue, as

g

is increased, a single collective mode (given by the eigenvector associated to the outlier eigenvalue) destabilizes, as illustrated in the central panel of Fig. 1. In this case, the system reaches a fixed steady state above the instability threshold. If, on the other hand, the maximum eigenvalue lies at the edge of the bulk, a continuum of modes may become destabilized as

g

is increased across the edge of instability, as illustrated in Fig. 1, right panel. In this case, the dynamics of the system turns out to be non-trivial and complex trajectories emerge. This difference is the core of the distinction between type-I and type-II criticalities [42], [43].

It is important to highlight that linear stability analysis does not enable us to characterize the two types of criticality in terms, for example, of order parameters and critical exponents, nor does it help us to locate the transition between type-I and type-II phases, away from the quiescent state. Thus, to tackle analytically this general problem one needs to resort to a more sophisticated dynamical mean-field approach as developed in the following section.

Before concluding this section, it is worth remarking that a linear stability analysis of the symmetrized matrix (

H = \frac{1}{2} (J + J^{T}) - I

) [44], [45], [46], [47] enables us to discern a regime within the linearly stable phase where the dynamics is “reactive”. This implies that perturbations are amplified before returning to the steady state. Reactive dynamics has been extensively discussed in neuroscience, and significant functional advantages have been attributed to it [48], [49]. In SM Sec. S.IV we show that the model under study exhibits such regimes of reactivity.

4. Path-integral formalism

To analyse the full phase diagram one can resort to a “dynamic mean-field” approach (DMF) as originally introduced by SCS [29] and later developed in a number of works (e.g, [50], [51], [52], [53]). Alternatively, one can employ a more systematic path-integral formalism (or generating functional analysis) as originally introduced by Martin, Siggia and Rose [54] and successfully applied to neural networks in [55], [56], [57], which allows to formally derive the DMF equations and quantify their limits of validity [50], [57]. The calculation is quite standard but, for the sake of completeness, here we reproduce the main steps.

The starting point of the method is the definition of the moment-generating functional of the random (vector) function

x (t)

(6)

Z [ψ] = 〈e^{i \sum_{i} \int d t x_{i} (t) ψ_{i} (t)}〉

where the average

〈 \dots 〉

is taken over the probability distribution

P [x (t)]

of trajectories

x (t)

, generated by Eq. (1), for a fixed (quenched) connectivity matrix

J

. From such a generating functional one can easily compute the chief quantities describing the system’s collective dynamics, such as the mean activity of a neuron

i

(7)

M_{i} (t) = 〈 x_{i} (t) 〉 = {\frac{δ Z [ψ]}{i δ ψ_{i} (t)}|}_{ψ = 0},

the two-time pairwise correlations (8)

C_{i j} (t, t^{'}) = 〈 x_{i} (t) x_{j} (t^{'}) 〉 = {\frac{δ Z [ψ]}{i δ ψ_{i} (t) i δ ψ_{j} (t^{'})}|}_{ψ = 0},

and the response functions (9)

R_{i j} (t, t^{'}) = \frac{δ 〈 x_{i} (t) 〉}{δ θ_{j} (t^{'})} = {\frac{δ^{2} Z [ψ]}{i δ θ_{j} (t^{'}) δ ψ_{i} (t)}|}_{ψ = 0} .

In order to compute disorder-averaged quantities, one needs to average the generating functional

Z [ψ]

over the distribution

P (J)

. This procedure leads — once the infinite network-size (

N \to \infty

) limit has been taken — to a self-consistent stochastic equation for the dynamics of a representative neuron, called dynamical-mean-field equation [29], [50], [58], see Supplemental Material (SM) Sec. S.I for further details: (10)

\dot{x} (t) = - x (t) + tanh [J_{0} g M (t) + γ J^{2} g^{2} \int d t^{'} R (t, t^{'}) x (t^{'}) + θ (t) + ϕ (t)] + ξ (t) .

In this equation,

ξ (t)

is a zero-mean Gaussian white noise with

〈 ξ (t) ξ (t^{'}) 〉 = 2 σ^{2} δ (t - t^{'})

ϕ (t)

is a random Gaussian field with zero mean and auto-correlation (11)

〈 ϕ (t) ϕ (t^{'}) 〉 = J^{2} g^{2} C (t, t^{'})

and the order parameters

M (t)

C (t, t^{'})

and

R (t, t^{'})

need to be calculated self-consistently, as averages over realizations of the effective single-neuron process: (12)

M (t) = 〈 x (t) 〉,

(13)

C (t, t^{'}) = 〈 x (t) x (t^{'}) 〉

(14)

R (t, t^{'}) = \frac{δ 〈 x (t) 〉}{δ θ (t^{'})} = 〈 \frac{δ x (t)}{δ ϕ (t^{'})} 〉,

where in the second equality of the last equation we have used that

θ

and

ϕ

have the same role in Eq. (10). The random function

ϕ (t)

can be interpreted as an “interference” term that quantifies the impact, on the dynamics of the representative neuron, of all the other neurons that interact with it, while the second term in the square brackets of Eq. (10) accounts for the so-called “retarded” self-interactions of the representative neuron with its own past, i.e. past values

x (t^{'})

influence

x (t)

at later times

t > t^{'}

. This is due to the fact that the representative neuron sends signal to its neighbours through links which are correlated with those used by the neighbours to send signal back to it. We note that (as usual) this term only arises when such correlations in the links are present (i.e.

γ \neq 0

) and it makes the effective dynamics non-Markovian. From here on, the external field is fixed to zero,

θ (t) = 0

In the following sections, we investigate separately the cases of (i) uncorrelated interactions (

γ = 0

) with noiseless dynamics (

σ = 0

); (ii) uncorrelated interactions (

γ = 0

) with noisy dynamics (

σ \neq 0

) and (iii) correlated interactions (

γ \neq 0

), focusing for simplicity on the noiseless case.

4.1. Uncorrelated interactions ( $γ = 0$ ) and noiseless dynamics ( $σ = 0$ )

For uncorrelated interactions (i.e.

γ = 0

), the equation of motion for the noiseless dynamics (

σ = 0

) takes a particularly simple form (15)

\dot{x} (t) = - x (t) + tanh [J_{0} g M (t) + ϕ (t)],

that can be explicitly analysed.

4.1.1. Fixed-point solutions

Let us start by assuming that the system reaches a single stable fixed point

x (t) = x

. In such a case,

〈 x 〉 = M

〈 x^{2} 〉 = C (0) = q

, and each realization of

ϕ (t)

becomes a static zero-mean Gaussian random variable with variance

J^{2} g^{2} q

[59], [60]. Therefore, imposing the stationarity condition,

\dot{x} = 0

, in Eq. (15) and averaging over

ϕ

one readily obtains: (16)

M = 〈 tanh (J_{0} g M + ϕ) 〉 = \int D ψ tanh (J_{0} g M + g J \sqrt{q} ψ)

where the short-hand notation

D ψ = e^{- ψ^{2} / 2} d ψ / \sqrt{2 π}

has been introduced. Similarly, the variance

q

at the fixed-point solution is (17)

q = 〈 {tanh}^{2} (J_{0} g M + ϕ) 〉 = \int D ψ {tanh}^{2} (J_{0} g M + g J \sqrt{q} ψ) .

An alternative, simpler derivation of Eqs. (16)–(17) — not generalizable to the

γ \neq 0

case — is presented in SM (Sec. S.II).

Before proceeding, we note that Eqs. (16), (17) are identical to the equations for the magnetization and the Edward–Anderson order parameter, respectively, in the well-known Sherrington–Kirkpatrick (SK) model for spin-glasses [61], within a replica symmetric ansatz, with

1 / g

playing the role of the temperature (see e.g. Eqs.(2.28)-(2.30) in [3]). Such an equivalence, at the level of order parameters, between the SK model and the neural-network rate model under study, is surprising as the two models are defined in rather different ways. In particular, the SK model is defined for Ising (i.e.

\pm 1

) variables with symmetric interactions while here Eqs. (16), (17) have been derived for continuous state variables and non-symmetric interaction matrices.

From the analogy with spin glasses, it is well-known that Eqs. (16), (17) have a trivial disordered solution (referred to as “paramagnetic” in the context of magnetic systems) with

(M, q) = (0, 0)

, which is stable at sufficiently high temperatures (corresponding here to small values of the coupling constant

g

). From the disordered phase, two different types of order may emerge, depending on the relative values of the mean

J_{0}

and variance

J

of the coupling distribution.

In particular, for

J_{0} > J

, there is a bifurcation at

g J_{0} = 1

from the disordered phase to an ordered (or “ferromagnetic”) phase in which the up-down symmetry is broken, i.e., there is a non-vanishing magnetization,

(M \neq 0, q \neq 0)

. On the other hand, for

J_{0} < J

, a bifurcation occurs at

g J = 1

from such a disordered phase to a “spin-glass” phase

(M = 0, q \neq 0)

[61]. Therefore, the disordered phase loses its stability at the critical value

g_{c}

, given by Eq. (4), in agreement with results from linear-stability analysis. In particular, type-I criticality discussed in Section 3 can be identified with the transition from the paramagnetic to the ferromagnetic phase, while type-II criticality corresponds to the transition from the paramagnetic to the spin-glass phase. Such criticalities can therefore be fully characterized by the order parameters

M

and

q

, which correspond to mean neural activity and mean-squared activity, respectively. In addition, the critical line separating the ferromagnetic from the spin-glass phase can be obtained by analysing the limit of linear stability of the

M = 0

solution, which — as derived in the SM (Sec. S.III) — leads to (18)

\frac{1}{g J} = \frac{J_{0}}{J} (1 - q^{*})

where

q^{*}

is the solution of (19)

q^{*} = \int D ψ {tanh}^{2} (g J \sqrt{q^{*}} ψ) .

The resulting phase diagram, in the parameter space (

J_{0} / J, 1 / g J

), is shown in Fig. 2, where the solid white lines in panels (A) and (B) represent the critical lines separating the paramagnetic (P), ferromagnetic (F) and spin-glass (SG) phases.

Expanding the self-consistency Eqs. (16)–(17) to the first sub-leading order, it is possible to obtain the critical exponents

α

and

β

such that

M \sim {(g - g_{c})}^{α}

and

q \sim {(g - g_{c})}^{β}

(see SM, Sec. S.III). For the ferromagnetic transition

J_{0} > J

, one has

α = 1 / 2

and

β = 1

(as in the Ising model) whereas for the SG transition, i.e. for

J_{0} < J

, the exponent

α

is trivially

0

and

β = 1

In addition to fixed-point solutions, a rich set of possible non-fixed-point solutions is expected to appear with non-symmetric interactions, including limit cycles and chaotic trajectories [29], [51], [57]. Among these, stationary solutions are characterized by a time-independent average

M (t) = M

and a time-translation invariant two-time correlator

C (t, t^{'}) = C (t - t^{'})

, which implies a time-independent equal-time correlator

C (t, t) = C (0)

(note that for fixed-point solutions

C (0) = q

, however for non-fixed-point solutions

C (0) \neq q

). Whenever the steady state is not a fixed point, node activities

x_{i} (t)

fluctuate in time, even at stationarity. One can therefore ask whether the phase diagram derived for fixed-points gives reliable information on the phase transition behaviour of the system, whose dynamics is not necessarily at a fixed point.

To address this question, we ran

S

independent computer simulations of the microscopic dynamics and, assuming that the system eventually reaches a steady-state, we computed the stationary magnetization and equal-time (or zero-lag) correlator for each simulated trajectory

x_{s} (t)

, with

s = 1, \dots, S

, as the sample average of the mean activity¹ (20)

{\hat{M}}_{s} = \frac{1}{t_{m}} \sum_{t = t_{0}}^{t_{0} + t_{m}} \frac{1}{N} \sum_{i = 1}^{N} x_{i, s} (t)

and mean-squared activity (21)

{\hat{C}}_{s} (0) = \frac{1}{t_{m}} \sum_{t = t_{0}}^{t_{0} + t_{m}} \frac{1}{N} \sum_{i = 1}^{N} x_{i, s}^{2} (t)

respectively, where

t_{m} ≫ 1

and

t_{0}

is a sufficiently long time to allow the dynamics to relax to stationarity. Fig. 2 shows a heat map of their averages

\hat{M} = S^{- 1} \sum_{s = 1}^{S} | {\hat{M}}_{s} |

(panel (A)) and

\hat{C} (0) = S^{- 1} \sum_{s = 1}^{S} {\hat{C}}_{s} (0)

(panel (B)), where the absolute value

| . . . |

is used to avoid the cancellation between positive and negative states of the system. Simulations are performed with

S = 1000

realizations, for a system with size

N = 1000

, for different values of the parameters

(J_{0} / J, 1 / g J)

. Notably, the phase diagram derived for fixed-point solutions captures very well the phase transition behaviour of

\hat{C} (0)

and

\hat{M}

, which vanish in the disordered (paramagnetic) phase and become non-zero in the ordered (spin-glass and ferromagnetic) phases. In addition, their asymptotic behaviour close to criticality matches the scaling theoretically predicted for fixed-point solutions (see Fig. 2 panels (C) and (D) and SM, Sec. S.III for details). This suggests that close to the instability line of the paramagnetic phase, trajectories are at a fixed point. This is consistent with bifurcation analysis, showing that the first non-trivial steady-state solutions to appear below criticality are fixed-point solutions (see SM, Sec. S.V B for details). In the next subsection, we study in detail non-fixed-point solutions at stationarity.

4.1.2. Non-fixed point steady states

Let us now focus on stationary solutions other than fixed points. These are characterized by a constant first moment

M (t) = M

and time-translation invariant correlator

C (t, s) = C (t - s)

, where

M

is given by (22)

M = \int D ψ tanh (J_{0} g M + g J \sqrt{C (0)} ψ)

with

C (0)

denoting the equal-time correlator, while the two-time correlator

C (τ)

satisfies (23)

(1 - \partial_{τ}^{2}) C (τ) = Ξ (C (τ), C (0), M)

with (24)

Ξ (C, C (0), M) = \int D_{ϕ ϕ^{'}} tanh [J_{0} g M + ϕ] tanh [J_{0} g M + ϕ^{'}]

where (25)

D_{ϕ ϕ^{'}} = \frac{d ϕ d ϕ^{'}}{2 π g^{2} J^{2} \sqrt{det C (τ)}} exp (- \frac{ϕ^{T} C^{- 1} (τ) ϕ}{2 J^{2} g^{2}}),

see SM, Sec. S.V A for details.

For

C

to be positive definite, we must have

C^{2} (0) \geq C^{2} (τ)

, i.e.

C (0) \geq | C (τ) |

, as expected for a stationary process, where the correlation function (at any time-lag

τ

larger than

0

) is bounded to be smaller than the variance. An expression for the variance

C (0)

is provided in the SM (Sec. S.V A), however, let us remark that it requires the knowledge of the time-dependent (or “non-persistent”) order parameter

C (τ)

, consistently with earlier findings in asymmetrically diluted recurrent neural networks [56]. Hence, the equations for the time-invariant (or “persistent”) order parameters

M

and

C (0)

are not a closed set (not even at stationarity), when the system is not at a fixed point. Only when the system is at a fixed point, i.e.

x (t) = x

, then

C (0) = q

and the equations for

M

and

q

form a closed set.

Following [57], we note that Eq. (23) can be cast in a gradient-descent equation (26)

\partial_{τ}^{2} C (τ) = - \frac{\partial V (C | C (0), M)}{\partial C}

on the potential (27)

V (C | C (0), M) = - \frac{C^{2}}{2} + \int_{0}^{C} d C^{'} Ξ (C^{'}, C (0), M)

which depends on the persistent order parameters

M

and

C (0)

Note that for

J_{0} = 0

, the previous potential reduces to the one governing the correlation function in the model of SCS [29]. Although such an equivalence might seem surprising — given that in the model of SCS, unlike ours, the sum over

j

in Eq. (1) appears outside the hyperbolic tangent — it can be rationalized by a simple argument as explained in Appendix A.

As already noticed in [57], solutions of Eq. (26) conserve the total energy (28)

E = \frac{1}{2} {\dot{C}}^{2} + V (C | C (0), M) .

Moreover, at stationarity,

C (τ) = C (- τ)

, so that

\dot{C} (0) = 0

and, hence, the initial kinetic energy vanishes. Thus, from energy conservation one has (29)

\frac{1}{2} {\dot{C}}^{2} (τ) + V (C | C (0), M) = V (C (0) | C (0), M) .

Physical solutions must be bounded, i.e.

| C | \leq C (0)

, this requires

V^{'} (C (0) | C (0), M) > 0

and

V^{'} (- C (0) | C (0), M) < 0

at the two boundaries

C = \pm C (0)

, respectively (since

\dot{C} (0) = 0

). It can be easily shown that, for any

M

V^{'} (C (0) | C (0), M) < 0

C (0) > q

, hence such initial conditions are unphysical: the system will have to select a stationary value of the equal-time correlator that is smaller or equal than the persistent parameter

q

Let us now study separately the cases (a)

M \neq 0

and (b)

M = 0

4.1.2.1. For $M \neq 0$

One can show (see SM S. V C) that the potential is a monotonic non-decreasing function for any initial condition

C (0) \leq q

and that there is only one stationary point, at the boundary of the physical region

C = C (0) = q

(see Fig. 3, panels (A)). Therefore, the only possible bounded steady-state solutions for

M \neq 0

are fixed-point solutions with

C = q

To verify this, we ran

S

simulations of the microscopic dynamics to obtain trajectories

x_{s} (t)

, with

s = 1 \dots S

and computed, for each realization

s

, the equal-time correlator at stationarity

{\hat{C}}_{s} (0)

, as defined in Eq. (21). The resulting distribution (30)

P (C_{0}) = \frac{1}{S} \sum_{s = 1}^{S} δ (C_{0} - {\hat{C}}_{s} (0))

is plotted in Fig. 3(B, C) as a heat map, for different values of

1 / g J

, two different network sizes

N

, and a fixed value of

J_{0} / J

(corresponding to the ferromagnetic region, i.e.

M \neq 0

). Results show that, as the network size

N

is increased, the probability density concentrates on the curve corresponding to

C (0) = q

(solid line), confirming that, indeed, the value of

C (0)

selected dynamically by the system coincides with

q

4.1.2.2. For $M = 0$

One can show that — as illustrated in Fig. 3 panel (D) —

V

has the shape of a double-well potential for any initial condition

C (0)

smaller than

q

but larger than a certain threshold value

C_{th}

solution of (31)

\int D ψ {tanh}^{2} (g J \sqrt{C_{th}} ψ) = 1 - \frac{1}{g J},

whereas below such a threshold (

C (0) < C_{th}

) it becomes a single-well potential (see SM S. V C for further details).

Therefore, any initial condition

0 \leq C (0) \leq q

leads to bounded solutions, but — depending on the value of

C (0)

that is dynamically selected by the system — different steady-state solutions may arise as summarized in the following:

•
For $0 < C (0) < C_{th}$ the potential has a single-well shape, so that the motion is periodic between $C (0)$ and $- C (0)$ .
•
For $C_{th} < C (0) < q$ , the potential has a double-well shape with a local maximum at $C = 0$ , and three possible types of dynamics can emerge:
- –
  (i) periodic motion in one of the single wells corresponding to motion below the separatrix curve ( $C (0) > C_{0}^{*}$ such that $V (C (0) | C (0), 0) < 0$ );
- –
  (ii) asymptotic motion towards $C = 0$ , corresponding to motion on the separatrix curve ( $C (0) = C_{0}^{*}$ such that $V (C_{0}^{*} | C_{0}^{*}, 0) = 0$ );
- –
  (iii) periodic motion between $C (0)$ and $- C (0)$ , i.e. covering the two wells, corresponding to motion above the separatrix curve ( $C (0) < C_{0}^{*}$ such that $V (C (0) | C (0), 0) > 0$ ).
•
For $C (0) = q$ , the system remains at the fixed point $C (τ) = q$ for all values of $τ$ , i.e. the solution is a fixed point.

In Fig. 3 (E, F), we show

C_{th}

C_{0}^{*}

and

q

as a function of

1 / g J

and

C (0)

(red, green and yellow solid lines, respectively), at a fixed value of

J_{0} / J

(corresponding to the spin-glass region,

M = 0

) and for two different system sizes (panel (E),

N = 100

; panel (F),

N = 1000

The question that remains to be answered is: what is the value of the equal-time correlation,

C (0)

, that is selected by the system at stationarity? To ascertain this, we plot, in the same figure, the probability distribution function Eq. (30), computed from

S = 1000

simulations, versus

1 / g J

, as a heat map, for two different system sizes. These results clearly reveal that, as the system size

N

increases, the probability density concentrates on the curve

C (0) = C_{0}^{*}

. This means that in the thermodynamic limit the time-average

C_{s} (0)

of any single instance

s

of the system (obtained for a given initial condition and a particular realization of the quenched disorder) settles precisely on the separatrix curve, i.e. at the boundary that delimits the basins of attraction of the two potential minima

\pm \bar{C} (C_{0}^{*})

. Therefore, the solution sits on an unstable stationary state, where any slight perturbation causes it to move away from the separatrix in either of the two adjacent states.

In particular, in any single network instance the instantaneous mean-squared activity after a transitory time

t_{0}

(32)

{\hat{q}}_{s} (t_{0} + τ) = \frac{1}{N} \sum_{i} x_{i, s}^{2} (t_{0} + τ)

fluctuates in time (

τ

) about the stationary value

{\hat{C}}_{s} (0)

(unless the system is at a fixed-point steady state). Therefore, in the case in which

{\hat{C}}_{s} (0)

is sufficiently close to

C_{0}^{*}

, such fluctuations can effectively move the system above and below the separatrix as illustrated in Fig. 4. This leads to a motion that is highly irregular and sensitive to perturbations and small changes in initial conditions, as shown in Fig. 5(A). As a result, the motion turns out to be chaotic. More specifically: the time-lagged correlator, computed as a sliding time window average over a single dynamical trajectory,² (33)

{\hat{C}}_{s} (τ) = \frac{1}{t_{m}} \sum_{t = t_{0}}^{t_{0} + t_{m}} \frac{1}{N} \sum_{i} x_{i, s} (t) x_{i, s} (t + τ)

does not exhibit stationary behaviour, but an irregular motion which is not time-translation invariant (TTI), i.e. it retains a dependence on

t_{0}

, see Fig. 5(A), bottom. Indeed, our analysis in Appendix B shows that, in the thermodynamic limit, stationary solutions are unstable against perturbations that break TTI, suggesting that such non-TTI steady states are selected in large networks.

On the other hand, when

C_{s} (0)

is far from

C_{0}^{*}

, small fluctuations cannot shift the system across the separatrix and the resulting motion is periodic, as illustrated in Fig. 5(B, C). As a consequence, in finite networks of size

N

, the system exhibits periodic dynamics in wide regions of the phase diagram, where

C_{s} (0)

is far from

C_{0}^{*}

(see Fig. 3(E)). In addition, Fig. 3(E) shows that, at finite

N

, there is a finite probability that the system exhibits fixed-point dynamics, at sufficiently high values of

1 / g J

To illustrate all this, in Fig. 6 we show a heat map of the probability of finding fixed-point solutions, versus

J_{0} / J, 1 / g J

, in a system of (relatively small) size

N = 100

, as well as typical trajectories exhibited by the system at representative points of the phase diagram. Chaotic trajectories are observed only at sufficiently low values of

1 / g J

, where the probability distribution function of

C (0)

becomes peaked around

C_{0}^{*}

, as shown in Fig. 3. Furthermore, we note that, while chaotic trajectories seem to arise, in finite networks, from fluctuations around the separatrix, this does not seem a sufficient condition to have chaotic motion. For example, when fluctuations of the mean-squared activity around the separatrix value are periodic and in sync with the natural oscillations of the two-time correlator, as shown in Fig. 5(D, E), the system is observed to settle on a periodic trajectory.

In this regard, the dynamics of the two-time correlator resembles that of a forced oscillator, where the external “force” is given by the fluctuations of the initial condition. In particular, if the driving force is periodic and it has the same frequency as the natural frequency of the oscillator, the system oscillates at the frequency of the applied force as shown in Fig. 5(B, D). If the frequency of the applied force is different from the natural frequency of the oscillator, the system may exhibit sub-harmonic oscillations i.e. a periodic motion with a frequency that is a fractional multiple of the input frequency, as shown in Fig. 5(C, E).

As discussed earlier, for large

N

C (0)

is peaked around

C_{0}^{*}

. Therefore, the scenario depicted in Fig. 5(B, C) cannot arise because dynamical fluctuations of the mean-squared activity about its time-averaged value move the system above and below the separatrix curve. However, periodic trajectories as described in Fig. 5 (D, E), can in principle still arise. In order to determine whether periodic motion survives for large

N

, it is therefore important to establish whether such trajectories remain stable in the thermodynamic limit. Our analysis in Appendix B suggests that they lose stability for large

N

, as stationary trajectories in the spin-glass region are unstable, in the thermodynamic limit, against perturbations that break TTI.

Numerical evaluation of the largest Lyapunov exponent (LLE)

Λ

, in simulated trajectories of the system (see Appendix D.3) shows that, upon increasing

N

, only chaotic orbits prevail inside the spin-glass region. Fig. 7 (A) shows that at small system sizes,

Λ

is positive for small values of

1 / g J

within the spin-glass region, and it decreases to zero as

1 / g J

increases. However, as the system size increases, the value of

Λ

1 / g J = 1

increases and eventually vanishes (see inset of Fig. 7, panel B), suggesting that, in this limit,

Λ

remains positive for any

1 / g J < 1

. On the other hand, in the ferromagnetic region

Λ

remains negative for any value of

1 / g J

, see Fig. 7 (C). The inset in panel (B) indicates that the approach to zero of

Λ

1 / g J = 1

is slow, suggesting that in large but finite systems, where

C (0)

has already converged to the separatrix value

C_{0}^{*}

, periodic motion as described in Fig. 5 (D, E) can still be observed. However, as explained earlier, the stability analysis of stationary trajectories carried out in Appendix B suggests that such a periodic motion eventually loses its stability.

In summary, we have found two types of criticality, one SG like (type-II), and the other ferromagnetic like (type-I). All across the SG phase the system exhibits chaotic behaviour in the thermodynamic limit, while in finite networks fixed points and periodic oscillations may emerge, especially close to the transition line.

4.2. Switching on fluctuations: Noisy dynamics ( $σ \neq 0$ ) with uncorrelated interactions ( $γ = 0$ )

In this section, we analyse the effect of additional sources of external variability, particularly the presence of an external noise (i.e.,

σ \neq 0

), on the firing-rate model for uncorrelated networks (

γ = 0

). In this case, there cannot be any fixed-point solution, as the rates

x_{i} (t)

fluctuate stochastically. Non-fixed point steady states are characterized by the same magnetization

M

as in the noiseless dynamics, given by Eq. (22), while the equation for the correlator needs to be modified to (see SM, Sec. S.V A) (34)

\partial_{τ}^{2} C = - \frac{\partial V (C | C (0), M)}{\partial C} - 2 σ^{2} δ (τ) .

Integrating over

τ \in [0, ϵ]

, with

ϵ ≪ 1

leads to (35)

\dot{C} (ϵ) = \dot{C} (0) - σ^{2} + O (ϵ)

and recalling that

C (τ) = C (- τ)

implies

\dot{C} (0) = 0

, one gets for

ϵ \to 0^{+}

that the initial velocity is

\dot{C} (0^{+}) = - σ^{2}

, as previously found in [50]. From now on, for the sake of simplicity, we focus on the case

M = 0

Let us remark that, as noticed earlier,

V^{'} (C (0) | C (0), 0) < 0

for

C (0) > q

, so that this led to un-physical solutions in the noiseless case, where

\dot{C} (0) = 0

. However, the situation is different in the presence of noise; given that the initial velocity is negative, a physical solution can still be obtained if the total energy is equal to zero, which corresponds to the separatrix curve, i.e. asymptotic motion towards

C = 0

(see Fig. 8(B)). This occurs for the initial condition

C (0) = C_{σ}^{*}

such that (36)

V (C_{σ}^{*} | C_{σ}^{*}, 0) = - \frac{1}{2} σ^{4} .

On the other hand, any

C (0) \leq q

is physically plausible as it leads to bounded motion. To ascertain the value of

C (0)

that is dynamically selected by the system, we plot in Fig. 8(C) the distribution

P (C_{0})

, defined in Eq. (30), computed from

S = 1000

simulations of the microscopic dynamics, versus

1 / g J

, as a heat map, as well as the value of

C_{σ}^{*}

and

q

, as defined in Eqs. (17), (36), respectively. Results show that

P (C_{0})

concentrates on the curve

C (0) = C_{σ}^{*}

. As we explained above, while for

C (0) > q

this is the only physical solution, for

C (0) \leq q

, there are in principle many possible solutions. A stability analysis of steady-state solutions (see Appendix B) shows however that such solutions are all unstable, except for

C (0) = C_{σ}^{*}

, which is marginally stable. Then, the system selects the one laying on the separatrix curve, which delimits different possible (unstable) motions. As any fluctuation leads to an instability, chaotic behaviour is expected to emerge in single instances for

C_{σ}^{*} < q

. The critical line

C_{σ}^{*} = q

is thus predicted to mark the transition to chaos. In Fig. 9(A) we plot

q

and

C_{σ}^{*}

versus

1 / g J

, for different values of

σ

. The intersection between

C_{σ}^{*}

and

q

, obtained from (37)

C_{σ}^{*} = \int D ψ {tanh}^{2} (g J \sqrt{C_{σ}^{*}} ψ)

gives a curve in the

(σ^{2}, 1 / g J)

-plane plotted in Fig. 9(B). A similar condition to Eq. (37) has been derived in [50] for the model originally defined by SCS in [29]. Let us remark that the variant of the model we consider here has a narrower chaotic region than the original model, in spite of the fact that their behaviour is identical in the absence of noise. This observation can be rationalized by a simple argument (see Appendix A) which shows that in our variant of the model the noise strength

σ^{2}

is effectively amplified by a factor

{(g J)}^{2}

, which is larger than

1

, in the spin glass region. In the inset of Fig. 9(B) we show the equivalent curve

g^{*} J

as a function of

σ g^{*} J

, for fixed

J

, recovering precisely the result obtained in [50].

The emergence of chaos in the noisy dynamics is confirmed by numerical analysis of the LLE, shown in Fig. 10. For

J_{0} / J = 0.5

, the LLE

Λ

is plotted as a function of

1 / g J

for different values of

σ^{2}

and fixed

N = 1000

. The critical values

1 / g^{*} J

at which the

Λ

becomes positive, i.e.

Λ (1 / g^{*} J) = 0

, are plotted in the inset, which are in agreement with the critical line defined in Eq. (37) (solid, red line).

Thus, the size of the chaotic phase is reduced by the presence of external noise, so that the boundary of the chaotic phase is shifted towards larger coupling values, in agreement with previously reported results.

4.3. Switching on correlations ( $γ \neq 0$ )

In this section, we analyse the dynamics when interactions are correlated (

γ \neq 0

). Any degree of correlations introduces a so-called “retarded” self-interaction in the equation for the effective single-neuron Eq. (10) i.e. a term non-local in time, which makes the dynamics non-Markovian and the analysis considerably harder. For this reason, here, we restrict our analysis to fixed-point solutions of the noiseless dynamics (

σ = 0

In the fixed-point regime, the response function is time-translation invariant

R (t, t^{'}) = R (t - t^{'})

and

C (t, t^{'})

becomes independent of

t, t^{'}

, so that it can be written as

q = C (t, t^{'})

. As before, each realization of

ϕ (t)

becomes a static Gaussian random variable with zero mean and variance

g^{2} J^{2} q

. Setting

ϕ = g J \sqrt{q} ψ

, where

ψ

is a zero-average random variable with unit variance, the fixed point of Eq. (10) satisfies, for each realization of

ψ

, (38)

x (ψ) = tanh [J_{0} g M + γ J^{2} g^{2} χ x (ψ) + g J \sqrt{q} ψ]

where

χ = \int_{0}^{\infty} d τ R (τ)

is the integrated response, which is found from Eq. (14) as (39)

χ = 〈\frac{\partial x (ϕ)}{\partial ϕ}〉 \equiv \frac{1}{g J \sqrt{q}} 〈\frac{\partial x (ψ)}{\partial ψ}〉 .

Averaging Eq. (38) over

ψ

, we obtain (40)

M = 〈 x (ψ) 〉 = \int D ψ tanh [J_{0} g M + γ J^{2} g^{2} χ x (ψ) + g J \sqrt{q} ψ],

and similarly (41)

q = 〈 x^{2} (ψ) 〉 = \int D ψ {tanh}^{2} [J_{0} g M + γ J^{2} g^{2} χ x (ψ) + g J \sqrt{q} ψ] .

As before, for

g = 0

the system is in the paramagnetic phase

M = 0

q = 0

and we expect this solution to remain stable below a critical value of

g

. However, in contrast with the case

γ = 0

, the equations for the moments

M

and

q

do not close, as they retain a dependence on the realizations

x

of the stochastic process. Nevertheless, progress can be made close to the paramagnetic instability, where

tanh (x)

can be linearized and equation closure can be achieved. In Appendix C we show that the critical line where the paramagnetic phase becomes unstable can be calculated explicitly as (42)

\frac{1}{g J} = max (γ + 1, \frac{J_{0}}{J} + γ \frac{J}{J_{0}})

Eq. (42) retrieves the line Eq. (4) for uncorrelated coupling (

γ = 0

) and Eq. (5) for reciprocal interactions (

γ = 1

), and it is consistent with results from random matrix theory. In particular, the first term corresponds to the largest eigenvalue, as derived in [40], while the second term recovers the expression for the outlier eigenvalue in matrices with a non-vanishing mean, derived in [41]. This underscores the versatility and significance of the derived expression across different analytical contexts. Also, as noted in [62], the largest eigenvalue of the bulk and the outlier are each given by the sum of their values at

γ = 0

and a term which is proportional to

γ

. In the context of our analysis, this term can be rationalized as the additional contribution that

M

and

q

receive from the susceptibility (i.e. the integrated response), around the instability line. This relation shows that positively correlated interactions (

γ > 0

) enlarge the regions of ferromagnetic and spin-glass order, with respect to the case of uncorrelated interactions (

γ = 0

), while negatively correlated interactions (

γ < 0

) shrink the extension of the ordered phases, with the spin-glass region disappearing at

γ = - 1

Finally, note that, owing to the non-linear and non-closed nature of equations, analyses away from criticality cannot be carried out exactly. Thus, we resort to numerical simulations.

In Fig. 11 we show a heat map of the time averaged

\hat{M}

and

\hat{C} (0)

(computed from numerical simulations) in the parameter space

J_{0} / J, 1 / g J

, together with the theoretically determined critical line Eq. (42), for different values of

γ

(computed from the stability condition of fixed-point solutions). Observe that, even if non-fixed-point solutions are expected to appear in the ordered phase when interactions are non-symmetric (i.e.

γ < 1

), the theoretical lines obtained for fixed points are seen to capture well the actual phase-transitions for all values of

γ

. The figure clearly confirms the theoretical analyses: the spin-glass phase is enlarged for positively correlated couplings and progressively shrinks for negatively correlated ones (disappearing in the limit case

γ = - 1

) and the theoretical lines correctly reproduce the actual phase transition in all cases.

Also, in Fig. 12 we show results for the LLE as a function of the coupling strength (

1 / g J

), computed for different values of

γ

. The curves show that, for all values of

γ

, the exponents become positive below some critical value (marked with diamond symbol) that is consistent with the critical line

1 / g_{c} J = 1 + γ

. Finite-size effects are quantified in the inset (for

γ = - 0.15

). This shows that chaotic behaviour extends to the whole spin-glass region for all values of

γ

Some remarks are in order:

•
The critical value (in terms of $1 / g J$ ) decreases as $γ$ is reduced, so that it is maximal for $γ = 1$ (in agreement with the theoretical prediction).
•
Smaller values of $γ$ have larger LLEs in the limit of large coupling values (small $1 / g J$ ’s), so that, in this regime, the dynamics becomes more chaotic for anti-correlated couplings than for positively correlated ones.
•
Indeed, as the value of $γ$ approaches $1$ , the LLE becomes very close to $0$ all across the spin-glass phase, so that, actually, the dynamics becomes marginally stable in such a limit. This last result resembles the recent observation of a phase of marginal stability in fully-correlated (i.e. symmetric) random networks in models for theoretical ecology [18]. Dynamics in such region can be regarded as “at the edge of chaos”, as defined in other contexts see e.g. [11], [63], [64], [65].

In summary, the introduction of correlations in the couplings quantitatively modifies the shape of the phase diagram with respect to the uncorrelated case: the spin-glass region is enlarged for positive correlations and shrinks for anticorrelated ones. On the other hand, the “strength” of chaos is larger for anticorrelations and tends to be reduced by positive correlations, leading to marginal stability in the limit of perfectly correlated couplings.

5. Discussion and conclusions

In this study we have considered a simple neural network rate model with random Gaussian non-reciprocal interactions, having non-vanishing mean and possible correlations. We have analysed the dynamics averaged over the network ensemble in the thermodynamic limit using a path integral formalism. These analytical studies have been complemented with extensive numerical simulations finite networks.

In general, much like in the physics of standard (reciprocal) disordered systems, fixed-point solutions are characterized by two order parameters, namely the mean (

M

) and the variance (

q

) of the neural activity. Thus, two different types of criticality emerge depending on whether the mean becomes non-zero (type-I) or not (type-II) at the transition, as the strength of the coupling is increased.

These criticalities correspond to the transition of the system from the paramagnetic phase to either a ferromagnetic (type-I) or spin-glass (type-II) phase, respectively.

Our analyses show that as soon as some degree of non-reciprocity is switched on —i.e., for any value of the coupling correlations, from nearly symmetric to anti-symmetric,

- 1 < γ < 1

— there is a region, within the paramagnetic phase, where the network is reactive so that initial perturbations can be amplified before relaxing back to the fixed point solution. Moreover, for all values of the coupling correlations, a chaotic behaviour emerges within the spin-glass phase (

M = 0

q \neq 0

). Furthermore, we have shown how the presence of correlations (i.e. some degree of reciprocity) in the couplings changes the phase diagram: negative correlations shrink the spin-glass phase (where chaotic dynamics emerges) while positive correlations enlarge it. On the other hand, negative correlations increase the strength of chaos, while positive correlations reduce it, with chaos disappearing for symmetric couplings. In this case, equilibrium is attained, and the whole spin-glass phase exhibits marginal stability, similar to what is observed in [18]. We have demonstrated that the extent of the chaotic phase diminishes with the introduction of external noise, indicating that noise tends to suppress chaos, consistent with findings from previous studies [50].

Chaotic phases observed in networks characterized by non-reciprocal couplings represent an intrinsic non-equilibrium phenomenon, yet they exhibit intriguing analogies with equilibrium spin-glass phases. For instance, our analysis demonstrates that within the chaotic phase, the dynamical equation governing the global two-time correlation function permits the existence of multiple steady states, encompassing fixed points and periodic orbits. However, our analyses reveal that the system dynamically selects the steady state that corresponds to motion along the separatrix curve, which delimits the basins of attractions of different types of solutions. This behaviour is reminiscent of equilibrium systems in spin-glass (fRSB) phases, where a multitude of steady states emerges and the system is observed to select marginally stable ones, i.e. saddles in the free-energy landscape, that are linked together by flat directions. Our analyses also suggest that chaotic motion is the manifestation, at the level of single network instances, of an ensemble-averaged dynamics that lies on the separatrix curve delimiting different possible steady states. Such a motion is not time-translation invariant (TTI), similarly to ageing dynamics in spin-glass phases.

This last observation led us to propose a method for identifying the initiation of chaotic behaviour. This method relies on a linear stability analysis of the steady-state solutions with respect to perturbations that disrupt time translation invariance. Our method is consistent with the two-replica approach described in [50], [57], but it is considerably simpler, and it shows that chaotic behaviour emerges whenever the equation for the global correlation function allows for multiple steady-states.

The sensitive dependence on initial conditions in the chaotic phase implies that two “replicas” — two copies of the system with the same disorder but initialized with different configurations — exhibit distinct dynamical behaviours, resembling replicas of an equilibrium system ending up in different states. Furthermore, we find that while the global correlator exhibits multiple steady states in the chaotic phase, the presence of both a signal (i.e., a non-vanishing mean connectivity) and noise drives the system to a phase where only one steady state exists. This, again, mirrors what happens in equilibrium settings, where signal and thermal noise are known to induce transitions from the spin-glass phase to the ferromagnetic and paramagnetic phase respectively.

A theoretical prediction for the transition line between the spin-glass and the ferromagnetic region, for general values of the coupling asymmetry as well as a derivation of the asymptotic behaviour of the order parameters near criticality are currently missing. In particular, as we have seen, when interactions are not uncorrelated (or fully asymmetric), the equations for the time-dependent moments of the stochastic process do not close. This prevents an exact analysis away from the paramagnetic instability, to address the previous questions. However, one may be able to use closure schemes (e.g. Gaussian closure) as an approximation to further characterize the phase diagram of fixed-point solutions. In addition, within such approximation schemes, one might be able to obtain a closed equation for the two-time correlation at stationarity, which can be formulated as a gradient-descent on a potential. This approach can extend the analysis conducted for fully asymmetric couplings to interactions with arbitrary asymmetry. Another outlook could involve adapting the methodologies outlined in [66], [67] to either numerically sample or analytically calculate fluctuations in the largest Lyapunov exponent. This approach could help identify model parameters where such fluctuations become unusually large, indicating the onset of chaos.

Finally, it would be important to further investigate whether chaotic dynamics and its links with spin-glass phases survive when different couplings, other than Gaussian, are chosen. For example, for applications in neuroscience, it would be interesting to carry out a similar analysis for couplings which satisfy Dale’s rule [37] or are non-Gaussian [68], whereas the choice of sparse couplings may be relevant to applications in biology (e.g. gene regulatory networks [11], [69]) and theoretical ecology (e.g. Lotka–Volterra models [18], [20], [43]).

CRediT authorship contribution statement

Carles Martorell: Writing – original draft, Validation, Software, Investigation, Formal analysis. Rubén Calvo: Validation, Software, Investigation, Formal analysis. Alessia Annibale: Writing – review & editing, Writing – original draft, Supervision, Methodology, Formal analysis, Conceptualization. Miguel A. Muñoz: Writing – review & editing, Supervision, Funding acquisition, Formal analysis, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

We acknowledge the Spanish Ministry and Agencia Estatal de investigación (AEI) through Project of I+D+i Ref. PID2020-113681GB-I00, funded by MICIN/AEI/10.13039/501100011033 for financial support. We also acknowledge the support from the Consejería de Conocimiento, Investigación Universidad, Junta de Andalucía and European Regional Development Fund , Project reference P20-00173. AA wishes to thank Carlos-I Institute of Theoretical and Computational Physics at the University of Granada, for the kind hospitality.

Appendix A. On the equivalence between our model and the classic one (SCS)

In this section we show that the model originally introduced in [29] can be mapped to the variant we introduce here, when interactions are uncorrelated and zero-averaged, upon suitably scaling the noise strength. Multiplying Eq. (1), (A.1)

{\dot{x}}_{i} = - x_{i} + tanh (g \sum_{j} J_{i j} x_{j}) + ξ_{i},

times

g J_{ℓ i}

and summing over

i

, we have, upon defining

y_{ℓ} = g \sum_{i} J_{ℓ i} x_{i}

and

ϕ_{ℓ} = g \sum_{i} J_{ℓ i} ξ_{i}

(A.2)

{\dot{y}}_{ℓ} = - y_{ℓ} + g \sum_{i} J_{ℓ i} tanh (y_{i}) + ϕ_{ℓ}

where (A.3)

〈 ϕ_{ℓ} (t) ϕ_{k} (t^{'}) 〉 = g^{2} \sum_{i j} J_{ℓ i} J_{k j} 〈 ξ_{i} (t) ξ_{j} (t^{'}) 〉 = 2 g^{2} σ^{2} \sum_{i j} J_{ℓ i} J_{k j} δ_{i j} δ (t - t^{'}) = 2 g^{2} σ^{2} \sum_{i} J_{ℓ i} J_{k i} δ (t - t^{'})

For uncorrelated, zero-mean Gaussian interactions

J_{i j} \sim N (0, J^{2} / N)

, we have (A.4)

〈 ϕ_{ℓ} (t) ϕ_{k} (t^{'}) 〉 = 2 g^{2} σ^{2} J^{2} δ_{k ℓ} δ (t - t^{'}),

hence the two models are fully equivalent when rescaling

σ

with

g J

. If, however, the

J_{i j}

’s are correlated or non-zero averaged, then the two models are no longer equivalent.

Appendix B. Stability of steady-state solutions with uncorrelated interactions ( $γ = 0$ )

In this section we study the linear stability of non-fixed points steady-states solutions. Away from stationarity, we can write the equation of motion for the two-time correlator

C (t, s)

as (B.1)

(1 + \partial_{t}) (1 + \partial_{s}) C (t, s) = F (C (t, s), C (t, t), C (s, s), C (s, t), M (t), M (s))

where (B.2)

F (C (t, s), C (t, t), C (s, s), M (t), M (s)) = \int D_{ϕ, ϕ^{'}} (t, s) tanh [J_{0} g M (t) + ϕ] tanh [J_{0} g M (s) + ϕ^{'}]

where the Gaussian kernel is defined as (B.3)

D_{ϕ, ϕ^{'}} (t, s) = \frac{d ϕ d ϕ^{'}}{2 π g^{2} J^{2} \sqrt{det M (t, s)}} exp (- \frac{ϕ^{T} M {(t, s)}^{- 1} ϕ}{2 J^{2} g^{2}})

and

M (t, s) = (\begin{matrix} C (t, t) & C (t, s) \\ C (s, t) & C (s, s) \end{matrix})

As noted above, in Sections 4.1 Uncorrelated interactions (, 4.2 Switching on fluctuations: Noisy dynamics (, chaotic motion is not time translation invariant, as trajectories are highly sensitive to their initial conditions; i.e. even small initial perturbations can lead to vastly different outcomes as time progresses.

In what follows, we examine the stability of stationary trajectories under small perturbations

δ (t, s)

that break time-translation invariance. Recalling that stationary solutions are characterized by a constant first moment

M (t) = M

and time-translation invariant correlator

C (t, s) = C (t - s)

, we can write (B.4)

C (t, s) = C (τ) + ϵ δ (t, s)

(B.5)

C (t, t) = C (s, s) = C (0)

(B.6)

M (t) = M

where

τ = t - s

, i.e. we consider a non-stationary regime where one-time quantities remain constant and two-time quantities break time-translation invariance — much as happens in glassy regimes.

Inserting Eqs. (B.4), (B.5), (B.6) in the equation for the two-time correlator

C (t, s)

, given by Eq. (B.1), Taylor expanding the right hand side for small

ϵ

to linear order and using

C (τ) = C (- τ)

and (B.7)

F (C (τ), C (0), C (0), C (- τ), M, M) \equiv Ξ (C (τ), C (0), M),

we get (i) to

O (ϵ^{0})

terms, Eq. (23); and (ii), to

O (ϵ)

, the following equation for the deviations from stationarity (B.8)

(1 + \partial_{t}) (1 + \partial_{s}) δ (t, s) = δ (t, s) {\frac{\partial Ξ (C, C (0), M)}{\partial C}|}_{C (τ)} .

Following [50], we set

T = t + s

and

τ = t - s

and denote

δ (t, s) = K (T, τ)

to rewrite the above as (B.9)

[{(\partial_{T} + 1)}^{2} - \partial_{τ}^{2}] K (T, τ) = K (T, τ) Ξ^{'} (C (τ), C (0), M)

The above partial differential equation can be solved via separation of variables, i.e. by searching for a solution in the factorized form (B.10)

K (T, τ) = ψ (τ) e^{κ T / 2}

Inserting this ansatz into Eq. (B.9), a Schrödinger’s equation for

ψ (τ)

is obtained (B.11)

[- \partial_{τ}^{2} - V^{''} (C (τ) | C (0), M)] ψ (τ) = [1 - {(\frac{κ}{2} + 1)}^{2}] ψ (τ)

with quantum mechanical potential

- V^{''} (C | C (0), M)

. As observed in [50] there will be a set of allowed energies

E_{0} < E_{1} < E_{2} \dots

, with

E_{n} = 1 - {(\frac{κ_{n}}{2} + 1)}^{2}

, and associated eigenfunctions

ψ_{n} (τ)

. A stationary solution

C (τ)

will be linearly stable if the perturbation

K (T, τ)

decays with

T

, i.e. if the largest value of

κ

κ_{0} = - 1 + \sqrt{1 - E_{0}}

, is negative. This requires the ground state energy

E_{0}

to be positive.

It is also useful to note that, as pointed out in [50], Eq. (23) implies (B.12)

[- \partial_{τ}^{2} - V^{''} (C (τ) | C (0), M)] \dot{C} (τ) = 0,

hence non-fixed point steady-state solutions

\dot{C} (τ) \neq 0

, correspond — when they exist — to eigenfunctions with energy

E = 0

As for an eigenfunction to exist, the energy must be greater or equal than the quantum-mechanical potential

- V^{''} (C | C (0), M)

; if the latter is positive, all the eigenfunctions, including the ground state, have positive energy, hence the motion is stable.

1.
For $M \neq 0$ , $V^{'} (C | C (0), M)$ is strictly decreasing for any $- C (0) \leq C \leq C (0)$ and $0 < C (0) \leq q$ , hence $V^{''} (C | C (0), M) < 0$ . This implies that for any $C (0) \leq q$ , $E_{n} > 0 \forall n$ and the only allowed solution of Eq. (B.12) is the fixed-point solution $C (τ) = q$ (non-trivial solutions with zero energy are not allowed). Furthermore, this must be stable, as $E_{0} > 0$ .
2.
Conversely, for $M = 0$ , the potential has at least one minimum, either at $C = 0$ or $C = \pm \bar{C} (C (0))$ , for any value of $C (0) \leq q$ , hence Eq. (B.12) allows a non-trivial solution, for any $C (0) < q$ , which is either a periodic orbit or the separatrix curve. These represent eigenfunctions with energy $E = 0$ . Periodic solutions have multiple nodes (as $\dot{C}$ vanishes at every turning point), whereas the separatrix curve has precisely one node in the noiseless case (where $\dot{C} (0) = 0$ ) and no node in the noisy case (where $\dot{C} (0^{+}) = - σ^{2}$ ). As noted in [50], in the noiseless case $\dot{C}$ has at least one node, hence it cannot be the ground state, but it must correspond to one of the excited states $ψ_{n} (τ)$ with $n \geq 1$ , implying that the ground-state energy is $E_{0} < 0$ . Recalling that, in the noiseless case, steady-state solutions are physical only for $C (0) \leq q$ , we have that all physical steady-state solutions are unstable in the noiseless case. On the other hand, in the presence of noise, solutions with $C (0) > q$ are allowed, only for $C (0)$ equating the separatrix value $C_{σ}^{*}$ , which solves Eq. (36). Since this is the only physical solution, it must be stable for $C_{σ}^{*} > q$ . Conversely, any $C (0) < q$ leads to a physical solution. As for the noiseless case, such solutions are all unstable, except the separatrix curve, which is marginally stable (as it has no node for $σ \neq 0$ ). However, as any fluctuation leads to an instability, chaotic behaviour is expected to emerge in single instances for $C_{σ}^{*} < q$ . The critical line $C_{σ}^{*} = q$ is thus predicted to mark the transition to chaos.

Appendix C. Phase diagram for $γ \neq 0$

In this Appendix we derive the phase diagram for the fixed-point solutions of Eq. (10), when interactions are correlated. We note that for

g = 0

, the solution of Eqs. (41)–(42) is

M = 0

q = 0

. It is expected that

(M, q) = (0, 0)

remains stable for

g

below a critical value

g_{c}

, where non-trivial solutions emerge. In this region, Eq. (38) reduces to (C.1)

x (ψ) = tanh (γ g^{2} J^{2} x (ψ))

so that

x (ψ)

is deterministic. We have

x (ψ) = 0

for

γ g^{2} J^{2} < 1

— either

γ > 0

and

{(g J)}^{- 1} > \sqrt{γ}

, or

γ < 0

. For

γ > 0

and

{(g J)}^{- 1} < \sqrt{γ}

, two non-zero solutions emerge

\pm x^{⋆}

, one positive and one negative. Since

x (ψ) = \pm x^{⋆}

implies

q \neq 0

, bifurcations in

x

imply bifurcations in

q

, hence we assume that

x

is small when

q

is small. To locate

g_{c}

, we expand Eqs. (41)–(42) for small

M

q

and

x (ψ)

. Starting with Eq. (41), we have (C.2)

M = \int D ψ [J_{0} g M + γ J^{2} g^{2} χ x (ψ) + g J \sqrt{q} ψ] = (J_{0} g + γ J^{2} g^{2} χ) M

To make progress, we need an expression for

χ

. Setting

ϕ = g J \sqrt{q} ψ

we can write Eq. (38) as (C.3)

x (ϕ) = tanh [J_{0} g M + γ J^{2} g^{2} χ x (ϕ) + ϕ]

Differentiating Eq. (C.3) with respect to

ϕ

we obtain (C.4)

\frac{\partial x}{\partial ϕ} = [1 - {tanh}^{2} (g J_{0} M + γ g^{2} J^{2} χ x + ϕ)] (γ g^{2} J^{2} χ \frac{\partial x}{\partial ϕ} + 1),

and re-arranging (C.5)

\frac{\partial x}{\partial ϕ} {1 - γ g^{2} J^{2} χ [1 - {tanh}^{2} (g J_{0} M + γ g^{2} J^{2} χ x + ϕ)]} = 1 - {tanh}^{2} (g J_{0} M + γ g^{2} J^{2} χ x + ϕ),

so that, writing in terms of

ψ

, we have (C.6)

\frac{1}{g J \sqrt{q}} \frac{\partial x}{\partial ψ} {1 - γ g^{2} J^{2} χ [1 - {tanh}^{2} (g J_{0} M + γ g^{2} J^{2} χ x + g J \sqrt{q} ψ)]} = 1 - {tanh}^{2} (g J_{0} M + γ g^{2} J^{2} χ x + g J \sqrt{q} ψ)

Close to the transition line

M

q

and

x

are small. A leading-order expansion in these quantities gives (C.7)

\frac{1}{g J \sqrt{q}} \frac{\partial x}{\partial ψ} {1 - γ g^{2} J^{2} χ} = 1 .

Finally, averaging over

ψ

and using Eq. (39) we get (C.8)

χ (1 - γ g^{2} J^{2} χ) = 1

The physical solution is (C.9)

χ = \frac{1 - \sqrt{1 - 4 γ g^{2} J^{2}}}{2 γ g^{2} J^{2}}

(as the other diverges when

g \to 0

). Substituting in Eq. (C.2), we have that

M

bifurcates at (C.10)

g J_{0} + γ g^{2} J^{2} \frac{1 - \sqrt{1 - 4 γ g^{2} J^{2}}}{2 γ g^{2} J^{2}} = 1

which gives (C.11)

\frac{1}{g J} = \frac{J_{0}}{J} + γ \frac{J}{J_{0}}

Similarly, expanding the equation for

q

, Eq. (42), and setting

M = 0

we obtain (C.12)

q = \int D ψ {[J_{0} g M + γ J^{2} g^{2} χ x (ψ) + g J \sqrt{q} ψ]}^{2} = γ^{2} J^{4} g^{4} χ^{2} q + g^{2} J^{2} q + 2 γ {(g J)}^{3} χ \sqrt{q} \int D ψ x (ψ) ψ = (γ^{2} J^{4} g^{4} χ^{2} + g^{2} J^{2} + 2 γ {(g J)}^{4} χ^{2}) q

where we have integrated by parts (C.13)

\int D ψ x (ψ) ψ = \int \frac{d ψ}{\sqrt{2 π}} (- \frac{\partial}{\partial ϕ} e^{- ψ^{2} / 2}) x (ψ) = \int \frac{d ψ}{\sqrt{2 π}} e^{- ψ^{2} / 2} \frac{\partial}{\partial ψ} x (ψ) = 〈\frac{\partial x}{\partial ψ}〉 = g J \sqrt{q} χ .

Substituting Eq. (C.9) in Eq. (C.12) we see that bifurcations in

q

occur at (C.14)

1 = \frac{γ + 2}{2 γ} (1 - \sqrt{1 - 4 γ g^{2} J^{2}} - 2 γ g^{2} J^{2}) + g^{2} J^{2}

which can be simplified to (C.15)

\frac{2}{g^{4} J^{4}} - (2 γ^{2} + 3 γ + 2) \frac{1}{g^{2} J^{2}} - γ {(γ + 1)}^{2} = 0 .

As the physical solution must be positive, it follows that (C.16)

\frac{1}{g J} = \frac{1}{2} \sqrt{(2 γ^{2} + 3 γ + 2) \pm \sqrt{{(2 γ^{2} + 3 γ + 2)}^{2} + 8 γ {(γ + 1)}^{2}}} .

With little algebra, the expression above can be simplified to (C.17)

\frac{1}{g J} = \frac{1}{2} \sqrt{2 + γ (2 γ + 3) \pm | 2 γ^{2} + 5 γ + 2 |} \equiv a (γ) .

2 γ^{2} + 5 γ + 2 \geq 0

for

γ \geq - 1 / 2

, we have (C.18)

a (γ) = \{\begin{matrix} \frac{1}{2} \sqrt{2 + γ (2 γ + 3) \pm (2 γ^{2} + 5 γ + 2)} & γ \geq - 1 / 2 \\ \frac{1}{2} \sqrt{2 + γ (2 γ + 3) \mp (2 γ^{2} + 5 γ + 2)} & γ < - 1 / 2 \end{matrix}

hence there are two possible solutions,

a_{+} (γ) = 1 + γ

and

a_{-} (γ) = \sqrt{- γ / 2}

, the latter existing only for

γ < 0

. As for

γ = - 1

, all the eigenvalues of the interaction matrix

J

are purely imaginary, no (non-trivial) fixed-point solution are expected to bifurcate from

x = 0

, as bifurcating solutions must be purely oscillatory. Hence, it is expected that

a (- 1) = 0

. This suggests that the solution

a_{+} (γ)

is the physical one for any

γ \in [- 1, 1]

. Combining Eq. (C.11) with the result

a (γ) = 1 + γ

, we obtain the critical line where the paramagnetic phase becomes unstable given in Eq. (42).

Appendix D. Numerical procedure

D.1. Computation of trajectories

In simulations, trajectories of Eq. (1) are computed numerically by a simple Midpoint Runge–Kutta method (see, for instance, [70]). The time step is fixed

h = 0.1

, except for the calculation of the LLE for which

h = 0.01

because of the numerical sensitivity of the Bennetin–Wolf procedure. In order to avoid the transient regime, trajectories are integrated up to

t_{m a x} = 2000

units of time. The initial step

x_{0}

and the matrix

J

are selected randomly for each instance of the simulation. In particular, each pair

(J_{i j}, J_{j i})

, for

1 \leq i < j \leq N

, is generated following a multivariate Gaussian distribution with average

μ

and covariance matrix

Σ

, as defined in (S13), respectively. The diagonal elements,

J_{i i}

, are each obtained from a Gaussian distribution with average

J_{0} / N

and variance

J^{2} / N

D.2. Shape of the potential $V (C | C (0), M)$

In order to numerically compute

V (C | C (0), M)

, we simplify the function

Ξ (C, C (0), M)

, defined in Eq. (24), as follows. The Gaussian kernel, Eq. (25), is divided into two parts using the conditional probability formula of a multivariate-Gaussian distribution, i.e. the probability of

ϕ

respect to fixed

ϕ^{'}

. In these terms, one can write (D.1)

\{\begin{matrix} Φ (ϕ, ϕ^{'}) & = g [J (\sqrt{C_{0} - C^{2} / C_{0}} ϕ + C / \sqrt{C_{0}} ϕ^{'}) + J_{0} M] \\ Ψ (ϕ) & = g [J \sqrt{C_{0}} ϕ^{'} + J_{0} M] \end{matrix}

then, by means of the Fubini’s theorem, function

Ξ (C, C (0), M)

becomes (D.2)

Ξ = \int \frac{d ϕ}{\sqrt{2 π}} e^{- ϕ^{2} / 2} tanh (Ψ (ϕ)) \{\int \frac{d ϕ^{'}}{\sqrt{2 π}} e^{- ϕ^{' 2} / 2} tanh (Φ (ϕ, ϕ^{'}))\} .

With this simple form, it is straightforward to compute

V (C | C (0), M)

directly by solving numerically the iterated integrals.

D.3. Computation of the largest Lyapunov exponent

In order to compute the LLE we use the Bennetin–Wolf algorithm [71], [72], that works as follows. Consider a dynamical system described by the vector

x (t) \in R^{N}

at time

t \geq 0

, which evolves following a differentiable dynamical system given by

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

. In order to integrate it numerically, we construct the map

x_{i + 1} = f (x_{i})

where

f

is the discretization of the vector function

F

(which depends on the integration algorithm used) and

{x_{i} : i = 0, 1, 2, \dots}

defines the discrete trajectory of the system. We use a simple integration algorithm, i.e. the Euler’s method: the vector function is defined as the linear evolution of state

x

in a time-step

h

f (x) \equiv x + F (x) h .

Upon defining the matrix

T (x) = J_{x} f

, where

J_{x} f

is the

N \times N

Jacobian matrix of the vector function

f

at state

x

, we have

T (x_{i + 1}) = I + h J_{x_{i + 1}} F

. Writing (D.3)

T_{x_{0}}^{n} = T (x_{n - 1}) \dots T (x_{1}) T (x_{0}),

as a matrix product, the LLE can be computed as the limit (D.4)

Λ = lim_{n ⟶ \infty} \frac{1}{n} log ‖ T_{x_{0}}^{n} u ‖

for a given vector

u

, where

‖ v ‖

denotes the norm of the vector

v

. To implement this method computationally, the Bennetin–Wolf correction [73], [74] needs to be used to avoid the divergence of

T_{x_{0}}^{n} u

for large

n

. The correction consists in normalizing this vector at each step. Hence, fixing a starting point

x_{0}

, and defining an initial vector

u_{0} = [1, \dots, 1] / \sqrt{N}

, the algorithm is described as follows: for

i = 1, \dots, T

1.
Obtain the new state at $i + 1$ , $\begin{matrix} t_{i + 1} = & t_{i} + h \\ x_{i + 1} = & x_{i} + h F (x_{i}) \\ {\tilde{u}}_{i + 1} = & T (x_{i + 1}) u_{i} \end{matrix}$
2.
Normalize the perturbation vector $u_{i + 1} = {\tilde{u}}_{i + 1} / ‖ {\tilde{u}}_{i + 1} ‖$
3.
Compute the LLE at time-step $i + 1$ , $L_{i + 1} = L_{i} + log (‖ {\tilde{u}}_{i + 1} ‖) / t_{i + 1}$

The algorithm gives a set

{L_{i} : i = 0, \dots, T}

of exponents which, in the limit, converge to the LLE of the system,

Λ

Consider now a stochastic differential equation,

\dot{x} (t) = F (x (t)) + ξ (t)

such that

ξ (t)

is a (multi-component) white-noise process with zero mean and variance equal to

2 σ^{2}

. Hence, the previous algorithm should be modified to introduce the stochastic integration: at step 1., the integration of the trajectory should be replaced by an Euler–Maruyama integration step [70], (D.5)

x_{i + 1} = x_{i} + F (x_{i}) h + \sqrt{2 h σ^{2}} u_{i}

where

u_{i}

, for

i = 1, 2, \dots

, is a set of i.i.d. random Gaussian variables with zero mean and unit variance.

For the dynamical system given by Eq. (1), it is sufficient to fix

h = 0.01

, and integrate the system up to

t_{m a x} = 200

. As an example, in Fig. 13 we show

S = 100

realizations of the algorithm for a system of size

N = 100

(grey, solid line) and

σ^{2} = 0

, each of them initialized with random conditions,

x_{0}

and

J

. The system is set at

J_{0} / J = 0.5

(i.e. inside the spin-glass phase). We show results for

1 / g J = 0.9

(left panel) and

1 / g J = 0.35

(right panel). The solid, black line indicates the average over different realizations at each time step, showing that the mean activity is stable (negative LLE) for

1 / g J = 0.9

and unstable (positive LLE) for

1 / g J = 0.35

. In the latter case, there exists a large amount of unstable trajectories, such that

L_{i}

is positive, and some stable trajectories, such that

L_{i}

negative, e.g. fixed points.

In order to study the convergence, we compute the LLE as a function of the system size

N

. In Fig. 14 we plot, as an example, the LLE

Λ

, defined as the asymptotic value of

L (t)

(i.e. the value at

t = 200

), averaged over

S = 20

realizations, for the particular case of

σ = 0

. Fixing

J_{0} / J = 0.5

, we study three different values of

1 / g J

in order to characterize the stability of the system in the limit of large system size:

1 / g J = 1.1

(solid green line) corresponding to the paramagnetic phase;

1 / g J = 1.01

(solid blue line) falling close to the critical line separating the spin-glass and the paramagnetic region;

1 / g J = 0.75

(solid orange line), lying in the spin-glass region, away from the critical line. Coloured area shows the standard deviation of

Λ

for each value of

N

. For

1 / g J > 1

(blue and green lines), the LLE is typically negative, indicating the existence of stable, fixed point trajectories. In this region, the standard deviation of the LLE is smaller compared to the

1 / g J < 1

case (orange line). The reason behind this is that, above this critical value of

1 / g J

, the system converges towards fixed point solutions; while below it, finite size effects in the network encourage the appearance of many different trajectories: chaotic, limit cycles and even, occasionally, fixed points, all of it contributing to increasing the standard deviation of the LLE. In practice, for large enough systems (i.e.

N > 1000

) the LLE stabilizes considerably and the standard deviation reduces, as it can be seen in Fig. 3.

Appendix E. Supplementary data

The following is the Supplementary material related to this article. Download: Download Acrobat PDF file (474KB)

MMC S1. The supplemental material contains details on the path-integral method and on the solutions of the resulting self-consistent equations. The reactivity bifurcation diagram is discussed.

Data availability

No data was used for the research described in the article.

References

[1]
Mézard M., Parisi G., Virasoro M.A.
Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications
World Scientific Publishing Company (1987)
Google Scholar
[2]
Castellani T., Cavagna A.
Spin-glass theory for pedestrians
J Stat Mech Theory Exp, 2005 (05) (2005), p. P05012
View in Scopus Google Scholar
[3]
Nishimori H.
Statistical physics of spin glasses and information processing: an introduction
Clarendon Press (2001)
Google Scholar
[4]
Charbonneau P., Marinari E., Parisi G., Ricci-tersenghi F., Sicuro G., Zamponi F., et al.
Spin glass theory and far beyond: Replica symmetry breaking after 40 years
World Scientific (2023)
Google Scholar
[5]
Aspelmeier T., Bray A.J., Moore M.A.
Complexity of ising spin glasses
Phys Rev Lett, 92 (2004), Article 087203
Google Scholar
[6]
Cavagna A., Giardina I., Parisi G.
Numerical study of metastable states in ising spin glasses
Phys Rev Lett, 92 (2004), Article 120603
Google Scholar
[7]
Annibale A., Gualdi G., Cavagna A.
Coexistence of supersymmetric and supersymmetry-breaking states in spherical spin-glasses
J Phys A: Math Gen, 37 (47) (2004), p. 11311
View at publisher
Crossref View in Scopus Google Scholar
[8]
Müller M., Wyart M.
Marginal stability in structural, spin and electron glasses
Annu Rev Condens Matter Phys, 6 (1) (2015), pp. 177-200
View at publisherCrossref Google Scholar
[9]
Herz A.V., Hopfield J.J.
Earthquake cycles and neural reverberations: Collective oscillations in systems with pulse-coupled threshold elements
Phys Rev Lett, 75 (1222) (1995), pp. 1222-1225
View in Scopus Google Scholar
[10]
Challet D., Chessa A., Marsili M., Zhang Y.-C.
From minority games to real markets
Quant Finance, 1 (1) (2001), pp. 168-176
View in Scopus Google Scholar
[11]
Stokić D., Hanel R., Thurner S.
Inflation of the edge of chaos in a simple model of gene interaction networks
Phys Rev E, 77 (2008), Article 061917
View at publisherCrossref View in Scopus Google Scholar
[12]
Beggs J.M., Timme N.
Being critical of criticality in the brain
Front Physiol, 3 (163) (2012)
Google Scholar
[13]
Mora T., Bialek W.
Are biological systems poised at criticality?
J Stat Phys, 144 (268) (2011)
Google Scholar
[14]
Muñoz M.A.
Colloquium: Criticality and dynamical scaling in living systems
Rev Modern Phys, 90 (3) (2018), Article 031001
View in Scopus Google Scholar
[15]
Moretti P., Muñoz M.A.
Griffiths phases and the stretching of criticality in brain networks
Nat Commun, 4 (1) (2013), pp. 1-10
Google Scholar
[16]
Martinello M., Hidalgo J., Maritan A., Di Santo S., Plenz D., Muñoz M.A.
Neutral theory and scale-free neural dynamics
Phys Rev X, 7 (4) (2017), Article 041071
View in Scopus Google Scholar
[17]
di Santo S., Villegas P., Burioni R., Muñoz M.A.
Landau–Ginzburg theory of cortex dynamics: Scale-free avalanches emerge at the edge of synchronization
Proc Natl Acad Sci, 115 (7) (2018), pp. E1356-E1365
View in Scopus Google Scholar
[18]
Biroli G., Bunin G., Cammarota C.
Marginally stable equilibria in critical ecosystems
New J Phys, 20 (8) (2018), Article 083051
View at publisher
Crossref View in Scopus Google Scholar
[19]
Altieri A., Franz S.
Constraint satisfaction mechanisms for marginal stability and criticality in large ecosystems
Phys Rev E, 99 (1) (2019), Article 010401
View in Scopus Google Scholar
[20]
Altieri A., Roy F., Cammarota C., Biroli G.
Properties of equilibria and glassy phases of the random Lotka-Volterra model with demographic noise
Phys Rev Lett, 126 (25) (2021), Article 258301
View in Scopus Google Scholar
[21]
Bunin G.
Ecological communities with Lotka-Volterra dynamics
Phys Rev E, 95 (2017), Article 042414
View at publisherCrossref View in Scopus Google Scholar
[22]
Pham T.M., Kaneko K.
Double-replica theory for evolution of genotype-phenotype interrelationship
Phys Rev Res, 5 (2) (2023), Article 023049
View in Scopus Google Scholar
[23]
Correale L., Leone M., Pagnani A., Weigt M., Zecchina R.
The computational core and fixed point organization in Boolean networks
J Stat Mech Theory Exp, 2006 (03) (2006), p. P03002
View in Scopus Google Scholar
[24]
Fyodorov Y.V., Khoruzhenko B.A.
Nonlinear analogue of the May-Wigner instability transition
Proc Natl Acad Sci, 113 (25) (2016), pp. 6827-6832
View at publisherCrossref View in Scopus Google Scholar
[25]
Hwang S., Folli V., Lanza E., Parisi G., Ruocco G., Zamponi F.
On the number of limit cycles in asymmetric neural networks
J Stat Mech Theory Exp, 2019 (5) (2019), Article 053402
View at publisher
Crossref Google Scholar
[26]
Fedeli S.B., Fyodorov Y.V., Ipsen J.
Nonlinearity-generated resilience in large complex systems
Phys Rev E, 103 (2) (2021), Article 022201
Google Scholar
[27]
Lacroix-A-Chez-Toine B., Fyodorov Y.V.
Counting equilibria in a random non-gradient dynamics with heterogeneous relaxation rates
J Phys A, 55 (14) (2022), Article 144001
View at publisher
Crossref View in Scopus Google Scholar
[28]
Ros V., Roy F., Biroli G., Bunin G., Turner A.M.
Generalized Lotka-Volterra equations with random, nonreciprocal interactions: The typical number of equilibria
Phys Rev Lett, 130 (2023), Article 257401
View in Scopus Google Scholar
[29]
Sompolinsky H., Crisanti A., Sommers H.-J.
Chaos in random neural networks
Phys Rev Lett, 61 (3) (1988), p. 259
View in Scopus Google Scholar
[30]
Ginzburg I., Sompolinsky H.
Theory of correlations in stochastic neural networks
Phys Rev E, 50 (1994), pp. 3171-3191
View in Scopus Google Scholar
[31]
Fruchart M., Hanai R., Littlewood P.B., Vitelli V.
Non-reciprocal phase transitions
Nature, 592 (7854) (2021), pp. 363-369
View at publisher
Crossref View in Scopus Google Scholar
[32]
Wainrib G., Touboul J.
Topological and dynamical complexity of random neural networks
Phys Rev Lett, 110 (2013), Article 118101
View at publisherCrossref View in Scopus Google Scholar
[33]
Dayan P., Abbott L.F.
Theoretical neuroscience: computational and mathematical modeling of neural systems
Cambridge, MIT Press, USA (2006)
Google Scholar
[34]
Van Mieghem P.
Graph spectra for complex networks
Cambridge University Press (2010)
Google Scholar
[35]
Dorogovtsev S.N., Goltsev A.V., Mendes J.F., Samukhin A.N.
Spectra of complex networks
Phys Rev E, 68 (4) (2003), Article 046109
Google Scholar
[36]
Tao T.
Outliers in the spectrum of iid matrices with bounded rank perturbations
Probab Theory Related Fields, 155 (1) (2013), pp. 231-263
Finding PDF…Crossref View in Scopus Google Scholar
[37]
Rajan K., Abbott L.
Eigenvalue spectra of random matrices for neural networks
Phys Rev Lett, 97 (18) (2006), Article 188104
Finding PDF…Crossref View in Scopus Google Scholar
[38]
Knowles A., Yin J.
The isotropic semicircle law and deformation of wigner matrices
Comm Pure Appl Math, 66 (11) (2013), pp. 1663-1749
Finding PDF…Crossref View in Scopus Google Scholar
[39]
Knowles A., Yin J.
The outliers of a deformed Wigner matrix
Ann Probab, 42 (5) (2014)
Google Scholar
[40]
Sommers H.J., Crisanti A., Sompolinsky H., Stein Y.
Spectrum of large random asymmetric matrices
Phys Rev Lett, 60 (19) (1988), pp. 1895-1898
View in Scopus Google Scholar
[41]
O’Rourke S., Renfrew D.
Low rank perturbations of large elliptic random matrices
Electron J Probab, 19 (43) (2014), pp. 1-65
arXiv:1309.5326 [math-ph]
View in Scopus Google Scholar
[42]
Dahmen D., Grün S., Diesmann M., Helias M.
Second type of criticality in the brain uncovers rich multiple-neuron dynamics
Proc Natl Acad Sci, 116 (26) (2019), pp. 13051-13060
Finding PDF…Crossref View in Scopus Google Scholar
[43]
Garcia E.R., Crumpton M.J., Galla T.
Niche overlap and Hopfield-like interactions in generalized random Lotka-Volterra systems
Phys Rev E, 108 (2023), Article 034120
View in Scopus Google Scholar
[44]
Neubert M.G., Caswell H.
Alternatives to resilience for measuring the responses of ecological systems to perturbations
Ecology, 78 (3) (1997), pp. 653-665
View in Scopus Google Scholar
[45]
Tang S., Allesina S.
Reactivity and stability of large ecosystems
Front Ecol Evol, 2 (2014), p. 21
View in Scopus Google Scholar
[46]
Yan G., Martinez N.D., Liu Y.-Y.
Degree heterogeneity and stability of ecological networks
J R Soc Interface, 14 (131) (2017), Article 20170189
Finding PDF…Crossref View in Scopus Google Scholar
[47]
Caravelli F., Staniczenko P.P.
Bounds on transient instability for complex ecosystems
PLoS One, 11 (6) (2016), Article e0157876
View at publisherCrossref View in Scopus Google Scholar
[48]
Murphy B.K., Miller K.D.
Balanced amplification: a new mechanism of selective amplification of neural activity patterns
Neuron, 61 (4) (2009), pp. 635-648
View PDF View article View in Scopus Google Scholar
[49]
Di Santo S., Villegas P., Burioni R., Muñoz M.A.
Non-normality, reactivity, and intrinsic stochasticity in neural dynamics: a non-equilibrium potential approach
J Stat Mech Theory Exp, 2018 (7) (2018), Article 073402
View at publisher
Crossref View in Scopus Google Scholar
[50]
Schuecker J., Goedeke S., Helias M.
Optimal sequence memory in driven random networks
Phys Rev X, 8 (4) (2018), Article 041029
View in Scopus Google Scholar
[51]
Kadmon J., Sompolinsky H.
Transition to chaos in random neuronal networks
Phys Rev X, 5 (4) (2015), Article 041030
View at publisherCrossref View in Scopus Google Scholar
[52]
Rajan K., Abbott L.F., Sompolinsky H.
Stimulus-dependent suppression of chaos in recurrent neural networks
Phys Rev E, 82 (1) (2010), Article 011903
View at publisherCrossref View in Scopus Google Scholar
[53]
Mastrogiuseppe F., Ostojic S.
Intrinsically-generated fluctuating activity in excitatory-inhibitory networks
PLoS Comput Biol, 13 (4) (2017), Article e1005498
View at publisherCrossref View in Scopus Google Scholar
[54]
Martin P.C., Siggia E.D., Rose H.A.
Statistical dynamics of classical systems
Phys Rev A, 8 (1973), pp. 423-437
View in Scopus Google Scholar
[55]
Coolen A.
Chapter 15 statistical mechanics of recurrent neural networks II — dynamics
Moss F., Gielen S. (Eds.), Neuro-informatics and neural modelling, Handbook of biological physics, vol. 4, North-Holland (2001), pp. 619-684
View PDF View article View in Scopus Google Scholar
[56]
Hatchett J.P.L., Coolen A.C.C.
Asymmetrically extremely dilute neural networks with Langevin dynamics and unconventional results
J Phys A: Math Gen, 37 (29) (2004), p. 7199
View in Scopus Google Scholar
[57]
Crisanti A., Sompolinsky H.
Path integral approach to random neural networks
Phys Rev E, 98 (6) (2018), Article 062120
View in Scopus Google Scholar
[58]
Helias M., Dahmen D.
Dynamic mean-field theory for random networks
Statistical field theory for neural networks, Springer (2020), pp. 95-126
View at publisherCrossref View in Scopus Google Scholar
[59]
Opper M., Diederich S.
Phase transition and 1/f noise in a game dynamical model
Phys Rev Lett, 69 (10) (1992), pp. 1616-1619
View in Scopus Google Scholar
[60]
Sidhom L., Galla T.
Ecological communities from random generalized Lotka-Volterra dynamics with nonlinear feedback
Phys Rev E, 101 (3) (2020), Article 032101
View in Scopus Google Scholar
[61]
Sherrington D., Kirkpatrick S.
Solvable model of a spin-glass
Phys Rev Lett, 35 (26) (1975), pp. 1792-1796
View in Scopus Google Scholar
[62]
Cure S., Neri I.
Antagonistic interactions can stabilise fixed points in heterogeneous linear dynamical systems
SciPost Phys, 14 (5) (2023), p. 093
Google Scholar
[63]
Hanel R., Pöchacker M., Thurner S.
Living on the edge of chaos: minimally nonlinear models of genetic regulatory dynamics
Phil Trans R Soc A, 368 (2010), pp. 5583-5596
View at publisher
Crossref View in Scopus Google Scholar
[64]
Morales G.B., Di Santo S., Muñoz M.A.
Quasiuniversal scaling in mouse-brain neuronal activity stems from edge-of-instability critical dynamics
Proc Natl Acad Sci, 120 (9) (2023), Article e2208998120
View in Scopus Google Scholar
[65]
Morales G.B., Muñoz M.A.
Optimal input representation in neural systems at the edge of chaos
Biology, 10 (8) (2021), p. 702
View at publisherCrossref View in Scopus Google Scholar
[66]
Laffargue T., Lam K.-D.N.T., Kurchan J., Tailleur J.
Large deviations of Lyapunov exponents
J Phys A, 46 (25) (2013), Article 254002
View at publisher
Crossref View in Scopus Google Scholar
[67]
Laffargue T., Sollich P., Tailleur J., van Wijland F.
Large-scale fluctuations of the largest Lyapunov exponent in diffusive systems
Europhys Lett, 110 (1) (2015), p. 10006
View at publisher
Crossref View in Scopus Google Scholar
[68]
Wardak A., Gong P.
Extended Anderson criticality in heavy-tailed neural networks
Phys Rev Lett, 129 (2022), Article 048103
View in Scopus Google Scholar
[69]
Torrisi G., Kühn R., Annibale A.
Percolation on the gene regulatory network
J Stat Mech Theory Exp, 2020 (8) (2020), Article 083501
Crossref View in Scopus Google Scholar
[70]
Toral R., Colet P.
Stochastic numerical methods: an introduction for students and scientists
John Wiley & Sons (2014)
Google Scholar
[71]
Eckmann J.P., Ruelle D.
Ergodic theory of chaos and strange attractors
Rev Modern Phys, 57 (3) (1985), pp. 617-656
View in Scopus Google Scholar
[72]
Pikovsky A., Politi A.
Lyapunov exponents: a tool to explore complex dynamics
Cambridge University Press, Cambridge (2016)
Google Scholar
[73]
Benettin G., Galgani L., Giorgilli A., Strelcyn J.-M.
Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. Part 2: Numerical application
Meccanica, 15 (1) (1980), pp. 21-30
Google Scholar
[74]
Wolf A., Swift J.B., Swinney H.L., Vastano J.A.
Determining Lyapunov exponents from a time series
Physica D, 16 (3) (1985), pp. 285-317
View PDF View article View in Scopus Google Scholar

Cited by (3)

Dynamical Mean-Field Theory of Complex Systems on Sparse Directed Networks
2025, Physical Review Letters
Dynamical theory for adaptive systems
2024, Journal of Statistical Mechanics: Theory and Experiment
Frequency-Dependent Covariance Reveals Critical Spatiotemporal Patterns of Synchronized Activity in the Human Brain
2024, Physical Review Letters

¹: We use the “hat” notation for quantities measured in simulations, to distinguish them from their theoretical counterparts, as obtained from the path integral approach. In particular, single-realizations are denoted by a subscript $s$ , e.g., ${\hat{M}}_{s}$ , whereas averages over many realizations are denoted with $\hat{M}$ .

²: For an ergodic system, where time averages are equivalent to ensemble averages, such quantity should equate the ensemble-averaged correlator $C (τ)$ , in the thermodynamic limit, if the system is at stationarity.

View Abstract

[1] [1]
Mézard M., Parisi G., Virasoro M.A.
Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications
World Scientific Publishing Company (1987)
Google Scholar

[2] [2]
Castellani T., Cavagna A.
Spin-glass theory for pedestrians
J Stat Mech Theory Exp, 2005 (05) (2005), p. P05012
View in Scopus Google Scholar

[3] [3]
Nishimori H.
Statistical physics of spin glasses and information processing: an introduction
Clarendon Press (2001)
Google Scholar

[4] [4]
Charbonneau P., Marinari E., Parisi G., Ricci-tersenghi F., Sicuro G., Zamponi F., et al.
Spin glass theory and far beyond: Replica symmetry breaking after 40 years
World Scientific (2023)
Google Scholar

[5] [5]
Aspelmeier T., Bray A.J., Moore M.A.
Complexity of ising spin glasses
Phys Rev Lett, 92 (2004), Article 087203
Google Scholar

[6] [6]
Cavagna A., Giardina I., Parisi G.
Numerical study of metastable states in ising spin glasses
Phys Rev Lett, 92 (2004), Article 120603
Google Scholar

[7] [7]
Annibale A., Gualdi G., Cavagna A.
Coexistence of supersymmetric and supersymmetry-breaking states in spherical spin-glasses
J Phys A: Math Gen, 37 (47) (2004), p. 11311
View at publisher
Crossref View in Scopus Google Scholar

[8] [8]
Müller M., Wyart M.
Marginal stability in structural, spin and electron glasses
Annu Rev Condens Matter Phys, 6 (1) (2015), pp. 177-200
View at publisherCrossref Google Scholar

[9] [9]
Herz A.V., Hopfield J.J.
Earthquake cycles and neural reverberations: Collective oscillations in systems with pulse-coupled threshold elements
Phys Rev Lett, 75 (1222) (1995), pp. 1222-1225
View in Scopus Google Scholar

[10] [10]
Challet D., Chessa A., Marsili M., Zhang Y.-C.
From minority games to real markets
Quant Finance, 1 (1) (2001), pp. 168-176
View in Scopus Google Scholar

[11] [11]
Stokić D., Hanel R., Thurner S.
Inflation of the edge of chaos in a simple model of gene interaction networks
Phys Rev E, 77 (2008), Article 061917
View at publisherCrossref View in Scopus Google Scholar

[12] [12]
Beggs J.M., Timme N.
Being critical of criticality in the brain
Front Physiol, 3 (163) (2012)
Google Scholar

[13] [13]
Mora T., Bialek W.
Are biological systems poised at criticality?
J Stat Phys, 144 (268) (2011)
Google Scholar

[14] [14]
Muñoz M.A.
Colloquium: Criticality and dynamical scaling in living systems
Rev Modern Phys, 90 (3) (2018), Article 031001
View in Scopus Google Scholar

[15] [15]
Moretti P., Muñoz M.A.
Griffiths phases and the stretching of criticality in brain networks
Nat Commun, 4 (1) (2013), pp. 1-10
Google Scholar

[16] [16]
Martinello M., Hidalgo J., Maritan A., Di Santo S., Plenz D., Muñoz M.A.
Neutral theory and scale-free neural dynamics
Phys Rev X, 7 (4) (2017), Article 041071
View in Scopus Google Scholar

[17] [17]
di Santo S., Villegas P., Burioni R., Muñoz M.A.
Landau–Ginzburg theory of cortex dynamics: Scale-free avalanches emerge at the edge of synchronization
Proc Natl Acad Sci, 115 (7) (2018), pp. E1356-E1365
View in Scopus Google Scholar

[18] [18]
Biroli G., Bunin G., Cammarota C.
Marginally stable equilibria in critical ecosystems
New J Phys, 20 (8) (2018), Article 083051
View at publisher
Crossref View in Scopus Google Scholar

[19] [19]
Altieri A., Franz S.
Constraint satisfaction mechanisms for marginal stability and criticality in large ecosystems
Phys Rev E, 99 (1) (2019), Article 010401
View in Scopus Google Scholar

[20] [20]
Altieri A., Roy F., Cammarota C., Biroli G.
Properties of equilibria and glassy phases of the random Lotka-Volterra model with demographic noise
Phys Rev Lett, 126 (25) (2021), Article 258301
View in Scopus Google Scholar

[21] [21]
Bunin G.
Ecological communities with Lotka-Volterra dynamics
Phys Rev E, 95 (2017), Article 042414
View at publisherCrossref View in Scopus Google Scholar

[22] [22]
Pham T.M., Kaneko K.
Double-replica theory for evolution of genotype-phenotype interrelationship
Phys Rev Res, 5 (2) (2023), Article 023049
View in Scopus Google Scholar

[23] [23]
Correale L., Leone M., Pagnani A., Weigt M., Zecchina R.
The computational core and fixed point organization in Boolean networks
J Stat Mech Theory Exp, 2006 (03) (2006), p. P03002
View in Scopus Google Scholar

[24] [24]
Fyodorov Y.V., Khoruzhenko B.A.
Nonlinear analogue of the May-Wigner instability transition
Proc Natl Acad Sci, 113 (25) (2016), pp. 6827-6832
View at publisherCrossref View in Scopus Google Scholar

[25] [25]
Hwang S., Folli V., Lanza E., Parisi G., Ruocco G., Zamponi F.
On the number of limit cycles in asymmetric neural networks
J Stat Mech Theory Exp, 2019 (5) (2019), Article 053402
View at publisher
Crossref Google Scholar

[26] [26]
Fedeli S.B., Fyodorov Y.V., Ipsen J.
Nonlinearity-generated resilience in large complex systems
Phys Rev E, 103 (2) (2021), Article 022201
Google Scholar

[27] [27]
Lacroix-A-Chez-Toine B., Fyodorov Y.V.
Counting equilibria in a random non-gradient dynamics with heterogeneous relaxation rates
J Phys A, 55 (14) (2022), Article 144001
View at publisher
Crossref View in Scopus Google Scholar

[28] [28]
Ros V., Roy F., Biroli G., Bunin G., Turner A.M.
Generalized Lotka-Volterra equations with random, nonreciprocal interactions: The typical number of equilibria
Phys Rev Lett, 130 (2023), Article 257401
View in Scopus Google Scholar

[29] [29]
Sompolinsky H., Crisanti A., Sommers H.-J.
Chaos in random neural networks
Phys Rev Lett, 61 (3) (1988), p. 259
View in Scopus Google Scholar

[30] [30]
Ginzburg I., Sompolinsky H.
Theory of correlations in stochastic neural networks
Phys Rev E, 50 (1994), pp. 3171-3191
View in Scopus Google Scholar

[31] [31]
Fruchart M., Hanai R., Littlewood P.B., Vitelli V.
Non-reciprocal phase transitions
Nature, 592 (7854) (2021), pp. 363-369
View at publisher
Crossref View in Scopus Google Scholar

[32] [32]
Wainrib G., Touboul J.
Topological and dynamical complexity of random neural networks
Phys Rev Lett, 110 (2013), Article 118101
View at publisherCrossref View in Scopus Google Scholar

[33] [33]
Dayan P., Abbott L.F.
Theoretical neuroscience: computational and mathematical modeling of neural systems
Cambridge, MIT Press, USA (2006)
Google Scholar

[34] [34]
Van Mieghem P.
Graph spectra for complex networks
Cambridge University Press (2010)
Google Scholar

[35] [35]
Dorogovtsev S.N., Goltsev A.V., Mendes J.F., Samukhin A.N.
Spectra of complex networks
Phys Rev E, 68 (4) (2003), Article 046109
Google Scholar

[36] [36]
Tao T.
Outliers in the spectrum of iid matrices with bounded rank perturbations
Probab Theory Related Fields, 155 (1) (2013), pp. 231-263
Finding PDF…Crossref View in Scopus Google Scholar

[37] [37]
Rajan K., Abbott L.
Eigenvalue spectra of random matrices for neural networks
Phys Rev Lett, 97 (18) (2006), Article 188104
Finding PDF…Crossref View in Scopus Google Scholar

[38] [38]
Knowles A., Yin J.
The isotropic semicircle law and deformation of wigner matrices
Comm Pure Appl Math, 66 (11) (2013), pp. 1663-1749
Finding PDF…Crossref View in Scopus Google Scholar

[39] [39]
Knowles A., Yin J.
The outliers of a deformed Wigner matrix
Ann Probab, 42 (5) (2014)
Google Scholar

[40] [40]
Sommers H.J., Crisanti A., Sompolinsky H., Stein Y.
Spectrum of large random asymmetric matrices
Phys Rev Lett, 60 (19) (1988), pp. 1895-1898
View in Scopus Google Scholar

[41] [41]
O’Rourke S., Renfrew D.
Low rank perturbations of large elliptic random matrices
Electron J Probab, 19 (43) (2014), pp. 1-65
arXiv:1309.5326 [math-ph]
View in Scopus Google Scholar

[42] [42]
Dahmen D., Grün S., Diesmann M., Helias M.
Second type of criticality in the brain uncovers rich multiple-neuron dynamics
Proc Natl Acad Sci, 116 (26) (2019), pp. 13051-13060
Finding PDF…Crossref View in Scopus Google Scholar

[43] [43]
Garcia E.R., Crumpton M.J., Galla T.
Niche overlap and Hopfield-like interactions in generalized random Lotka-Volterra systems
Phys Rev E, 108 (2023), Article 034120
View in Scopus Google Scholar

[44] [44]
Neubert M.G., Caswell H.
Alternatives to resilience for measuring the responses of ecological systems to perturbations
Ecology, 78 (3) (1997), pp. 653-665
View in Scopus Google Scholar

[45] [45]
Tang S., Allesina S.
Reactivity and stability of large ecosystems
Front Ecol Evol, 2 (2014), p. 21
View in Scopus Google Scholar

[46] [46]
Yan G., Martinez N.D., Liu Y.-Y.
Degree heterogeneity and stability of ecological networks
J R Soc Interface, 14 (131) (2017), Article 20170189
Finding PDF…Crossref View in Scopus Google Scholar

[47] [47]
Caravelli F., Staniczenko P.P.
Bounds on transient instability for complex ecosystems
PLoS One, 11 (6) (2016), Article e0157876
View at publisherCrossref View in Scopus Google Scholar

[48] [48]
Murphy B.K., Miller K.D.
Balanced amplification: a new mechanism of selective amplification of neural activity patterns
Neuron, 61 (4) (2009), pp. 635-648
View PDF View article View in Scopus Google Scholar

[49] [49]
Di Santo S., Villegas P., Burioni R., Muñoz M.A.
Non-normality, reactivity, and intrinsic stochasticity in neural dynamics: a non-equilibrium potential approach
J Stat Mech Theory Exp, 2018 (7) (2018), Article 073402
View at publisher
Crossref View in Scopus Google Scholar

[50] [50]
Schuecker J., Goedeke S., Helias M.
Optimal sequence memory in driven random networks
Phys Rev X, 8 (4) (2018), Article 041029
View in Scopus Google Scholar

[51] [51]
Kadmon J., Sompolinsky H.
Transition to chaos in random neuronal networks
Phys Rev X, 5 (4) (2015), Article 041030
View at publisherCrossref View in Scopus Google Scholar

[52] [52]
Rajan K., Abbott L.F., Sompolinsky H.
Stimulus-dependent suppression of chaos in recurrent neural networks
Phys Rev E, 82 (1) (2010), Article 011903
View at publisherCrossref View in Scopus Google Scholar

[53] [53]
Mastrogiuseppe F., Ostojic S.
Intrinsically-generated fluctuating activity in excitatory-inhibitory networks
PLoS Comput Biol, 13 (4) (2017), Article e1005498
View at publisherCrossref View in Scopus Google Scholar

[54] [54]
Martin P.C., Siggia E.D., Rose H.A.
Statistical dynamics of classical systems
Phys Rev A, 8 (1973), pp. 423-437
View in Scopus Google Scholar

[55] [55]
Coolen A.
Chapter 15 statistical mechanics of recurrent neural networks II — dynamics
Moss F., Gielen S. (Eds.), Neuro-informatics and neural modelling, Handbook of biological physics, vol. 4, North-Holland (2001), pp. 619-684
View PDF View article View in Scopus Google Scholar

[56] [56]
Hatchett J.P.L., Coolen A.C.C.
Asymmetrically extremely dilute neural networks with Langevin dynamics and unconventional results
J Phys A: Math Gen, 37 (29) (2004), p. 7199
View in Scopus Google Scholar

[57] [57]
Crisanti A., Sompolinsky H.
Path integral approach to random neural networks
Phys Rev E, 98 (6) (2018), Article 062120
View in Scopus Google Scholar

[58] [58]
Helias M., Dahmen D.
Dynamic mean-field theory for random networks
Statistical field theory for neural networks, Springer (2020), pp. 95-126
View at publisherCrossref View in Scopus Google Scholar

[59] [59]
Opper M., Diederich S.
Phase transition and 1/f noise in a game dynamical model
Phys Rev Lett, 69 (10) (1992), pp. 1616-1619
View in Scopus Google Scholar

[60] [60]
Sidhom L., Galla T.
Ecological communities from random generalized Lotka-Volterra dynamics with nonlinear feedback
Phys Rev E, 101 (3) (2020), Article 032101
View in Scopus Google Scholar

[61] [61]
Sherrington D., Kirkpatrick S.
Solvable model of a spin-glass
Phys Rev Lett, 35 (26) (1975), pp. 1792-1796
View in Scopus Google Scholar

[62] [62]
Cure S., Neri I.
Antagonistic interactions can stabilise fixed points in heterogeneous linear dynamical systems
SciPost Phys, 14 (5) (2023), p. 093
Google Scholar

[63] [63]
Hanel R., Pöchacker M., Thurner S.
Living on the edge of chaos: minimally nonlinear models of genetic regulatory dynamics
Phil Trans R Soc A, 368 (2010), pp. 5583-5596
View at publisher
Crossref View in Scopus Google Scholar

[64] [64]
Morales G.B., Di Santo S., Muñoz M.A.
Quasiuniversal scaling in mouse-brain neuronal activity stems from edge-of-instability critical dynamics
Proc Natl Acad Sci, 120 (9) (2023), Article e2208998120
View in Scopus Google Scholar

[65] [65]
Morales G.B., Muñoz M.A.
Optimal input representation in neural systems at the edge of chaos
Biology, 10 (8) (2021), p. 702
View at publisherCrossref View in Scopus Google Scholar

[66] [66]
Laffargue T., Lam K.-D.N.T., Kurchan J., Tailleur J.
Large deviations of Lyapunov exponents
J Phys A, 46 (25) (2013), Article 254002
View at publisher
Crossref View in Scopus Google Scholar

[67] [67]
Laffargue T., Sollich P., Tailleur J., van Wijland F.
Large-scale fluctuations of the largest Lyapunov exponent in diffusive systems
Europhys Lett, 110 (1) (2015), p. 10006
View at publisher
Crossref View in Scopus Google Scholar

[68] [68]
Wardak A., Gong P.
Extended Anderson criticality in heavy-tailed neural networks
Phys Rev Lett, 129 (2022), Article 048103
View in Scopus Google Scholar

[69] [69]
Torrisi G., Kühn R., Annibale A.
Percolation on the gene regulatory network
J Stat Mech Theory Exp, 2020 (8) (2020), Article 083501
Crossref View in Scopus Google Scholar

[70] [70]
Toral R., Colet P.
Stochastic numerical methods: an introduction for students and scientists
John Wiley & Sons (2014)
Google Scholar

[71] [71]
Eckmann J.P., Ruelle D.
Ergodic theory of chaos and strange attractors
Rev Modern Phys, 57 (3) (1985), pp. 617-656
View in Scopus Google Scholar

[72] [72]
Pikovsky A., Politi A.
Lyapunov exponents: a tool to explore complex dynamics
Cambridge University Press, Cambridge (2016)
Google Scholar

[73] [73]
Benettin G., Galgani L., Giorgilli A., Strelcyn J.-M.
Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. Part 2: Numerical application
Meccanica, 15 (1) (1980), pp. 21-30
Google Scholar

[74] [74]
Wolf A., Swift J.B., Swinney H.L., Vastano J.A.
Determining Lyapunov exponents from a time series
Physica D, 16 (3) (1985), pp. 285-317
View PDF View article View in Scopus Google Scholar

Outline

Cited by (3)

Figures (14)

Extras (1)

Chaos, Solitons & Fractals

Highlights

Abstract

1. Introduction

2. Rate model

3. Linear-stability analysis

4. Path-integral formalism

4.1. Uncorrelated interactions ( $γ = 0$ ) and noiseless dynamics ( $σ = 0$ )

4.1.1. Fixed-point solutions

4.1.2. Non-fixed point steady states

4.1.2.1. For $M \neq 0$

4.1.2.2. For $M = 0$

4.2. Switching on fluctuations: Noisy dynamics ( $σ \neq 0$ ) with uncorrelated interactions ( $γ = 0$ )

4.3. Switching on correlations ( $γ \neq 0$ )

5. Discussion and conclusions

CRediT authorship contribution statement

Declaration of competing interest

Acknowledgements

Appendix A. On the equivalence between our model and the classic one (SCS)

Appendix B. Stability of steady-state solutions with uncorrelated interactions ( $γ = 0$ )

Appendix C. Phase diagram for $γ \neq 0$

Appendix D. Numerical procedure

D.1. Computation of trajectories

D.2. Shape of the potential $V (C | C (0), M)$

D.3. Computation of the largest Lyapunov exponent

Appendix E. Supplementary data

Data availability

References

Cited by (3)

Dynamical Mean-Field Theory of Complex Systems on Sparse Directed Networks

Dynamical theory for adaptive systems

Frequency-Dependent Covariance Reveals Critical Spatiotemporal Patterns of Synchronized Activity in the Human Brain

Dynamically selected steady states and criticality in non-reciprocal networks

Highlights

Abstract

1. Introduction

2. Rate model

3. Linear-stability analysis

4. Path-integral formalism

4.1. Uncorrelated interactions (γ=0) and noiseless dynamics (σ=0)

4.1.1. Fixed-point solutions

4.1.2. Non-fixed point steady states

4.1.2.1. For M≠0

4.1.2.2. For M=0

4.2. Switching on fluctuations: Noisy dynamics (σ≠0) with uncorrelated interactions (γ=0)

4.3. Switching on correlations (γ≠0)

5. Discussion and conclusions

CRediT authorship contribution statement

Declaration of competing interest

Acknowledgements

Appendix A. On the equivalence between our model and the classic one (SCS)

Appendix B. Stability of steady-state solutions with uncorrelated interactions (γ=0)

Appendix C. Phase diagram for γ≠0

Appendix D. Numerical procedure

D.1. Computation of trajectories

D.2. Shape of the potential V(C|C(0),M)

D.3. Computation of the largest Lyapunov exponent

Appendix E. Supplementary data

Data availability

References

4.1. Uncorrelated interactions ( $γ = 0$ ) and noiseless dynamics ( $σ = 0$ )

4.1.2.1. For $M \neq 0$

4.1.2.2. For $M = 0$

4.2. Switching on fluctuations: Noisy dynamics ( $σ \neq 0$ ) with uncorrelated interactions ( $γ = 0$ )

4.3. Switching on correlations ( $γ \neq 0$ )

Appendix B. Stability of steady-state solutions with uncorrelated interactions ( $γ = 0$ )

Appendix C. Phase diagram for $γ \neq 0$

D.2. Shape of the potential $V (C | C (0), M)$