Correlations between synapses in pairs of neurons slow down dynamics in randomly connected neural networks
Daniel Martí Nicolas Brunel Srdjan Ostojic
q-bio.NC





Laboratoire de Neurosciences Cognitives,

Inserm UMR No. 960, Ecole Normale Supérieure,

PSL Research University, Paris, France

Department of Statistics and Department of Neurobiology,

University of Chicago, Chicago, IL 60637, USA
Department of Neurobiology and Department of Physics,

Duke University, Durham, NC 27710, USA

Laboratoire de Neurosciences Cognitives,

Inserm UMR No. 960, Ecole Normale Supérieure,

PSL Research University, Paris, France




 Networks of randomly connected neurons are among the most popular models in
 theoretical neuroscience. The connectivity between neurons
 in the cortex is however not fully random, the simplest and
 most prominent deviation from randomness found in experimental data
 being the overrepresentation of bidirectional connections among
 pyramidal cells. Using numerical and analytical methods, we
 investigated the effects of partially symmetric connectivity on
 the dynamics in networks of rate units. We consider the two dynamical
 regimes exhibited by random neural networks: the weak-coupling
 regime, where the firing activity decays to a single fixed point
 unless the network is stimulated, and the strong-coupling or chaotic
 regime, characterized by internally generated fluctuating firing
 rates. In the weak-coupling regime, we compute analytically, for an
 arbitrary degree of symmetry, the autocorrelation of network
 activity in the presence of external noise. In the chaotic regime, we
 perform simulations to determine the timescale of the intrinsic
 fluctuations. In both cases, symmetry increases the characteristic
 asymptotic decay time of the autocorrelation function and therefore
 slows down the dynamics in the network.


Correlations between synapses in pairs of neurons slow down dynamics in randomly connected neural networks
    Srdjan Ostojic
    December 30, 2023
==========================================================================================================



 

§ INTRODUCTION


 The dynamics and function of a network of neurons is to a large
 extent determined by its pattern of synaptic connections. In the
 mammalian brain, cortical networks exhibit a complex connectivity
 that to a first approximation can be regarded as random. This
 connectivity structure has motivated the study of networks of
 neurons connected through a random synaptic weight matrix with
 independent and identically distributed (i.i.d.) entries, which have become a
 central paradigm in theoretical
 neuroscience <cit.>.
 Randomly connected networks of firing-rate units exhibit a chaotic phase <cit.>, which can be exploited as a
 susbstrate for complex
 computations <cit.>.
 Networks of randomly connected spiking neurons also exhibit rich dynamics that can account for
 the highly irregular spontaneous activity observed in
 the cortex in
 vivo <cit.>.
 Importantly, these models are to a large extent amenable to a mathematical
 analysis, which allows for a thorough understanding of the
 mechanisms underlying their dynamics.

 Detailed analyses of experimental data on cortical connectivity have however
 identified patterns of connectivity that strongly deviate from the
 i.i.d. assumption <cit.>. The
 most prominent of such deviations is the overrepresentation of
 reciprocal
 connections <cit.>,
 and the fact that synapses of bidirectionally connected pairs of
 neurons are on average stronger than synapses of unidirectionally
 connected pairs. These observations are consistent with a partially
 symmetric connectivity structure, intermediate between full symmetry
 and full asymmetry. How partial symmetry in the connectivity impacts
 network dynamics is not yet understood, in part because such
 partial symmetry renders the mathematical analyses more
 challenging <cit.>. Here we study the
 impact of partial symmetry in the connectivity structure on the dynamics
 of a simple network model consisting of interacting rate
 units. Depending on the overall strength of coupling, such a network
 can display either a stable or a chaotic regime of activity, as in
 the random asymmetric case <cit.>. We examined how the
 degree of symmetry in the network influences the temporal dynamics
 in both regimes. For the stable regime, we exploited recent results
 from random matrix theory <cit.> to
 derive analytical expressions for the autocorrelation
 functions. These expressions demonstrate that increasing the
 symmetry in the network leads to a slowing down of the
 dynamics. Numerical simulations in the chaotic regime show a similar
 effect, with time scales increasing far more substantially with
 symmetry than in the fixed point regime. Altogether, our results indicate that symmetry in the connectivity can act as an additionnal source of slow dynamics, an important ingredient for implementing computations in networks of neurons <cit.>.

 

§ DESCRIPTION OF THE MODEL


 We consider a network of N fully connected neurons, each described 
 by an activation variable (synaptic current) x_i, i=1,…,N, obeying
 
    x̣_i/ṭ = -x_i + g ∑_j=1^N
     J_ijϕ(x_j),

 where g is a gain parameter that modulates the strength of
 recurrent connections, and where ϕ(·) is the input-output
 transfer function that transforms activations x_i into firing
 rates. This transformation is non-linear and we model it as
 ϕ(x)=tanh(x) for mathematical convenience (see <cit.> for studies of network models with
 different choices of ϕ). The elements J_ij of the
 connectivity matrix are drawn from a Gaussian distribution with zero
 mean, variance 1 / N, and correlation
 
    [J_ij J_ji]_J = η / N,

 with the square brackets [·]_J denoting an average over
 realizations of the random connections. The parameter η is the
 correlation coefficient between the two weights connecting pairs of
 neurons, and quantifies the degree of symmetry of the
 connections. For η=0 the elements J_ij and J_ji are
 independent and the connectivity matrix is fully asymmetric; for
 η=1 the connectivity matrix is fully symmetric; For η=-1 it
 is fully antisymmetric.  In Secs. I–IIIA we study the full range
 η∈[-1,1], while in Secs IIIB–IV we focus on η∈[0,1].

 

§ DYNAMICAL REGIMES OF THE NETWORK


 For fully asymmetric matrices, previous work has shown that the
 network activity described by (<ref>) undergoes a phase
 transition at g=1 in the limit of large N <cit.>. For g < 1 the
 activity for all units decays to 0, which is the unique stable fixed point of the
 dynamics <cit.>, while for g > 1 the activity is chaotic. Such a
 transition can be partially understood by assessing the stability of
 the fixed point at x_i=0 for i=1, …, N. If we linearize
 Eq. (<ref>) around this fixed point we obtain the
 stability matrix, with components
 
    M_ij = -δ_ij + gJ_ij.

 The eigenvalues of M_ij are therefore those of the matrix J_ij, scaled
 by the gain g and shifted along the real axis by -1. In the limit
 N→∞, for a connectivity matrix J_ij whose entries are
 i.i.d. Gaussian random variables of zero mean and variance 1/N, eigenvalues
 are uniformly distributed in the unit disk of the complex
 plane <cit.>.
 This implies that the eigenvalues of the stability matrix have a negative real
 part as long as g < 1, and therefore that the fixed point at 0 is stable in
 that range.

 An analogous transition occurs when connections are partially
 symmetric. The presence of correlations among weights deforms the
 spectrum of eigenvalues into an ellipse, elongating its major radius
 by a factor of 1 + η and shortening the minor radius by a
 factor 1 - η <cit.>
 [Fig. <ref>a]. This property is usually referred to as the
 elliptic law. For the network described by (<ref>) such a
 deformation causes the fixed point at x_i = 0 for i=1,…,N
 to lose its stability at g = 1 / (1 + η)
 [Figs <ref>b and <ref>c]. In other words, symmetry lowers the
 critical coupling.



 Our goal is to characterize how the degree of symmetry in the
 connections affects the network activity on each side of the
 instability: the relaxation response of the network at low gains and
 the chaotic self-generated activity observed at strong
 gains. Our description of the network activity will be based on the
 average autocorrelation function,
 
    C(τ) = 1/N∑_i=1^N[x_i(t) x_i(t + τ)]_J,

 where the average is over both the population and the realizations of the
 connectivity matrix <cit.>, and where
 we are assuming for now that the system is stationary.
 





§ DYNAMICS IN THE FIXED-POINT REGIME




 §.§ Derivation of the autocorrelation function



In the fixed-point regime, the activity decays to zero unless the
network is stimulated by external inputs. To characterize the dynamics
of the network in this regime, we induce network activity by feeding
each neuron with independent Gaussian white noise <cit.>. The amplitude of
this noise is assumed to be small enough so that the synaptic
activation of all neurons lies within the linear range of their
input-to-rate transfer function [see Fig. <ref>f for the range of validity of that approximation]. Under these conditions, ϕ(x) can
be approximated by its first order Taylor expansion ϕ(0) +
ϕ'(x)|_x=0 x = x, and the dynamical equations become
 
    x̣̣⃗(̣ṭ)̣/ṭ = (-1 + g𝐉)x⃗(t) + σξ(t),

 where x⃗(t)=(x_1(t), …, x_N(t))^T, 1 is the
 identity matrix, J⃗ is the connectivity matrix, and
 ξ(t) = (ξ_1(t), …, ξ_N(t))^T is a vector
 of independent white noise sources of zero mean and unit variance:
 ξ_i(t) = 0, ξ_i(t) ξ_j(t') = δ_ijδ(t - t'), with angular brackets representing averages over
 noise realizations. The parameter σ is the standard deviation
 of the white noise injected into neurons.

 The time scales displayed by a linear system like (<ref>)
 are strongly affected by the real part of the eigenvalues of the
 system's stability matrix and, in particular, they get longer as
 eigenvalues get closer to the imaginary axis. To disentangle this
 type of slowing down from the effects due to symmetry alone, we vary
 the parameter η while keeping the spectral gap fixed. By
 spectral gap we mean the distance between the spectrum of
 eigenvalues of the stability matrix M_ij,
 Eq. (<ref>), and the imaginary axis [Fig. <ref>a]. From the elliptic law, the eigenvalue of the
 stability matrix with the largest real part is z = -1 + g(1 + η),
 and we can keep the spectral gap at δ by setting the gain to
 g = (1 - δ) / (1 + η).





The system described by (<ref>) is linear and can be solved
by diagonalizing the connectivity matrix. The matrix J⃗ admits
a set of right eigenvectors {R⃗_1, …, R⃗_N} that
obey J⃗R⃗_i = λ_i R⃗_i for
i=1,…,N. These eigenvectors are in general complex-valued and,
except for the symmetric case η=1, not orthogonal to one another,
which implies that J⃗ cannot be diagonalized through a unitary
transformation. Matrices of this kind are called non-normal and
do not commute with their transpose conjugate: J⃗J⃗^†≠J⃗^†J⃗ <cit.>. Even if
non-normal matrices cannot be diagonalized by an orthogonal set of
eigenvectors, it is always possible to form a biorthogonal
basis by extending the set of right eigenvectors with the set of left
eigenvectors, which obey L⃗_i^†J⃗ = λ_i
L⃗^†_i. This extended basis is biorthogonal in the sense
that L⃗^†_i R⃗_j = δ_ij. We can summarize all
these properties in a compact way by defining the square matrices
R⃗ and L⃗ that result from adjoining in columns the set
of, respectively, right and left eigenvectors, and by introducing the
diagonal matrix Λ that contains the eigenvalues
λ_i of J⃗ in its diagonal entries. In this notation the
biorthogonality condition is L⃗^†R⃗ = 1⃗ and
the eigenvalue equations for the right and left eigenvectors read
J⃗R⃗ = R⃗Λ and L⃗^†J⃗ =
ΛL⃗^†.

We can now write the formal solution of (<ref>):

    x⃗(t)     = σ∫_-∞^t^(-1⃗ + gJ⃗)(t -
    s)ξ(s) ṣ
        = σR⃗∫_-∞^t^(-1⃗ + g Λ)(t -
    s)R⃗^-1ξ(s) ṣ ,

where in the last equality we used the basis of right eigenvectors to
write J⃗ = R⃗ΛR⃗^-1 and we
implicitly expanded the exponential in its power series to obtain the
final result. From this expression we can derive the
population-average autocorrelation for a particular realization of the
connectivity:

    C_J(τ)     = 1/Nx⃗^†(t)
     x⃗(t + τ) =
     1/Nx⃗(t + τ)
     x⃗^†(t)
        = σ^2/N∫_0^∞^-2 u - τ{R⃗^†R⃗ ^gΛ(u + τ)L⃗^†L⃗ ^gΛ^†u} ụ.

In the second equality we used the cyclicity of the trace, and in the
last line we changed the integration variable to u=t - s and we used
the biorthogonality condition to write R⃗^-1 =
L⃗^†. The average over noise amounts to applying the
identity ⟨ξ(t)ξ^†(t')
⟩ = σ^2 1δ(t - t'). Note that the
σ^2 appears as an overall factor, so we can set σ=1
without loss of generality.

We can simplify (<ref>) by introducing the so-called overlap matrix, with components

    O_ij = (L⃗^†L⃗)_ij (R⃗^†R⃗)_ji,

and which characterizes the correlations between left and right eigenvectors <cit.>. Equation (<ref>) then becomes

    C_J(τ) = 1/N∫_0^∞^-2 u - τ∑_i=1^N∑_j=1^N^g λ_i(u + τ)
     O_ij^g λ̅_ju ụ.

If the connectivity matrix were normal, the overlap would be the
identity matrix and the autocorrelation (<ref>)
would just be a sum of independent contributions—one per
eigenvalue. These contributions are coupled for non-normal matrices.

We can make further analytical progress by studying the
autocorrelation (<ref>) in the limit N
→∞, in which the differences of C_J(τ) across
realizations of the connectivity matrix disappear, and the autocorrelations of all units become close to the population average [Fig. <ref>g]. In that limit sums
over indices are replaced with integrals over eigenvalues, while the
overlap matrix is replaced with the local average of the overlap,
defined as

    D(z_1, z_2) = lim_N→∞4/N[ ∑_i=1^N∑_j=1^N O_ijδ^2(z_1 - λ_i) δ^2(z_2 -
    λ_j)]_J.

Here z=x +  y are complex numbers and we defined the complex Dirac
delta as δ^2(z) ≡ (1/2) δ(x) δ(y) so that it
satisfies the normalization condition ∫δ^2(z) z≡∫δ^2(z) ẓẓ̣̅ = 1. After taking the
limit N→∞, Eq. (<ref>)
becomes

    C(τ) = ∫_0^∞^-2 u - τ A(u, τ) ụ.

where we defined

    A(u, τ) = 1/4∬^g(z_1 + z̅_2) u + g z_1 τ D(z_1,
     z_2) z_1z_2.

Each of the integrals in Eq. (<ref>) is over complex values,
and involves the expression of D(z_1, z_2) for the ensemble of
Gaussian random matrices with partial symmetry, which was derived
using diagrammatic techniques in <cit.> and
whose functional form can be found in Appendix <ref>. We
used the result of <cit.> to evaluate the
double complex integral A(u, τ) in Eq. (<ref>). The
details of the evaluation are given in Appendix <ref>, and
the result is

    A(u, τ) = A_1(u, τ) + A_2(u, τ),

with

    A_1(u, τ)     = (1 + η^2) I_0(g ψ(u, τ; η) )
        - 2 η(1 + 2(1 -η)^2τ^2/ψ(u, τ; η)^2)
     I_2(g  ψ(u, τ; η)), 
    
     A_2(u, τ)     = -1/g^2 u (u + τ)∑_k=1^∞η^k k^2 I_k(2g√(η)u)
       × I_k(2g√(η)(u + τ)),

where I_k(·) is the modified Bessel function of order k, and
where in Eq. (<ref>) we defined

    ψ(u, τ;η) = 2 √((1 + η)^2 u (u + τ) + ητ^2).

The autocorrelation is finally computed from Eq. (<ref>),
integrating numerically over u.

Expressions. (<ref>), (<ref>),
(<ref>) are valid for full range -1 ≤η≤ 1. For negative
η we replace √(η) with √(|η|) and apply the
identity I_ν( z) = ^ν J_ν(z), which is valid for integer
ν.

The analytical prediction given by Eqs. (<ref>) and
Eqs. (<ref>)–(<ref>) matches with the
autocorrelation estimated from numerical simulations
[Fig. <ref>b], although for long time lags the numerical
estimate becomes noisy due to finite-size effects. To check the
validity of our prediction also at long time lags, we compared our
analytical prediction with three alternative derivations
[Fig. <ref>c]. One such derivation consists of estimating
the autocorrelation for large but finite N, by computing numerically
the eigenvalues and eigenvectors of randomly generated matrices,
evaluating the time integral of Eq. (<ref>),
which gives

    C_J(τ) = -1/N∑_i=1^N∑_j=1^NO_ij^-(1 - gλ_i) τ/2 + g(λ_i + λ̅_j),

and then by averaging C_J(τ) over multiple realizations of the
connectivity matrix. Another derivation is based on dynamical
mean-field
theory <cit.>,
which gives rise to a set of integro-differential equations involving
C(τ) that can be solved numerically
(Appendix <ref>). Finally, we numerically computed the
inverse Fourier transform of the power spectrum derived in
<cit.> for this same system.
<cit.> used a perturbative method to
derive the system of integro-differential
equations (<ref>) and (<ref>), which they solved for
the correlation and reponse functions by using a Laplace transform.
All derivations yield the same result, except for the deviations we
observe when applying Eq. (<ref>) at long τ,
which are caused by finite-size effects.

Our results show that an increase in symmetry tends to spread
autocorrelations toward longer time lags, and that this effect gets
larger the closer the system gets to the onset of chaos [Fig. <ref>e]. An intuitive
explanation for this slowing down is that the deformation of the
eigenspectrum caused by symmetry increases the density of eigenvalues
with small imaginary parts, thereby enlarging the contribution of
low-frequency modes.




 §.§ Behavior at long time lags



While equations (<ref>)–(<ref>) are exact, they provide little analytical insight into how the autocorrelation depends on parameters. A more explicit dependence can be obtained by evaluating C(τ) in the limit of long τ. We relegate the details of the calculation to Appendix <ref> and summarize the main results here. The analysis shows that, in the fixed point regime, there exist two subregimes of activity that differ in how the asymptotic decay rate of the autocorrelation depends on the symmetry parameter η and the spectral gap δ. For small values of η and δ, the autocorrelation decays as a pure exponential at long τ (regime I),

    C(τ) = F_I(η, δ) ^-τ G_I(η, δ),

with

    F_I(η, δ)     = δ^-1/2(1 - η)^2/2√(2) (1 - δ),
    
     G_I(η, δ)     = 1 - η/1 + η√(2δ - δ^2).

Conversely, for sufficiently large values of η and δ autocorrelation for long τ can be approximated by a power multiplied by an exponential decay (regime II):

    C(τ) = τ^-3/2 F_II(η, δ) ^-τ G_II(η, δ),

with

    F_II(η, δ)     = 1/4√(π)( 1 + η/1 - δ)^3/2[  2 η^-1/4(1 + η^2)/δ(1 + η) - [1 - √(η)]^2aa
       - η^5/6 (1 + η)/(1 - √(η))^2 + 2 √(η)δ]
    
    G_II(η, δ)     = (1 - √(η))^2 + 2δ√(η)/1 + η.

A comparison between the asymptotic expression in
Eq. (<ref>) and the full expression for the
autocorrelation function reveals, however, that the power law is not
observed in practice because the range below the cutoff falls below
the values of τ where the asymptotic approximation starts
matching the exact expression.

Figure <ref>a shows the exact parameter region of each asymptotic regime, after transforming the spectral gaps into gains. In both regimes the autocorrelation's asymptotic decay rate matches the exact result for time lags longer than a few time units [see Fig. <ref>a, lateral panels]. It seems therefore reasonable to associate the time scale of the autocorrelation with the inverse of G_I,II(η,δ) [see Eqs. (<ref>) and (<ref>)], where the subindex I, II is chosen according to the subregime found at the parameter values (η, δ). The asymptotic time scale of the autocorrelation increases monotonically with symmetry regardless of the subregime the network operates in [Fig. <ref>b], although this dependence is convex in the exponential subregime and concave in the power-law-with-cutoff regime [in Fig. <ref>b see the curves split by the red dots, which mark the boundary between subregimes]. Note also that as the spectral gap δ shrinks to 0 the system enters the exponential regime and timescales diverge as δ^-1/2, according to Eq. (<ref>).



 §.§ Effect of overlaps


As shown in Eq. (<ref>), the autocorrelation function in general depends on two factors, the full eigenspectrum of the connectivity matrix and the overlaps between eigenvectors, both of which are modified when η is changed. To disentangle the effects of the changes in the eigenspectrum and the changes in eigenvector overlaps, in this section we compare the autocorrelation we
derived in Sec. <ref> with the autocorrelation we
would obtain if we assumed that the eigenvectors of J⃗ were
orthogonal. If that were the case, the autocorrelation is
computed as a sum of decoupled contributions associated with the
different eigenvalues and, in particular, the distribution of
eigenvectors would play no role in the result (see
Appendix <ref> for
details). Figure <ref>a shows the predicted autocorrelations,
both including and excluding the contribution from the overlap (<ref>). As expected, both predictions coincide for
η=1 and they increasingly depart from each other for decreasing
values of η. To better characterize this difference we show the
variance C(0) as a function of the symmetry parameter for several
values of the spectral gap [Fig. <ref>b]. Remarkably, the
variance decreases with symmetry, but the opposite occurs when we
remove the contribution from eigenvector overlaps [see 'with' and
'without' overlap curves in Fig. <ref>b]. In both cases the
variance increases as spectral gaps get smaller, which is consistent
with the fact that the restoring drive towards the fixed point gets
weaker as the spectral gap gets smaller. This effect is, however, much
subtler when overlaps are not taken into
account <cit.>.

The overlap also contributes to the overall time scale of the
autocorrelation, which we define by the quantity τ̂ = ∫_0^∞
t C(t) ṭ / ∫_0^∞ C(t) ṭ. This definition
guarantees that for an exponential autocorrelation C(τ) ∝exp(-|τ|/τ_0) the overall time scale is exactly τ_0,
and provides a rough estimate of a natural time scale for autocorrelations
with more complex dependences. The numerical evaluation of T shows
that the overall time scale is systematically smaller if the
contribution of the eigenvalues is removed
[Fig. <ref>c]. Unsurprisingly, either with or without the
overlap contribution the timescale gets longer as the spectral gap
gets smaller. Note also that the overall time scale T varies
non-monotonically with the symmetry parameter [Fig. <ref>c],
unlike the asymptotic dependence shown in Fig. <ref>b.





§ DYNAMICS IN THE CHAOTIC REGIME


In the chaotic regime, the network generates its own fluctuating
activity without the need for external noise. Recall that chaotic
activity emerges as soon as the largest of the real parts of the
eigenspectrum, given by -δ and usually called spectral
abscissa, becomes positive. We follow the strategy of the preceding
section and we keep the spectral abscissa fixed while we vary the
symmetry parameter η.



The evolution of firing activities shown in Fig. <ref>a suggests that in the chaotic regime the self-generated fluctuations get slower as η increases. This slowing is accompanied by an increasing tendency of firing rates to linger around the extreme values of their dynamical range, as reflected by an increasingly bimodal distribution of currents x and rates ϕ(x) when η increases [Fig. <ref>b]. We quantified the slowing down of the fluctuations with the population-average autocorrelation. For η=0 the autocorrelation can be derived self-consistently in the limit of infinitely large networks, using the dynamical mean-field approach (, see also Appendix <ref> for a general derivation). Unfortunately, this method does not lead to a closed-form solution for the autocorrelation as soon as η >0 (see Appendix <ref>) and we have to resort to numerical estimates, summarized in Fig. <ref>c for several values of η. For completeness we also include the autocorrelation functions for fixed gain, rather than fixed spectral abscissa [Fig. <ref>d].

The numerical estimates show that the time scale associated with the autocorrelation increases strongly as a function of η and is considerably longer than in the fixed-point regime [Fig. <ref>e]. Such a slowing is rather insensitive to whether we fix the spectral abscissa or the gain, despite the fact that the variance C(0) varies far more strongly when gain is fixed [Fig. <ref>d].


Quite strikingly, for η=1 fluctuations become slower as time goes by, and our initial assumption that the activity is stationary does not hold. The population-averaged autocorrelation C(t, t+τ)=[x(t) x(t + τ)]_J,N at different points in time shows that the characteristic timescale of the autocorrelation grows with t [Fig. <ref>d], a signature of aging dynamics <cit.>. For lower values of η, the dependence on the autocorrelation on the two timescales is less clear. Due to strong finite-size effects, it is difficult to determine from simulations alone whether aging appears also when the connectivity is not fully symmetric.



§ DISCUSSION


In this work we examined the effect of partially symmetric connectivity on the
dynamics of randomly connected networks composed of rate units. We have derived
an analytical expression for the autocorrelation function in the regime of
linear fluctuations around the fixed point, and shown that increasing the
symmetry of the connectivity leads to a systematic slowing-down of the
dynamics. Numerical simulations confirm that a similar phenomenon takes place
in the chaotic regime of the network.

The impact of the degree of symmetry of the connectivity matrix on the
dynamics of neural networks has been a long-standing question in
theoretical neuroscience. Theorists initially focused on fully
symmetric networks of binary spin-like neurons <cit.> for
which tools from equilibrium statistical mechanics could be readily
applied <cit.>. After these initial studies, the realization
that brain networks are not symmetric led physicists to investigate
the dynamics of networks whose connectivity matrix has a random
antisymmetric component. It was found that departures from full
symmetry destroys spin-glass states, while retrieval states in
associative memory models were found to be robust to the presence of
weak asymmetry <cit.>.



Theorists also studied fully asymmetric networks, using rate models
<cit.>, networks of binary neurons
<cit.> and networks of spiking neurons
<cit.>. In all these models, chaotic states were shown to be
present for sufficiently strong coupling. In networks of spiking
neurons, chaotic states are characterized by strongly irregular
activity of the constituent neurons, with self-generated fluctuations
that evolve on fast time scales. Motivated by experimental findings,
recent studies have considered synaptic connectivity matrices where
bidirectionally connected pairs are overrepresented with respect to a
random network. In contrast with our model, in which no structure
exists beyond the level of pairs of neurons, these studies have
considered structured connectivity matrices in which partial symmetry
is a consequence of a larger-scale
structure. <cit.> considered a
connectivity clustered into groups of highly connected neurons and
demonstrated that clustered connectivity could lead to slow
firing-rate dynamics generated by successive transitions between up
and down states within individual clusters.
An overrepresentation bidirectional connections can also arise in
networks with broad in- and out-degree distributions, which affect the
dynamics and the stability of asychronous states in such networks
<cit.>. Other works have considered connectivities
with non-trivial second-order connectivity statistics, and studied the
resulting network dynamics. <cit.> analyzed how
the presence of connectivity patterns involving two connections (not
only bidirectionally connected pairs) affected the tendency for a
neuronal network to synchronize, while
<cit.> focused on the oscillatory
activity generated by partially antisymmetric, delayed
interactions. Taking a completely different approach,
<cit.> showed that maximizing the number of patterns
stored in a network entails an overrepresentation of bidirectionnally
connected pairs of neurons, which suggests that partially symmetric
connectivity may be a signature of optimal information storage.

An important ingredient in our analysis is the fact that partially
symmetric interaction matrices are non-normal, i.e., they are not
diagonalizable by a set of mutually orthogonal eigenvectors. The
influence of non-normal connectivity on network dynamics has recently
received a considerable attention in the neuroscience
community <cit.>. Particularly
relevant to our study is the work by
<cit.>, who quantified the effects
of non-normality on the amplitude of the autocorrelation function in
random networks. Here we extend their results by studying the full
temporal shape of the autocorrelation function and by characterizing
how this shape is affected by the partial symmetry of connections.





The present work is also related to models of disordered systems and
spin glasses <cit.>. Most
studies in that field were inspired by physical phenomena and
considered fully symmetric interaction matrices. In that context, a
major result has been the discovery of aging, the phenomenon by
which dynamics become slower the longer the system
evolves <cit.>. This
phenomenon has been observed in a broad class of complex systems
characterized by configuration spaces with extremely rugged energy
landscapes, composed of many local minima surrounded by high
barriers. In these systems a random initial condition is very likely
to set the system far from a stationary state and initiate a very slow
relaxation towards a fixed point. The relaxation takes infinitely long
for N →∞ because, loosely speaking, the longer the
system evolves, the deeper it wanders in the valleys of the energy
landscape, and the harder it becomes for it to find configurations of
lower energy <cit.>.

Whether fully symmetric interactions are necessary to observe aging
does not seem to be entirely understood, as to the best of our
knowledge only a few works seem to have considered partially symmetric
coupling <cit.>.
Fully asymmetric networks have received more attention, but they do
not exhibit any aging phenomena. Here we interpolated between fully
asymmetric and fully symmetric networks, and have been able to obtain
mathematical results only in linear networks, in the non-chaotic
regime. Interestingly, we found that the partially symmetric case is
mathematically more complex than the symmetric or asymmetric
limits. This can be seen in the form of autocorrelation function
(<ref>), which simplifies considerably when
η=0 or η=1, but also in the Dynamical Mean Field Theory
(Appendix <ref>), where a coupling between the
autocorrelation and the response function appears for η>0. This
additional complexity results from the fact that the influence of a
single neuron's activity on all the other neurons is fed back through
couplings that are correlated with the neuron's activity, due to the
partial symmetry of the connections. More specifically, the inputs
received by neuron i are given by terms ∑_j J_ijϕ(x_j),
which are themselves influenced by the activity of neuron i. As a
result, neuron i influences its own activity by an amount
proportional to the sum ∑_j J_ijJ_jiϕ(x_i), a random
number of mean ηϕ(x_i). The effect of this feedback loop is
that the individual input terms exhibit correlated fluctuations. When
η=0, the inputs received by neurons are uncorrelated and their
sum can be approximated by a Gaussian random variable whose mean and
variance can be determined
self-consistently <cit.>. At
the other extreme, when η=1, the inputs received by neurons are
correlated, but the dynamics of the network can be described as a
relaxation of an energy function and the standard machinery of
statistical mechanics can be used. For other values of η, none of
these analytical strategies can be applied and the analysis becomes
more complex. Demonstrating analytically whether aging dynamics are
present in partially symmetric, non-linear networks seems an
outstanding open problem.

Our results on the autocorrelation function in the linear network are
closely related to recent results published by
<cit.>, who used a different set of
methods to compute the power spectrum of the network activity, i.e.,
the Fourier transform of the autocorrelation function of the same
model we investigated. Unlike <cit.>, we
obtained the autocorrelation directly in real time, although our
results are fully consistent with theirs in that we obtain the same
two regimes with the same asymptotic timescales, depending on the
symmetry and the gain (or leak, in their case. cfr. Fig. 1 of
<cit.> with Fig. <ref>a).

Our work provides a potential bridge between two seemingly unrelated
observations in neuroscience. The first is the observation of strong
correlations between the synaptic strengths in pairs of cortical
pyramidal cells, the main excitatory neuronal type in cerebral cortex,
by multiple groups using in vitro electrophysiological recordings
<cit.>. These
correlations are a consequence of two features of the connectivity:
First, there exists an overrepresentation of bidirectionally connected
pairs, compared to a Erdős-Rényi network with the same connection
probability.  For instance, <cit.> found a connection
probability of c=0.116 in pairs of neurons whose somas are less than
100 μm apart, while the probability that a pair of such neurons
are connected bidirectionally is approximately 4c^2. This degree of
overrepresentation has been found in multiple cortical areas, except
in barrel cortex where no such overrepresentation exists
<cit.>. Second, synaptic connections in bidirectionally
connected pairs are on average stronger than those in unidirectionally
connected pairs, and are significantly correlated
<cit.>. These observations lead to estimates
of η∼ 0.5, a value that, according to our model, would lead to a
significant increase in autocorrelation time scales compared to a
random asymmetric connectivity.

The second is the observation of long time scales in the
autocorrelations of neuronal activity from in vivo
electrophysiological recordings (see
e.g. <cit.>). Interestingly, the time scales of these
autocorrelations increase from sensory to higher level areas such as
the prefrontal cortex. Several mechanisms have been proposed to
account for this phenomenon: differences in the level of expression of
slow NMDA receptors <cit.>, or an increase in the strength of
recurrent connectivity <cit.> which could in particular lead to the presence of multiple fixed points that can slow down the dynamics <cit.>. Our results suggest that this
increase in time scale could also be due to an increase in the degree
of symmetry of cortical connectivity. This would be consistent with
the study of <cit.>, who showed that the
overrepresentation of bidirectionally connected pairs of neurons is
significantly stronger in prefrontal cortex than in visual cortex.

From a neuroscience point of view, the model considered here is an
extremely simplified model of cortical networks because it lacks the
fundamental constraint that neurons are either excitatory or
inhibitory, and because it does not constrain firing rates to be
positive. These simplifications were made for the sake of mathematical
tractability. A few recent studies have investigated how these two
constraints influence the dynamics of such
networks <cit.>. Extending
those works to connectivity with segregated excitation and inhibition
and partial symmetry is an important direction for future work that
might be facilitated by recent developments in random matrix theory
<cit.>.



  §.§.§ Acknowledgments


We thank Johnatan Aljadeff for his comments on a previous version of the manuscript.
The research leading to these results has received funding from the People Programme (Marie Curie
Actions) of the European Union's Seventh Framework Programme FP7/2007–2013/ under REA grant
agreement 301671. This has also been funded by the Programme Emergences of City of Paris, and the program
“Investissements d'Avenir” launched by the French Government and
implemented by the ANR, with the references ANR-10-LABX-0087 IEC and
ANR-11-IDEX-0001-02 PSL* Research University.

The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.




§ DERIVATION OF THE DOUBLE COMPLEX INTEGRAL


We summarize here the derivation of the double complex integral of Eq. (<ref>). Before doing so, we sketch the derivation of the local density of the overlap done in <cit.>, as this will let us introduce some notation and pave the way for our calculation.

For a complex variable z=x +  y, with x and y real and with
conjugate z̅=x -  y, we define the Wirtinger derivatives
∂/∂ z = (∂ / ∂ x - ∂/
∂ y) / 2 and ∂/∂z̅ = (∂ / ∂ x +
∂/ ∂ y) / 2, which obey ∂ z / ∂ z =
∂z̅ / ∂z̅ = 1 as well as ∂z̅/
∂z = ∂ z / ∂z̅ = 0. The complex
differential is defined to be z≡ẓ ẓ̣̅ = 2 x̣ ỵ,
where the factor 2 comes from the Jacobian. We also define the complex Dirac delta
so that it obeys the relation ∫δ^2(z) z≡ 1, which
implies δ^2(z) = (1/2) δ(x)δ(y) given our
convention for the complex differential. Two useful identities for the δ function in the complex plane are

    δ^2(z) = 1/2π∂/∂z̅1/z = 1/2π∂/∂ z1/z̅,

which can be checked by integrating over z and applying the
complex version of Green's theorem:

    ∫(∂ v/∂ z + ∂v̅/∂z̅) z = ∮ (v dz̅ - v̅ dz),

where v and v̅ are to be considered independent functions.

The resolvent, defined for any matrix 𝐉 as (z1 - 𝐉)^-1, is a key quantity in the analysis of random matrices because it can be ensemble-averaged using standard methods and can be related to quantities of interest. So, for example, the empirical density of eigenvalues of a given 𝐉,

    ρ_J(z) = 1/N∑_i=1^Nδ(x - λ_i)
    δ(y - λ_i) = 2/N∑_i=1^Nδ^2(z -
    λ_i),

can be expressed thanks to the identities (<ref>) as

    ρ_J(z) = 1/π∂/∂z̅1/N∑_i=1^N1/z - λ_i = 1/π∂/∂z̅1/N (z1 - 𝐉)^-1.

This quantity is hard to compute for any particular realization at
finite N, but it becomes easier to handle in the limit of large N,
where all empirical densities converge to the average density

    ρ(z) = [ ρ_J(z) ]_J = 1/π∂/∂z̅[ 1/N (z1 - 𝐉)^-1]_J.

Deriving the average density, therefore, amounts to computing the function G(z) = [ (z1 -
𝐉)^-1 / N]_J in the large N limit.

The local
density of the overlap can be derived in a similar manner using the
spectral decomposition (z1 - 𝐉)^-1 = ∑_i=1^N𝐑_i (z - λ_i)^-1𝐋^†_i, where
R⃗_i and L⃗_i are the right and left eigenvectors
of J⃗, respectivelyt. If we substitute the definition of the overlap
matrix into Eq. (<ref>) and we use the
identities (<ref>) we obtain
 
    D(z_1, z_2)     = 4/N[ ∑_i,j=1^N(L⃗^†L⃗)_ij
    (R⃗^†R⃗)_jiδ^2(z_1 -
    λ_i) δ^2(z_2 - λ_j) ]_J
        = 1/Nπ^2∂/∂z̅_1∂/∂ z_2[ ∑_i,j=1^N R⃗_i
     1/z_1 - λ_iL⃗^†_iL⃗_j1/z̅_2 - λ̅_jR⃗^†_j
     ]_J
        = 1/π^2∂/∂z̅_1∂/∂ z_2[ 1/N 1/z_11 - 𝐉1/z̅_21 - 𝐉^†]_J,

and the problem reduces to computing the quantity

    G(z_1, z_2) = [ 1/N 1/z_11 - 𝐉1/z̅_21 - 𝐉^†]_J.

The expression of G(z_1, z_2) for the ensemble of Gaussian random matrices with partial symmetry was derived by
 <cit.>. The basic idea behind their calculation is to expand resolvents in power series, average over the disorder term by term, and organize the sums so that a recursive relation can be established and ultimately solved (for a thorough description of the method, see also <cit.>).
The result is a complex function that takes the value

    G(z_1, z_2) =
    1/1 - η^2
    ×( (1 - η^2)^2 + η (z_1^2 + z̅_2^2) - (1 + η^2) z_1 z̅_2/|z_1 - z_2|^2 - 1)

when both z_1 and z_2 lie inside the ellipse centered at the origin and which has major and minor radii 1 + η and 1 - η, respectively. We will call this ellipse E_η for later convenience. When z_1 and z_2 lie outside E_η we have instead

    G(z_1, z_2) = h_1 h̅_21 - h_1 h̅_2,

where

    h_1 = z_1 - √(z_1^2 - 4η)/2η,    h̅_2=z̅_2 - √(z̅_2^2 - 4η)/2η.

Right on the ellipse E_η, |h_i| = 1. When both z_1 and z̅_2 lie outside the ellipse, the function G(z_1, z_2) is analytic on z_1 and z̅_2. This analyticity implies, from (<ref>), that the local density of the overlap vanishes outside the ellipse.

We now proceed to compute A(u, τ) [Eq. (<ref>)]. Inserting the identity (<ref>) into (<ref>) leads to

    A(u, τ) = 1/4π^2∬^g(z_1 + z̅_2) u + g z_1 τ∂/∂z̅_1∂/∂ z_2 G(z_1, z_2) z_1z_2.

Because the exponential prefactor is analytic in z_1 and z̅_2,
it commutes with the two partial derivatives. We can therefore
apply Green's theorem twice to obtain

    A(u, τ) = 1/4π^2∮_E∮_E^g (z_1 + z̅_2) u + g z_1
    τ G(z_1, z_2) ẓ_̣1̣ẓ̣̅_̣2̣,

where both contour integrals are around the ellipse E, at whose boundary G(z_1, z_2) stops being analytic. To compute A(u, τ) we follow the approach of <cit.> and use the linear transformation w=(z - ηz̅) / (1 - η^2) (or, equivalently, z = w + ηw̅) to reshape the contour of integration from the ellipse E_η into the unit circle. Applying this transformation to both z_1 and z_2, the surface integrals in (<ref>) become countour integrals on the unit circle |w|^2=w w̅ = 1. On this contour we can replace every w̅ in the integrand by w^-1 and we can use the standard tools of complex analysis to carry out the integrals. We describe in more detail our derivation in the following.

We start by performing the integral over z̅_2, expressing Eq. (<ref>) in terms of w_1, w̅_1, and w_2, and replacing all w̅_2 by w_2^-1. After some simplifications, we obtain

    G(z_1(w_1), z_2(w_2)) = -1 + 1/η(α_+ - α_-)
    ×( 1 - η w_1 α_+/w_2 - α_+ - 1 - η w_1 α_-/w_2 - α_-),

where we defined the poles

    α_±(w_1) = w̅_1 + η w_1 ±√((w̅_1 + η w_1)^2 - 4η)/2 η.

These poles depend on w_1 and and can be shown to map the unit disk onto an annulus of inner radius 1 and outer radius 1/|η| [Fig. <ref>].
This information will be relevant when we use residue calculus.




The double integral (<ref>) has to be regularized because the integrand diverges at z_1 = z_2. Our regularization consists of first integrating w_2 on the unit circle while constraining w_1 to be on a concentric circle of smaller radius |w_1| = 1 - ϵ, with ϵ>0 small. Once the integral over w_2 is done, we take the limit ϵ→ 0 and perform the second integral over w_1.

Under this regularization, we decompose the double integral (<ref>) as

    A(u, τ) = lim_ϵ→ 01/4π^2∮_|w_1|=1 - ϵ^g (w_1 + ηw̅_1)(u + τ)𝒜(w_1, u)
    × (ẉ_1 + ηẉ̣̅_1)

where we used z_1 = w_1 + ηw̅_1 and we defined

    𝒜(w_1, u) = ∮_|w_2|=1^g (w^-1_2 + η w_2)u
     G(z_1(w_1), z_2(w_2))
    × (-1/w_2^2 + η) ẉ_̣2̣.

Note that here we used w̅_2 = w_2^-1 to express the integrand and the differential ẓ̣̅_2 = ẉ̣̅_2 + ηẉ_2 in terms of w_2 only. The integrand of Eq. (<ref>) contains one singularity inside the contour of integration, at w_2=0. This singularity is associated with the essential singularity from the exponential, ^1/w_2, as well as with the pole of second order 1 / w_2^2. Because we are assuming that |w_1| < 1, the poles of G(z_1(w_1), z_2(w_2)) at w_2=α_± lie outside the contour and therefore do not contribute to the integral. We are thus left with the task of computing the residue at the origin. We do that by expanding the integrand in Laurent series around w_2 = 0, using the relations
 
    ^(z^-1 + η z)t = ∑_k=-∞^∞( z √(η))^k I_k(2√(η)t),  for  z ≠ 0,
    
     (z - z_0)^-1 = -1/z_0∑_k=0^∞(z/z_0)^k,  for  |z| < |z_0|,

 with I_k(z) being the modified Bessel function of order k. We use the last power series to expand the terms (w_2 - α_±)^-1 in G(z_1(w_1), z_2(w_2)) [Eq. (<ref>)]. This power series converges because |w_2| < |α_±| when |w_1|<1, as we assume in our regularization scheme. After expanding, applying Cauchy's residue theorem, and taking the limit ϵ→ 0, we obtain
 
    𝒜(w_1, u) = -2π/gu∑_k=0^∞ w_1^-kη^-k/2k  I_k (2g√(η)u),

The final step is to compute the integral in (<ref>) with the same strategy we used for 𝒜(w_1, u, τ). In this case we express the integrand in terms of w_1 only and we expand the exponential factor in (<ref>) with the identity

    ^(z+ η z^-1)t = ∑_k=-∞^∞( z/√(η))^k I_k(2√(η)t).

We then pick the residue from the expansion and apply Cauchy's theorem. The result is

    A(u, τ) = A_1(u, τ) + A_2(u, τ),

with

    A_1(u, τ) =
    ∑_k=-∞^∞η^k/2 I_k(2g√(η)τ) [ (1 + η^2) I_k(2g (1 + η) u)
     - η(I_k-2(2g (1 + η) u) + I_k+2(2 g(1 +
    η) u)) ]

and A_2(u, τ) given by Eq. (<ref>).

The expression (<ref>) for A_1(u, τ) can be further simplified with the identity <cit.>

    ∑_k=-∞^∞^ kα J_k(w) J_k+ν (z) = ( z - w^-α/z - w^α)^ν/2
    × J_ν(√(w^2 + z^2 - 2wzcosα)),

which we can transform into a more convenient expression for our problem, using J_ν( z) = ^ν I_ν(z) and taking α = π - ( / 2) lnη so that
^α = -√(η). The identity (<ref>) then becomes

    ∑_k=-∞^∞η^k/2 I_k(w) I_k+ν (z) = ( z + wη^-1/2/z + w η^1/2)^ν/2
    × I_ν(√(w^2 + z^2 + wz(η^1/2 + η^-1/2))),

which allows us to rewrite Eq. (<ref>) as the final expression (<ref>).

The series A_2(u, τ) does not seem to have a closed expression for general η. For η=1, however, we can
exploit the identity

    ∑_k=-∞^∞ I_k(w) I_n - k(z) = I_n(w + z)

to conclude that

    A_2(u, τ) = -I_0(2g(2u + τ)) + I_2(2g(2u + τ))   (for η=1).




§ EVALUATION OF THE TIME INTEGRAL FOR LONG Τ


The exact average autocorrelation is given in (<ref>) as a time integral that can be decomposed as

    C(τ) = ∫_0^∞^-2u - τ[A_1(u, τ) + A_2(u, τ)] ụ,

where A_1(u, τ) and A_2(u, τ) are defined in Eqs. (<ref>) and (<ref>). To gain more analytical insight we will evaluate C(τ) when τ is sufficiently large, a limit that allows us to invoke Laplace's method and approximate the integral with a closed-form expression <cit.>. Before applying the limit it is convenient to express this integral in terms of a new variable ξ≡ u / τ, which is well defined for τ > 0. With this definition

    C(τ)    = ∫_0^∞^-τ (2ξ + 1)[A_1(u(ξ), τ) + A_2(u(ξ), τ)] τξ̣,
       ≡ C_1(τ) + C_2(τ).

which we split into the two terms composing the integrand. We start with the asymptotic dependence of C_1(τ), ignoring for the moment C_2(τ). The integrand of C_1(τ) contains I_0(·) and I_2(·), whose argument is large in the long-τ limit. We can therefore use the asymptotic expansion of the modified Bessel functions of order ν,

    I_ν(x) = ^x/√(2π x)[1 - 4ν^2 - 1/8x + O(x^-2)]  for x ≫ 1.

At this order, defining

    ψ(ξ;η) ≡ψ(u(ξ), τ; η) / τ = 2 √((1 + η)^2 ξ (ξ + 1) + η),

we obtain

    C_1(τ) = √(τ/2π)∫_0^∞exp{ -τ[ (2ξ + 1) - g ψ(ξ;η) ]}/√(g ψ(ξ; η))
    [ (1 + η^2)(1 + 1/8 τ g ψ(ξ; η))
     -2 η(1 + 2(1 -η)^2/ψ(ξ; η)^2)
    (1 - 15/8 τ g ψ(ξ; η)) ]  ξ̣.

with g = (1 - δ) / (1 + η).
This integral is of the form

    A(τ) = ∫_0^∞ f(ξ) ^τ b(ξ) ξ̣.

In the limit of large τ, only the infinitessimal interval around the maximum of b(ξ) contributes to the integral because the contribution of the remaining intervals is exponentially suppressed. In our particular case the maximum of b(ξ) is at

    ξ^* = 1/2( -1 + 1 - η/(1 + η)√(2δ - δ^2)).

We distinguish two cases. For large enough values of η and δ, ξ^* is negative [Fig. <ref>], which means that in the integration range [0, ∞) of Eq. (<ref>) the maximum value of b(ξ) is at ξ=0. Conversely, for low values of η and δ, the maximum of b(ξ) occurs within (0, ∞). These two cases will lead to different time dependences and will be studied separately in the following.


For values of η and δ such that ξ^* < 0, the limit τ→∞ of Eq. (<ref>) can be approximated by <cit.>

    A(τ)    ∼lim_ϵ→ 0∫_0^ϵ f(0)^τ[b(0) + b'(0) s] ṣ
       ∼∫_0^∞ f(0) ^τ[b(0) + b'(0) s] ṣ∼ -f(0)^τ b(0)/τ b'(0),

where b'(0) denotes the derivative of b(ξ) evaluated at ξ^*=0.
Applying this approximation to (<ref>), we obtain

    C_1(τ) = τ^-3/2η^-1/4/2√(π)1 + η^2/δ(1 + η) - [1 - √(η)]^2(1 + η/1 - δ)^3/2
    ×exp{-τ(1 - √(η))^2 + 2δ√(η)/1 + η}

Conversely, if η and δ are such that ξ^* > 0, the maximum ξ^* falls within the integration region and we can approximate the large-τ limit of the integral (<ref>) by <cit.>

    A(τ) ∼lim_ϵ→ 0∫_ξ^* - ϵ^ξ^* + ϵ f(ξ^*)^τ(b(ξ^*) + (s - ξ^*)^2 b”(ξ^*) / 2) ṣ
    ∼∫_-∞^∞ f(ξ^*)^τ(b(ξ^*) + (s - ξ^*)^2 b”(ξ^*) / 2) ṣ
    ∼√(2π) f(ξ^*)^τ b(ξ^*)/√(-τ b”(ξ^*)).

with b”(ξ^*) denoting the second derivative of b(ξ) evaluated at ξ^*.
In this case Eq. (<ref>) is approximately given by

    C_1(τ) = δ^-1/2(1 - η)^2/2√(2) (1 - δ)exp{-τ1 - η/1 + η√(2δ - δ^2)}.


We now turn to the asymptotic dependence of C_2(τ).
The original form for C_2(τ) is

    C_2(τ) = ∫_0^∞{^-τ (2ξ + 1)∑_k=0^∞η^k+2 (k+1)^2
     
    ×I_k + 1(2 g τ√(η)ξ)
    I_k+1(2 g τ√(η)(ξ + 1))/ητ g^2 ξ (ξ + 1)}ξ̣.

The integrand contains a product of Bessel functions that grows exponentially with τ (see Eq. (<ref>)), but this growth is kept in check by the exponential prefactor. To see this, we introduce
 the scaled modified Bessel function,

    I_ν^*(x) ≡^-x I_ν(x),

in terms of which (<ref>) becomes
 
    C_2(τ) = ∫_0^∞exp[-τ(1 - √(η))^2 + 2δ√(η)/1 + η]
     ∑_k=0^∞η^k+1 (k+1)^2 
    ×I^*_k + 1(2 gτ√(η)ξ)
     I^*_k+1(2 g τ√(η)(ξ + 1))/τ g^2 ξ(ξ+1) ξ̣.

 The integral converges because the exponential decays to zero for large τ and trumps the power-law decay of I^*_k(·). Equation (<ref>)
 also has the same form as (<ref>) and contains exactly the same exponent as in Eq. (<ref>); because the exponential attains its maximum at ξ=0, we can use the approximation (<ref>). To evaluate at ξ=0 we use the power expansion I^*_k + 1(2 g τ√(η)ξ) ξ∼ 0= (g τ√(η)ξ)^k+1, and identify f(ξ) in (<ref>) with

    f(ξ) ≡∑_k=0^∞η^3(k + 1)/2/ξ + 1 (k+1)^2 g^k-1ξ^kτ^k
     I^*_k+1(2g τ√(η) (ξ + 1)).

All the terms of f(ξ) with k>0 vanish at ξ=0, and Eq. (<ref>) reads in this case

    C_2(τ)
     ∼τ^-3/2η^5/6/4√(π)(1 + η)^5/2/(1 - δ)^3/21/(1 - √(η))^2 + 2 √(η)δ
    ×exp{ -τ(1 - √(η))^2 + 2δ√(η)/1 + η}

Summing this contribution to that in Eq. (<ref>) leads to the result reported in Eqs. (<ref>)–(<ref>).



§ AUTOCORRELATION WITHOUT OVERLAPS


Here we compute the autocorrelation ignoring the effect of the overlaps between eigenvectors. In this case, we need to compute the individual contribution of a single eigenvalue to the autocorrelation, and then sum over the contributions of all eigenvalues. We start with the one-dimensional version of Eq. (<ref>)

    x̣/ṭ = α x + σξ(t),

where the parameter α would be the (single) eigenvalue of the system, assumed to have negative real part to prevent x(t) to grow unbounded, and where ξ(t) is a source of standard Gaussian white noise. The solution of (<ref>) is

    x(t) = σ∫_-∞^t^α (t - s)ξ(s) ṣ,

from which we can derive the autocorrelation:

    x(t) x(t + τ)    = σ^2 ∫_-∞^t∫_-∞^t + τ^α (2t + τ - s - u)ξ(s)ξ(u) ṣ ụ
        = - σ^2^ατ/2α≡σ^2 C_α(τ).

The eigenvalue α determines both the time scale and the amplitude of the autocorrelation. We set the overall factor σ^2 to 1, without loss of generality.

The average autocorrelation for the high-dimensional system in the absence of overlaps is the sum of (<ref>) over all the eigenvalues. In the large-N limit we would have

    C(τ) = 1/N∑_i=1^N C_α_i(τ) ∫ C_α(τ)ρ(α) α̣α̣̣̅,

where ρ(α) is the probability density of eigenvalues and the integral is on the complex plane. For the system (<ref>) and for the connectivity matrices we consider, the density of the eigenvalues α is uniform and has support on an ellipse centered at z=-1 with major radius g(1 + η) and minor radius g(1 - η). The integral (<ref>) can be computed in that case and reads

    C(τ) = -1/π g^2 (1 - η^2)∫_E^ατ/2α α̣α̣̣̅,

where we used Eq. (<ref>) and the prefactor is the constant value that ρ(α) takes on the elliptic support E. To evaluate the integral we use the parametrization

    α = -1 + r (1 + η) cosθ +  r (1 - η) sinθ

and integrate over r ∈ [0, g] and θ∈ [0, 2π]. Noting that

    ∫_Eα̣α̣̣̅ = (1 - η^2) ∫_0^2π∫_0^g r ṛ θ̣,

Eq. (<ref>) becomes

    C(τ) = 1/π g^2∫_0^2π∫_0^gexp{-τ[1 - r ψ (θ) ]}/2 [1 - r ψ(θ)] r ṛ θ̣.

where for convenience we defined

    ψ(θ) ≡ (1 + η)cosθ +  (1 - η) sinθ.

The integral (<ref>) is hard to compute, but we can make progress by taking the derivative of C(τ) with respect to τ,
 
    C'(τ) = -^-τ/2π g^2∫_0^2π∫_0^g^τψ(θ)r r ṛ θ̣,

 which is easier to evaluate. Equation (<ref>) can be integrated over r by parts, yielding an integral over θ only that, excluding prefactors, reads
 
    B(g, τ) ≡∫_0^2π{g ^τψ(θ) g/τψ(θ) - ^τψ(θ) g - 1/τ^2 ψ^2(θ)} θ̣.

Again, this integral is hard to compute but we can use the same trick we used before, noting that the partial derivative of A(g, τ) with respect to g simplifies considerably,

    ∂ B(g, τ)/∂ g = g ∫_0^2π^τψ(θ) g θ̣ = 2 π g I_0(2gτ√(η)),

where in the last equation we used <cit.>. We recover the expression for A(g,τ) by integrating along g, with the initial condition A(0,τ)=0:

    B(g,τ) = 2π∫_0^g x I_0(2x τ√(η)) x̣ = 2 π/4 τ^2 η∫_0^2 g τ√(η) y I_0(y) ỵ.

The last integral can be computed with the help of the recurrence relation zI_0(z) = z I'_1(z) + I_1(z). An integration by parts of the term zI_1'(z) leads to the final identity ∫ x I_0(x) x̣ = x I_1(x) and therefore to

    B(g,τ) = π g/τ√(η) I_1(2 τ g √(η)).

Equation (<ref>) then reads

    C'(τ) = -^-τ/2τ g √(η) I_1(2 τ g √(η)),

which we have to integrate to recover C(τ). Such an integration is subject to the initial condition C(0):

    C(0)    = 1/2π g^2∫_0^2π∫_0^gr/1 - r ψ(θ) ṛ θ̣.
       = -1/2π g^2∫_0^2π{g/ψ(θ) + 1/ψ^2(θ)ln[1 - g ψ(θ)]} θ̣,

which can be evaluated numerically.





§ SUMMARY OF THE DYNAMIC MEAN FIELD DERIVATION


The starting point of the calculation is the moment generating functional for the state variables x_i(t) obeying Eq. (<ref>). We consider the more general case where the activation variable is driven by recurrent inputs as well as independent external white noise:

    ẋ_i(t) = -x_i(t) + g ∑_j=1^N J_ij r_j(t) + σξ_i(t),   i=1,…,N

where we defined r_j(t) ≡ϕ(x_j(t)) to simplify the notation. The white noise sources ξ_i(t) have zero mean and unit variance.
The moment generating functional for such a system can be shown to be 

    Z[l, l̃; 𝐉] = ∫𝒟x(t) 𝒟x̃(t) exp(-S[x, x̃; 𝐉]
     + ∑_i=1^N∫l̃_i(t) x_i(t) ṭ + ∑_i=1^N∫ l_i(t) x̃_i(t) ṭ, )

where 𝒟x(t) 𝒟x̃(t) = ∏_i=1^N𝒟x_i(t) 𝒟x̃_i(t) is the functional measure for all possible paths for all variables and we introduced the action

    S[x, x̃; 𝐉] = ∑_i=1^N∫x̃_i(t)
    ×{ẋ_i(t) + x_i(t) - g ∑_j=1^N J_ij r_j(t) - σ^2/2x̃_i(t) } ṭ

In this definition we assume that the auxiliary fields x̃(t) are purely imaginary, so we do not have to write explicit imaginary units all along. By construction the generating functional satisfies the normalization condition Z[0,0;𝐉] = 1. The fact that Z[0,0; 𝐉] does not depend on 𝐉 allows us to compute the quenched average directly on Z <cit.>,

    Z[l, l̃] ≡∫Z[l, l̃; 𝐉]/Z[0, 0; 𝐉] P̣(̣𝐉̣)̣ = ∫ Z[l, l̃; 𝐉] P̣(̣𝐉̣)̣,

which simplifies considerably the average, now reduced to computing [exp(-S[x,x̃, 𝐉])]_J. To do so, we use the decomposition of partially symmetric connectivity matrices

    J_ij = J_ij^s + k J_ij^a,

where J_ij^s = J_ji^s, J_ij^a = -J_ji^a, and where both J_ij^s and J_ij^a are Gaussian random variates with zero mean and variance

    [ (J_ij^s)^2 ]_J = [ (J_ij^a)^2 ]_J = 1/N1/1 + k^2

so that [J_ij^2]_J = J^2 / N. With these matrix decompositions, the correlation between bidirectional weights is <cit.>

    [J_ij J_ji]_J = 1/N1 - k^2/1 + k^2,

which must equal η / N by our definition of η. This leads to the relation k^2 = (1 - η)/(1 + η). To integrate over the disorder we use the Gaussian measures:

    P̣(𝐉^s)     = ∏_i ≤ jP̣(J_ij^s) ∝exp{-N/1+η∑_i ≤ j (J^s_ij)^2}𝐉̣^̣ṣ, 
    P̣(𝐉^a)     = ∏_i < jP̣(J_ij^a) ∝exp{-N/1+η∑_i < j (J^a_ij)^2}𝐉̣^̣ṣ,

with 𝐉̣^̣ṣ = ∏_i ≤ jJ̣_̣ịj̣^̣ṣ and 𝐉̣^̣ṣ = ∏_i < jJ̣_̣ịj̣^̣ạ.
We will ignore the contribution of diagonal elements of the synaptic matrix because it is negligible in the limit of large N. We can now integrate out the terms linear in J_ij that appear in Eq. (<ref>), by separating symmetric and antisymmetric components. Excluding prefactors and time integrals, these terms are of the form

    L(𝐉, t) ≡∑_i, j
    i≠ jx̃_i(t) J_ij r_j(t) =
     ∑_i, j
    i≠ jx̃_i(t) [J^s_ij + k J^a_ij] r_j(t)
    
     = ∑_i < j{J^s_ij[x̃_i(t) r_j(t) + x̃_j(t) r_i(t)]
     
     + k J^a_ij[x̃_i(t) r_j(t) - x̃_j(t) r_i(t)]}

so that

    ∫exp{g ∫ L(𝐉, t) ṭ} P̣(𝐉^s) P̣(𝐉^a)
    
     = exp{g^2/2N∑_i, j 
     i ≠ j∬{[x̃_i(t) r_j(t) x̃_i(t') r_j(t')]
    
     + η[x̃_i(t) r_j(t) x̃_j(t') r_i(t')]} ṭ ṭ'̣},

where we used the property that, for a Gaussian variable z of zero mean and variance σ^2, the expected value of exp(λ z) is exp(λ z)_z = exp(λ^2σ^2/2), which can be checked by completing the square in the exponential.

Putting back together all the pieces, the average generating functional (<ref>) is therefore

    Z[l, l̃] =
     ∫𝒟x(t) 𝒟x̃(t) exp(-S_0[x(t), x̃(t)]
     + σ^2/2x̃·x̃ + l̃· x_i + l ·x̃
    
    + g^2/2N∑_i, j
     i ≠ j∬{[x̃_i(t) r_j(t) x̃_i(t') r_j(t')]
    
     + η[x̃_i(t) r_j(t) x̃_j(t') r_i(t')]} ṭ ṭ'̣)

where we defined the free action

    S_0[x, x̃] ≡∑_i=1^N∫x̃_i(t)[ẋ_i(t) + x_i(t)] ṭ.

and we introduced the notation

    f · g ≡∑_i=1^N∫ f_i(t) g_i(t) ṭ.


As a result of averaging out the disorder, we obtained a coupling involving four fields with different indices and at different times. To proceed it is convenient to introduce auxiliary fields that involve terms local in space (i.e., with the same index) but not in time:

    q_1(t, t')    = g^2/N∑_i=1^Nx̃_i(t) x̃_i(t'),    
     q_2(t, t')    = g^2/N∑_i=1^N r_i(t) r_i(t'),
    
     q_3(t, t')    = g^2/N∑_i=1^Nx̃_i(t) r_i(t'),    
     q_4(t, t')    = g^2/N∑_i=1^N r_i(t) x̃_i(t'),

so Eq. (<ref>) now reads

    Z[l, l̃] =
    ∫𝒟x(t) 𝒟x̃(t) (∏_α=1^4N/g^2 𝒟q_α)
    ×δ(N/g^2 q_1 - ∑_i=1^Nx̃_i(t) x̃_i(t')) δ(N/g^2 q_2 - ∑_i=1^N r_i(t) r_i(t')) 
    ×δ(N/g^2 q_3 - ∑_i=1^Nx̃_i(t) r_i(t')) δ(N/g^2 q_4 - ∑_i=1^N r_i(t) x̃_i(t'))
    ×exp(-S_0[x, x̃] + σ^2/2x̃·x̃ + l̃· x_i + l ·x̃
     
    + N/2 g^2∬{ q_1(t, t') q_2(t, t') + η  q_3(t, t') q_4(t, t') } ṭ ṭ'̣)

We now express the Dirac functionals in their integral representation. The first Dirac functional appearing in Eq. (<ref>) can be written as

    δ(N/g^2 q_1 - ∑_i=1^Nx̃_i(t) x̃_i(t')) = 1/2π∫𝒟q̂_1(t, t') 
    ×exp{∬q̂_1(t, t') [N/g^2 q_1(t, t') - ∑_i=1^Nx̃_i(t) x̃_i(t') ] ṭ ṭ' },

where the integral over q̂ is understood to be along the imaginary axis. The other Dirac functionals in Eq. (<ref>) are rewritten analogously. Equation (<ref>) then becomes

    Z[l, l̃] =
    ∫𝒟X 𝒟Q exp(-S_0[x, x̃] + σ^2/2 x̃·x̃ + l̃· x + l ·x̃
    
    + N/g^2∬{∑_α = 1^4q̂_α(t, t') q_α(t, t')
    
    + 1/2[ q_1(t, t') q_2(t, t') + η  q_3(t, t') q_4(t, t')] 
    
    - g^2/N∑_i=1^N[ q̂_1(t, t') x̃_i(t) x̃_i(t') +
    q̂_2(t, t') r_i(t) r_i(t')
     + q̂_3(t, t') x̃_i(t) r_i(t') + q̂_4(t, t') r_i(t) x̃_i(t')]} ṭ ṭ'̣)

where we introduced the shorthand notation

    𝒟 Q    ≡∏_α=1^41/2πN/g^2𝒟 q_α𝒟q̂_α,
    𝒟 X    ≡𝒟x(t) 𝒟x̃(t) = ∏_i=1^N𝒟x_i(t) 𝒟x̃_i(t).


Equation (<ref>) can now be expressed as <cit.>

    Z[l, l̃] = ∫𝒟Q ^N f(q, q̂, x, x̃),

where

    f(q, q̂, x, x̃)    ≡ G(q, q̂) + 1/Nlog∫𝒟X exp[ℒ(q, q̂, x, x̃)],
    
     G(q, q̂)    ≡1/g^2∬{∑_α=1^4 q_αq̂_α + 1/2[ q_1 q_2 + η  q_3 q_4 ] } ṭ ṭ'̣,
    
     ℒ(q, q̂, x, x̃)    ≡ -S_0[x, x̃] + σ^2/2 x̃·x̃ + l̃· x + l ·x̃
        - ∑_i=1^N∬[ q̂_1(t, t') x̃_i(t) x̃_i(t') +
    q̂_2(t, t') r_i(t) r_i(t') 
        + q̂_3(t, t') x̃_i(t) r_i(t') + q̂_4(t, t') r_i(t) x̃_i(t')] ṭ ṭ'̣.

In the limit of large N we can apply the saddle-point method to Eq. (<ref>), which amounts to making the following approximation
 
    Z[l, l̃] = ∫𝒟Q ^N f(q, q̂, x, x̃)≈^N f(q^0, q̂^0, x, x̃)

 where q^0 and q̂^0 are the values that extremize f. Requiring δ f/δq̂_α = 0 leads to
 
    q^0_1(t, t')    = g^2/N∑_i=1^N⟨x̃_i(t) x̃_i(t')⟩_ℒ , 
     
     q^0_2(t, t')    = g^2/N∑_i=1^N⟨ r_i(t) r_i(t')⟩_ℒ ,
    
     q^0_3(t, t')    = g^2/N∑_i=1^N⟨x̃_i(t) r_i(t')⟩_ℒ ,
    
     q^0_4(t, t')    = g^2/N∑_i=1^N⟨ r_i(t) x̃_i(t')⟩_ℒ ,

 with the average ⟨·⟩_ℒ defined as

    𝒪_ℒ≡∫𝒪(X) exp[ℒ(X)] X̣∫exp[ℒ(X)] X̣.

Similarly, from the saddle-point conditions for q_α we obtain

    q̂^0_1(t, t')    = -1/2 q^0_2(t, t'),    q̂^0_2(t, t')    = -1/2q^0_1(t, t'),
    q̂^0_3(t, t')    = - η/2 q^0_4(t, t'),    q̂^0_4(t, t')    = - η/2 q^0_3(t, t').


Now the right hand side of Eq. (<ref>) reads

    Z[l, l̃] =
    ∫𝒟X exp(-S_0[x, x̃] - S_int[x, x̃] + l̃· x + l ·x̃)

with

    S_int[x, x̃] ≡1/2∬{N/g^2[q^0_1(t, t') q^0_2(t, t') 
    η  q^0_3(t, t') q^0_4(t, t') ]
    
     - ∑_i=1^N{[q^0_2(t, t') + σ^2 δ(t - t')]x̃_i(t) x̃_i(t')
    
    + q^0_1(t, t') r_i(t) r_i(t') η  q^0_4(t, t') x̃_i(t) r_i(t')
    η  q^0_3(t, t') r_i(t) x̃_i(t')}} ṭ ṭ'̣.


The auxiliary fields defined in Eqs. (<ref>)–(<ref>) are related to physically observable quantities.
 First, q_2^0(t, t') is related to the population-averaged autocorrelation function

    C(t, t') ≡1/N∑_i=1^N⟨ r_i(t) r_i(t') ⟩ ,

by q_2^0 (t, t') = g^2 C(t, t').

Second, the auxiliary fields q_3^0(t, t') and q_4^0(t, t') are related to the so-called response function, which characterizes the response of the system when it is perturbed by a weak field. More specifically, in our context the response function at site i would be

    G(t, t') ≡. δ⟨ r_i(t) ⟩/δ h_i(t')|_h=0,

where h_i(t') is a time-dependent external field, and angular brackets denote the average over the effective action S[x,x̃] = S_0[x,x̃] + S_int[x,x̃] that appears in Eq. (<ref>). Note that from the definition of response function G(t, t') has to be 0 whenever t < t', due to causality. To see the link between G(t, t') and q_3(t,t') and q_4(t,t'), we add an external field h_i(t) for each neuron in Eq. (<ref>), and evaluate (<ref>). With the new field the action becomes S_h[x,x̃] = S[x,x̃] -∑x̃_i(t) h_i(t) and

    . δ⟨ r_i(t) ⟩/δ h_i(t')|_h=0    = δ/δ h(t')∫𝒟X  r_i(t) exp(-S_h[x, x̃]) |_h=0
        = -⟨ r(t) . δ S_h/δ h_i(t')|_h=0⟩ = ⟨ r_i(t) x̃_i(t') ⟩.

Defining the population-averaged response function as

    G(t, t') ≡1/N∑_i=1^N⟨ r_i(t) x̃_i(t') ⟩ ,

we obtain q_4^0(t, t') = q_3^0(t', t) = g^2 G(t, t').

As for q_1^0(t,t'), it can be shown that the presence of vertices like
r_i(t) r_i(t') in the action necessarily leads to violation of causality
<cit.>. We thus need to impose q_1^0(t,t')=0
to obtain a physical solution.

We can finally write the interacting action in Eq. (<ref>) in terms of the physical quantities

    S_int[x,x̃] = -∑_i=1^N∬{1/2Γ(t, t') x̃_i(t) x̃_i(t') 
    + η g^2  G(t, t') x̃_i(t) r_i(t') } ṭ ṭ'̣,

where we ignored the term containing G(t, t')G(t, t'), which vanishes due to causality, and where we defined

    Γ(t, t') ≡ g^2 C(t, t') + σ^2 δ(t - t') .

Note that the final action involves only interactions that are local in space, which implies that all units are equivalent. This equivalence comes as no surprise, because all units are equivalent once we average over all realizations of the connectivity matrix. We can thus drop the irrelevant indices i and focus on the single relevant dynamical variable x(t).



 §.§ Equation of motion for the average activity

The local action S[x, x̃] = S_0[x, x̃] +S_int[x, x̃], with S_0 and S_int given by
Eqs. (<ref>) and (<ref>), has the form

    S[x, x̃] = ∬{x̃(t) G^-1_F(t - t') x(t') - η g^2 x̃(t) G(t, t') r(t') 
    
     - 1/2Γ(t, t') x̃(t)x̃(t')} ṭ ṭ'̣

where for later convenience we have introduced the inverse of the free propagator, G^-1_F(t - t'). The free propagator G_F(t - t') is just the Green's function associated with the operator /ṭ + 1, 
    (/ṭ + 1) G_F(t - t') = δ(t - t'),

and is related to its inverse through ∫ G_F^-1(t - s)G_F(s - t') ṣ = δ( t - t'), which is automatically satisfied if

    G_F^-1 (t - t') = δ(t - t') (/ṭ'̣ + 1).


From the original stochastic system (<ref>) and its associated
Martin-Siggia-Rose-Janssen-de Dominicis (MSRJD)
action (<ref>), we infer that the equation of motion
associated with the action (<ref>) is

    ẋ(t) = - x(t) + η g^2 ∫_-∞^t G(t, s) r(s) ṣ + φ(t),

where φ(t) is a source of noise with autocorrelation

    ⟨φ(t) φ(t') ⟩ = Γ(t, t') = g^2 C(t, t') + σ^2 δ(t - t').

This relation has to be consistent with the dynamics generated by Eq. (<ref>), that is, the noise φ(t) has to be such that the firing activity r(t) has autocorrelation C(t, t').

We can go further and write a self-consistent relation involving the two-point functions C(t,t') and G(t, t'). A starting point to derive them are the identities

    δ x(t)/δ x(t') = δx̃(t)/δx̃(t') = δ(t - t'),  δx̃(t)/δ x(t') = δ x(t)/δx̃(t') = 0

from which we can obtain relations such as

    ⟨δ x(t)/δ x(t')⟩   ≡∫𝒟X δ x(t)/δ x(t')exp{-S[x(t), x̃(t)]}
        = ⟨ x(t) δ S/δ x(t')⟩ = δ(t - t').

Other relations follow analogously:

    ⟨ x(t) δ S/δ x(t')⟩    = δ(t - t'),    ⟨ x(t) δ S/δx̃(t')⟩    = 0,
    ⟨x̃(t) δ S/δx̃(t')⟩    = δ(t - t'),
        ⟨x̃(t) δ S/δ x(t')⟩    = 0,


We now apply the identities (<ref>) and (<ref>) for the action (<ref>). In particular, we use the identities involving

    δ S/δx̃(t') = ẋ(t) + x(t) - η g^2 ∫_-∞^t G(t,s) r(s)  ṣ - ∫Γ(t, s) x̃(s) ṣ

and we define the autocorrelation and response function of the activation field x(t)

    Δ(t, t') ≡x(t) x(t'),   R(t, t') ≡x(t) x̃(t').

The last equation in (<ref>) and the first equation in (<ref>) then become, respectively,
 
    ∂/∂ tΔ(t, t')     = -Δ(t, t') + σ^2 R(t', t)
     η g^2 ∫^t_0 G(t, s) r(s)x(t') ṣ
            + g^2 ∫^t'_0 R(t',s) C(t,s) ṣ,
    ∂/∂ t R(t, t')     = -R(t, t') + δ(t - t')
            + η g^2 ∫^t_t' G(t, s) G(s,t') ṣ,

where in (<ref>) we have used x̃(t)x̃(t') = 0. It can be shown that the remaining identities in Eqs. (<ref>) and (<ref>), which involve δ S / δ x, do not provide additional information <cit.>. Note that Δ(t,t') has a cusp at t=t' due to the term σ^2 R(t', t), which from (<ref>) we know it must be of the form R(t,t') ∝Θ(t - t'), with Θ(t) being the step function. More specifically,

    [∂/∂ tΔ(t, t') ]_t'= t^+^t'= t^- = σ^2 [R(t',t)]_t'= t^+^t'= t^- = -σ^2.

Moreover, the symmetry of Δ(t, t') around t=t' implies lim_t'→ t^-∂_t Δ(t, t') = - lim_t'→ t^+∂_t Δ(t, t'), which leads to the relation lim_t'→ t^-∂_t Δ(t, t') = -σ^2 / 2. The amplitude of external noise thus determines the slope of the autocorrelation of x(t) at zero time lag. This is the only dependence on σ^2 of the solutions of (<ref>) and (<ref>).

Equations (<ref>) and (<ref>) cannot be solved in a closed-form except for η=0 <cit.>, but perturbative solutions can be found by expanding the nonlinearity r(t)=ϕ(x) in power series of x(t) and then solving the resulting hierarchy of equations, which involve correlations and response functions of increasingly larger order. The problem becomes unwieldly except for the linear case where r(t)=x(t). In that case, C(t,t')=Δ(t,t'), G(t,t')=R(t, t'), and Eqs. (<ref>) and (<ref>) form a closed system of integro-differential equations:

    ∂/∂ tΔ(t, t')     = -Δ(t, t') + σ^2 R(t', t) +
    η J^2 ∫^t_0 R(t, s) Δ(s, t') ṣ
        + g^2 ∫^t'_0 R(t',s) Δ(t,s) ṣ,
    ∂/∂ t R(t, t')     = -R(t, t') + δ(t - t')
        + η g^2 ∫^t_t' R(t, s) R(s, t') ṣ.



apsrev4-1




