Is the cortical dynamics ergodic? A numerical study in partially symmetric networks of spiking neurons

Is the cortical dynamics ergodic? A numerical study in partially symmetric networks of spiking neurons | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Is the cortical dynamics ergodic? A numerical study in partially symmetric networks of spiking neurons View ORCID Profile Ferdinand Tixidre , Gianluigi Mongillo , Alessandro Torcini doi: https://doi.org/10.1101/2025.07.28.666831 Ferdinand Tixidre 1 Laboratoire de Physique Théorique et Modélisation, UMR 8089, CY Cergy Paris Université, CNRS , Cergy-Pontoise, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Ferdinand Tixidre For correspondence: ferdinand.tixidre{at}cyu.fr Gianluigi Mongillo 2 Sorbonne Université, INSERM, CNRS, Institut de la Vision , F-75012 Paris, France 3 Centre National de la Recherche Scientifique , Paris, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Alessandro Torcini 1 Laboratoire de Physique Théorique et Modélisation, UMR 8089, CY Cergy Paris Université, CNRS , Cergy-Pontoise, France 4 CNR - Consiglio Nazionale delle Ricerche - Istituto dei Sistemi Complessi , via Madonna del Piano 10, 50019 Sesto Fiorentino, Italy 5 INFN Sezione di Firenze , Via Sansone 1, 50019 Sesto Fiorentino, Italy Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abstract Full Text Info/History Metrics Preview PDF Abstract Cortical activity in-vivo displays relaxational time scales much longer than the membrane time constant of the neurons or the deactivation time of ionotropic synaptic conductances. The mechanisms responsible for such slow dynamics are not understood. Here, we show that slow dynamics naturally and robustly emerges in dynamically-balanced networks of spiking neurons. This requires only partial symmetry in the synaptic connectivity, a feature of local cortical networks observed in experiments. The symmetry generates an effective, excitatory self-coupling of the neurons that leads to long-lived fluctuations in the network activity, without destroying the dynamical balance. When the excitatory self-coupling is suitably strong, the same mechanism leads to multiple equilibrium states of the network dynamics. Our results reveal a novel dynamical regime of the collective activity in spiking networks, where the memory of the initial state persists for very long times and ergodicity is broken. I. INTRODUCTION In the active cortex, neurons exhibit significant correlations in their spike counts, even when these counts are separated by lags as long as 1 second or more [ 1 – 3 ]. The dynamics of the membrane depolarization and the activation/deactivation of synaptic conductances have, instead, typical time scales on the order of 10ms [ 4 ]. Hence, these processes cannot account for the long-lived fluctuations of spiking activity observed in experiments. Theoretically, there are two “minimal” scenarios to account for them. In one scenario, the long time scales are simply inherited from some slow underlying biophysical processes. For instance, both neurons [ 5 , 6 ] and synapses [ 7 ] exhibit dynamics on multiple time scales, some of the order of 1-10s. We are not aware of any systematic investigation of this scenario. We note, however, that the relationship between the time scales of the microscopic (neuron/synapse) and the macroscopic (network) dynamics is far from trivial; slow microscopic dynamics does not automatically result in slow macroscopic dynamics [ 8 ]. In the other scenario, the long time scales emerge from the interaction of fast biophysical processes. In the following, we focus on this scenario. The theory of balanced networks [ 9 – 12 ] elegantly accounts for many non-trivial features of the cortical activity. According to this theory, in an active cortical network, the excitatory and the inhibitory drive to single neurons are both large, but the level of activity in the network dynamically adjusts so that the net drive is around threshold. This leads to a dynamical state—the balanced regime—characterized by temporally-irregular spiking at the single-neuron level with low, heterogeneous firing rates [ 13 – 16 ] and by weak pairwise correlations in the spiking of different neurons [ 12 , 17 ], reproducing many characteristics of the patterns of activity observed in the cortex. Long relaxation times, however, are not a generic feature of random networks operating in the balanced regime; additional mechanisms have to be postulated. There are presently three main proposals to explain the emergence of long time scales in cortical networks. According to one proposal, cortical networks operate close to the edge of stability [ 18 – 20 ]. Depending on the level of variability in the synaptic efficacies (synaptic gain), random networks of rate neurons display either a unique, stable fixed point or a chaotic attractor. Adjusting the synaptic gain suitably close to its critical level, one can achieve arbitrarily slow relaxational dynamics. As the network keeps operating in the balanced regime on both sides of the transition, pairwise correlations remain weak. However, such a transition does not exist in spiking networks [ 18 , 19 ]. More precisely, a “sharp” transition is recovered only when the synaptic time constant becomes infinitely long, so that spiking neurons effectively operate as rate elements [ 18 , 19 , 21 ]. Considering typical firing rates (∼1Hz) and synaptic time constants (10 − 100ms) in the cortex, one expects a significant “smoothing” of the transition; in particular, the autocorrelation time remains finite regardless of the level of the synaptic gain [ 19 ]. In fact, numerical simulations show that the autocorrelation time remains of the order of the synaptic time constant when the network parameters, in particular synaptic time constants, are chosen within ranges of biological credibility [ 18 ]. According to another proposal, cortical networks are metastable [ 22 – 24 ]; see Ref. [ 25 ] for a general overview. With suitable synaptic structuring, model networks feature multiple steady states of activity, differing in the spatial distributions of “active” (high-rate) and “inactive” (low-rate) neurons. Noise , then, induces transitions between these states. Thus, single neurons exhibit fluctuations in their firing rate on time scales comparable with the lifetime of network’s steady states. For suitably low levels of noise, these lifetimes can be significantly longer than the microscopic time constants. Furthermore, if pairs of neurons are co-active only in a small fraction of the steady states (sampled by the dynamics during the typical duration of the experiments), measured pairwise correlations remain weak [ 23 , 24 ]. Achieving levels of noise suitable to match the experimental observations, however, requires significant finetuning; see Ref. [ 26 ] for a discussion of this point. The reason is easy to understand. Life-times much longer than the microscopic time scales require the existence ofquasi-equilibrium states a la Kramers [ 27 ]. Accordingly, the transition rate between the different states will, to a good approximation, depend exponentially on the level of noise; see, e.g., Ref [ 23 ]. The level of noise, in turn, is controlled by the number of “active” neurons in a steady state, the strength of the recurrent synaptic couplings and the level of the external drive to the network. It is presently unclear, if neural and/or synaptic plasticity mechanisms can control these parameters with the required, exponential, precision; see, e.g., Ref [ 28 ] for a proposal based on synaptic plasticity. A more recent and, hence, much less investigated proposal is that long time scales result from the symmetry in the synaptic connectivity [ 29 – 31 ]. In fact, local connectivity (i.e., on spatial scales of the order of 100 µ m) in the cortex is not random. Rather, reciprocal synaptic connections are both over-represented, as compared to random networks, and stronger than average [ 32 – 36 ]. On first principles, then, one would expect a slowing down of the dynamics because of the effective positive selfcoupling (i.e., the so-called reaction term ; see Sec. II) generated by the symmetry. Indeed, the slowing down is observed in model networks, both of rate [ 29 , 31 , 37 ] and of spiking neurons [ 30 ]. A detailed investigation of the effects of symmetry on the dynamics of spiking networks, however, is still lacking. We provide here such an investigation. We consider networks of inhibitory leaky-integrate-and-fire neurons, with arbitrary levels of symmetry in the synaptic connectivity and operating in the balanced regime, as a minimal model of a local cortical network. At increasing levels of symmetry, and for suitably large synaptic efficacies, we find that the neurons display large and long-lived fluctuations in their spiking rate, similarly to what is observed in clustered model networks and in experiments [ 22 – 24 ]. As the symmetry level is increased, the lifetime of the fluctuations becomes comparable to the duration of our simulations (∼ 10 4 seconds). Notably, this behavior is independent of the synaptic time constant and occurs in partially symmetric networks. In the parameter region where these exceptionally long-lived fluctuations exist, interpolating between spiking and rate neurons by progressively increasing the synaptic time constant, we find that the network dynamics exhibit multiple, stable fixed points. Taken together, our results suggest that, for suitable levels of synaptic symmetry and of synaptic efficacies, the dynamics of a balanced network of spiking neurons is non-ergodic. In this parameter region, different initial conditions lead to different asymptotic states of activity. II. RESULTS As a paradigmatic example of a network operating in the balanced regime, we consider N randomly and sparsely connected inhibitory neurons receiving a constant, excitatory input from the outside [ 18 , 19 , 38 , 39 ]. Single-neuron dynamics is described by the standard leaky integrate-and-fire (LIF) mechanism [ 4 ]. Each neuron randomly makes K (≪ N ) synaptic contacts, on average, with the other neurons. In order to obtain a dynamically-balanced dynamics [ 10 , 12 ], we assume that the efficacy of the synaptic contacts is with g > 0, so that the recurrent synaptic input is proportional to and negative, and that the external excitatory input to each neuron is with h ext > 0. We control the level of symmetry in the network by manipulating the probability of reciprocal connection between pairs of neurons [ 30 ]. The excess reciprocity, η , defined as the ratio between the probability of finding recip-rocally connected neurons and the probability of reciprocal connections in a random network (i.e., directed Erdős-Rényi with the same probability of connection K/N ) is given by (see Sec. IV for details) where 0 ≤ q ≤ 1. For q = 0, the network is simply a random network ( η = 1), while for q = 1 the network is fully symmetric ( η = N/K ; all connections are reciprocal). Experimental estimates of the probability of local connections (0.1-0.3) and of the excess reciprocity (3-5) in the cortex [ 36 ] lead to 0.4 ≤ q ≤0.9. In particular, we note that the inhibitory-to-inhibitory network in the cortex is close to fully symmetric, i.e., 0.8 ≤ q ≤0.9. We start by considering the case of instantaneous synaptic currents (see Sec.IV A for details), i.e., when a neuron fires a spike, it instantaneously reduces the volt-age of its postsynaptic targets by . In the absence of pairwise correlation in the synaptic connectivity (i.e., q = 0), our model network operates in the balanced regime when K ≫ 1. In this regime, the large excitatory drive (the external input) is dynamically offset by the comparably large inhibitory drive (the recurrent input), resulting in the unsaturated, temporally irregular, and asynchronous spiking of the neurons in the network [ 10 , 12 ]. This regime is illustrated in Fig. 1 . In Fig. 1(a) we show the distribution of firing rates in the network; as can be seen it is right-skewed and long-tailed, and, it can be well approximated by a log-normal distribution (not shown) [ 14 , 16 ]. Moreover, the population-averaged firing rate, ⟨ ν ⟩ = 2.1Hz, is in good agreement with the value, ν B , predicted by the balance condition in the limit K → ∞, that is, ν B = h ext / ( τ m g ) = 2 Hz where τ m = 10 ms is the membrane time constant [ 38 , 39 ]. As expected in the balanced regime, spiking is irregular, as assessed by the coefficient of variation (CV) of the inter-spike intervals of single neurons ( Fig. 1(d) ), and asynchronous, as assessed by the zero-lag cross-correlation of the spike counts of different neurons ( Fig. 1(g) ) [ 12 , 17 ]. Download figure Open in new tab FIG. 1: (a)–(c) PDFs of the single-neuron firing rates { ν i }. Dashed lines denote the mean and the dotted lines represent the first and third quartiles. The insets display the activity of thirty randomly selected neurons. (d)–(f) PDFs of the coefficients of variation { CV i } of individual neuron firing rates. The dotted vertical lines indicate the means of the PDFs. (g)–(i) PDFs of the pairwise correlation coefficients ρ ij of the neuronal firing activity. The dotted horizontal lines indicate the full widths at half maximum of the PDFs. The solid black curves in (d)–(f) and (g)–(i) represents Gaussian fits to the data, the corresponding mean and standard deviation are denoted by µ and σ , respectively. The top row (blue) corresponds to q = 0.0; the middle row (magenta) to q = 0.25 and the bottom one (green) to q = 0.5. The PDFs were obtained by aggregating data from 10 independent replicas. Simulation parameters: N = 4000, K = 1200, g = 5, and h ext = 0.10. Increasing the levels of symmetry up to q = 0.5 does not affect these features qualitatively; spiking activity remains unsaturated, temporally irregular and asynchronous. There are some quantitative changes, however. First, the distribution of firing rates becomes more skewed while the population-averaged firing rate remains approximately constant, as can be seen by comparing Fig. 1(b-c) with Fig. 1(a) : ⟨ ν ⟩ varies from 2.1 Hz for q = 0 to 2.3 Hz for q = 0.5. Second, the temporal irregularity of the spike trains increases, as can be seen in Fig. 1(e-f) ; the population-averaged coefficient of variation, ⟨ CV ⟩, increases from 1 for q = 0 to 1.4 for q = 0.5. However, there are no noticeable changes in the level of “asynchrony”, as can be seen in Fig. 1(h-i) . Thus, the network continues to operate in the balanced regime. As we shall explain shortly, these quantitative changes are a consequence of the effective positive feedback generated by the symmetry in the synaptic connectivity (see also Ref. [ 30 ]). First, let us discuss the mechanism generating the positive feedback. Consider neuron i ; every time it fires a spike, it delays spiking in its postsynaptic targets. Indeed, upon spiking, neuron i reduces the voltage of its postsynaptic targets by and, hence, transiently reduces their instantaneous probability of spiking by a factor (for K ≫ 1). This is simply a consequence of the fact that the neurons are inhibitory. In the presence of symmetry, many of the postsynaptic targets of neuron i are also its presynaptic neurons, that is, a fraction q (on average). A reduction in the probability of spiking of one of these neurons, then, leads to a reduction ∼ g/K in the (average) inhibitory input to neuron i . As there are ∼ qK of them, the total reduction is ∼ 1 (i.e., independent of K , when K ≫ 1). As a result, the net synaptic input received by neuron i following its spike is transiently increased, because of the (transiently) reduced probability of spiking of its inhibitory presynaptic neurons. The typical size of this positive feedback is ∼ 1 and, hence, the balanced regime is unaffected by it [ 30 ]. The presence of such a positive feedback is readily revealed by the analysis of the spike-triggered average input, STA( τ ) (see Sec. IV B for details), where τ is time relative to spike emission ( τ = 0). Let us consider first the case q = 0, where no positive feedback is expected. The corresponding STA( τ ) is the blue curve in Fig. 2(a) . As can be seen, STA( τ ) increases above baseline briefly before the spike (i.e., τ 0), indicating that there is no consistent effect of the spike on subsequent inputs. As can be seen in the inset of Fig. 2(a) , this is the same STA( τ ), for all practical purposes, that one obtains from a LIF neuron driven by white noise [ 40 ] or by shot noise (i.e., a Poissonian spike train inducing finite-amplitude postsynaptic potentials), where, by construction, there is no correlation between spikes and synaptic inputs. Download figure Open in new tab FIG. 2: (a) Spike-triggered average input STA( τ ) for q = 0.0 (blue), q = 0.25 (magenta) and q = 0.5 (green). In the inset of (a) we compare the STA obtained from the network simulations (grey), with a single neuron’s dynamics under the diffusive approximation [ 40 ] (sky blue) and subject to Poissonian shot noise (orange). (b) Firing rates of the neurons versus the time averaged (mean) input currents for different levels of symmetry in the network q = 0.0 (blue), q = 0.25 (magenta), q = 0.5 (green); the Population Frequency-Input curves are shown in the inset and obtained by fitting the logistic function to the raw data (details in [ 41 ]). (c) Cumulative distribution functions (CDFs) of the mean input currents. The curves reported in (a) were obtained by averaging over 200, 000 spike events, while data in (b) and (c) has been obtained by averaging over a time interval of 250 seconds. Same parameters used as in Fig.1 As expected, the situation changes when the level of symmetry is increased. For q = 0.25, we see that STA( τ ) remains significantly above the baseline for about 20 ms after the spike (magenta curve in Fig. 2(a) ). This “afterspike transient” in the synaptic input is further amplified, both in amplitude and duration, when q = 0.5 (green curve in Fig. 2(a) ). It is convenient to quantify the strength of the positive feedback, R , by integrating the spike-triggered average input for positive τ [ 42 ], i.e., For q = 0, we find R ≃ 0.1 ms, while R ≃ 6 ms for q = 0.25 and R ≃ 24 ms for q = 0.5. With Fig. 2(a) in mind, it becomes easier to understand the effects of the symmetry on the patterns of spiking activity, illustrated in Fig. 1 . Let us start with the temporal structure of the spike trains. The presence of the positive feedback makes spiking a self-exciting process: firing a spike increases, transiently, the probability of firing another spike shortly. The strength of the self-excitation (i.e., R defined above) increases with increasing levels of symmetry (see Fig. 2(a) ) which, in turn, leads to an increasing tendency to fire spikes in “bursts”. As a result, the mean coefficient of variation ⟨ CV ⟩ increases with q . Next, let us consider the effect on the distribution of firing rates. At parity of synaptic input, the firing rate of a neuron increases with the level of symmetry q (and hence R ). This is illustrated in Fig. 2(b) , where we plot the firing rates of the neurons in the network as a function of their time-averaged synaptic input, for the same values of q (and using the same color coding) as in Fig. 1 . In the inset, we plot the population Frequency-Input (F-I) curves, obtained by fitting the data in the main panel with a logistic function [ 41 ]. However, the population-averaged firing rate ⟨ ν ⟩ cannot change dramatically. This is due to the fact that ⟨ ν ⟩ is essentially fixed (i.e., within to ν B by the balance condition, and the value of ν B is independent of q . To “solve” this apparent contradiction, the distribution of the firing rates adjusts so as to reduce the time-averaged synaptic inputs while (approximately) preserving the population-averaged firing rate. This reduction is illustrated in Fig. 2(c) . The reduction of the time-averaged synaptic inputs leads to an increase in the fraction of both low-firing and high-firing neurons with increasing q , as clearly visible in the inset of Fig. 2(b) . As a consequence of the increased tendency of neurons to fire spikes in bursts, we expect large, positive correlations between the synaptic inputs (to the same neuron) at time lags longer than the membrane time constant ( τ m = 10 ms). We also expect these correlations to persist over longer lags as q increases. This is confirmed in Fig. 3 , where we plot the population-averaged autocorrelation function of the synaptic inputs, C ( t ), for different values of q (see Sec. IV B for details). Download figure Open in new tab FIG. 3: The autocorrelations are shown as a function of the time lag for different levels q of symmetry in the synaptic connections: q = 0.0 (blue), q = 0.25 (magenta), q = 0.5 (green), q = 0.8 (orange) and q = 0.9 (red). The inset shows a zoomed-in view for q = 0.9, where different shades of the red color denotes different initial conditions. The autocorrelations were averaged over 500 selected neurons and over 3 replicas. The same parameters as in Fig.1 were used. As can be seen, consistent with the expectations, C ( t ) is non-negative, and the time it takes to decay to 0 gets longer with increasing q . For levels of symmetry up to q = 0.8, C ( t ) is well fitted by a single exponential decay, i.e., for t ≥ 0. The time constant τ C , extracted from the fit, provides a simple quantification of the decay time. For q = 0 (blue curve), we find τ C ≃ 20 ms, while for q = 0.8 (orange curve), we find τ C ≃ 100 ms. At q = 0.9 (red curve), however, the behavior of C ( t ) changes dramatically. As can be seen in the inset of Fig. 3 , C ( t ) remains significantly different from 0 at t = 50 s; three orders of magnitudes larger than τ m . Also, the decay of C ( t ) is clearly different from a single exponential decay; indeed, it features at least two time scales. Finally, C ( t ) displays noticeable fluctuations even at short lags when re-estimated (i.e., in the same network) starting from different initial conditions (lighter-red curves). For q ≤0.8, the same re-estimation procedure produces virtually identical C ( t ) (not visible in Fig. 3 since they are essentially overlapping). Thus, the analysis of C ( t ) reveals a qualitative change in the network dynamics, occurring between q = 0.8 and q = 0.9. For q = 0.8, the fluctuations in the activity decay exponentially at long time lags, suggesting that the network dynamics has reached a stationary state. For q = 0.9, the fluctuations in the activity do not decay exponentially; rather, they persist for a time comparable with the duration of the simulations, suggesting that the network dynamics has not reached a stationary state. To better understand the dynamical mechanism responsible for such long memory in the fluctuations, we simulated the dynamics in the same network starting from different initial conditions, that is, different values of neurons’ membrane potentials at time t = 0. In Fig. 4(a) , we plot the time course of the firing rate of a randomly chosen neuron for 5 different initial conditions, at q = 0.8. The firing rate is estimated in consecutive, non-overlapping bins of duration 10s. The 5 corresponding time series appear “stationary”, consistent with the analysis of C ( t ) (see Fig. 3 ). Indeed, this is expected if the spike trains for the different initial conditions were simply different realizations of the same “underlying”, constant (i.e., on time scales longer than τ C ) firing rate. To qualify this statement, for each initial condition, we have estimated the firing rate by counting all the spikes occurring between t = 0 and t = T . As can be seen in Fig. 4(b) , these estimates converge to the same value at increasing T , regardless of the initial condition. Download figure Open in new tab FIG. 4: (a),(c) The “instantaneous” firing rate for the same neuron and replicas, using time bins of Δ = 10 seconds, for 5 different replicas for q = 0.8 (orange) and q = 0.9 (red). Different shades of the main color represent different replicas. (b),(d) Running averages of the firing rate of a single selected neuron versus time. All other parameters are the same as in Fig. 1 . The same analysis for q = 0.9 is shown in Fig. 4(c) . In this case, the time series for the different initial conditions clearly do not resemble different realizations of the same “underlying”, constant firing rate. One can see large (up to ∼ 100 Hz) and very long-lived (up ∼ 1000 s) fluctuations in the firing rates (estimated in 10 s bins). It is not even clear that one can meaningfully define the “firing rate” of this neuron. Indeed, the procedure of estimating the firing rate, as above, by counting spikes in increasingly long time windows returns different values depending on the initial condition ( Fig. 4(d) ). By definition, the dynamics is ergodic when the longtime average of any (sufficiently well-behaved) observable asymptotically converges to the same value, independently of the initial condition. In our case, the observable is the number of spikes in a time window T , and its long-time average is the “firing rate”, that is, the probability of observing a spike per unit of time. As just shown, this property is satisfied for q = 0.8, and it is not satisfied for q = 0.9. The foregoing analysis, thus, reveals that the network dynamics is not ergodic at q = 0.9. We use the asymptotic convergence of firing rates to systematically study ergodicity breaking as a function of the level of symmetry, q , and of the synaptic efficacy, g . Specifically, we define a function of T, D ( T ), which quantifies the difference among the firing-rate estimates, over a time window T , obtained from different initial conditions (see Sec. IV B for details). We classify the dynamics as ergodic when D ( T ) decays as 1 /T for T → ∞ [ 43 ], and as non-ergodic otherwise. Note that changes in g at fixed h ext , the strength of the external drive, result in changes in the population-averaged firing rate, ⟨ ν ⟩. In the region of the ( q, g ) plane that we have explored numerically, ⟨ ν ⟩ ranges between about 1 Hz (for g = 7) and about 5 Hz (for g = 2). Consistent with the results illustrated in Figs. 3 and 4 , D ( T ) → 0 as 1 /T for q ≤0.8, while for q ≥ 0.9, D ( T ) does not decay, even when firing rates are estimated over time intervals as long as T = 10 4 s (see Fig. 5(a) ). In fact, as illustrated in Fig. 5(b) , ergodicity breaking occurs in a whole region of the ( q, g ) plane. As can be seen, ergodicity breaking requires large values of q , the level of symmetry in the synaptic connectivity. However, symmetry alone is not sufficient; synaptic efficacies have to be large enough. For instance, for g < 1.25, the network dynamics is ergodic at q = 1.0 ( Fig. 5(b) ). Download figure Open in new tab FIG. 5: (a) D as a function of time for different values of q : q = 0 (blue), q = 0.5 (green), q = 0.7 (yellow), q = 0.8 (orange), q = 0.9 (red), q = 1.0 (black). Five realizations (replicas) of the same network have been considered to estimate D . (b) Bidimensional phase diagram in the ( q, g ) plane, showing the ergodic and non-ergodic regimes in the red and blue regions, respectively. (c) Percentage of neurons with a firing rate lower than 0.01 ⟨ ν ⟩ (i.e., almost silent). The black dots in (b-c) indicate the parameter values for which the reported quantities in the plane have been measured. All other parameters are fixed as in Fig. 1 . The breaking of ergodicity is accompanied by the appearance of “silent” neurons. This is illustrated in Fig. 5(c) , where we report the fraction of neurons firing at less than 0.01 ⟨ ν ⟩ as a function of q and g . These neurons would be hardly detectable in experiments [ 44 ]. As can be seen, for q ≥ 0.9 and g ≥ 3, more than half of the neurons in the network are “almost silent” (up to 87% at q = 1.0). By contrast, there are no silent neurons at q = 0.7. The identity of the almost silent neurons and hence, of the active ones changes depending on the initial condition. In other words, in the region where ergodicity is broken, the network dynamics exhibits multiple, (very) long-lived states of activity differing in their spatial distribution of “active” and “silent” neurons. The fact that we were able to make such an observation by using a small number of initial conditions (i.e., 5) suggests that the network dynamics can, indeed, exhibit a very large number of such states. What is the mechanism responsible for ergodicity breaking in our model network? Generally speaking, ergodicity breaking results from the fact that the dynamics is for very long times, or even permanently, restricted to different regions of the phase space depending on the initial conditions. Perhaps the simplest scenario for the emergence of such “trapping regions” is the co-existence of multiple, stable equilibria, e.g., different fixed points of the microscopic dynamics. However, our network does not have fixed points stricto sensu , i.e., states in which the voltage of all neurons is constant in time; this would correspond to a completely silent network. To assess the relevance of the above scenario in our model network, therefore, we introduce a synaptic time constant, τ s (see Sec. IV A for details). Roughly speaking, this leads to synaptic inputs proportional to the “average”, rather than the “instantaneous”, presynaptic activity, effectively reducing their (fast) temporal fluctuations. The time window over which the presynaptic activity is averaged is controlled by τ s . The network we have studied so far corresponds to τ s = 0. In the opposite limit, i.e., τ s → ∞, the fast fluctuations due to spiking disappear, and the dynamics can be described by a rate model [ 18 , 19 ]. Note that the time-averaged input to a neuron does not depend on τ s . We start by investigating how the ergodic properties of the network dynamics are affected by τ s . To this end, we compute the ergodic distance, D ( T ), for different values of τ s and q , for g = 5. For this value of g , ergodicity is broken between q = 0.8 and q = 0.9 in the network with τ s = 0 (see Fig. 5 ). To properly compare the behavior of D ( T ) at different values of τ s > 0, we measure time in units of τ s , i.e., . If the breaking of ergodicity in our model network is associated with the appearance of stable fixed points in the “corresponding” rate model (i.e., the same network in the limit τ s → ∞), then one would expect that the region in the parameter space in which ergodicity is broken does not shrink with increasing τ s ; rather, it could increase. This is indeed what we observe. As can be seen in Fig. 6 , the dynamics remain non-ergodic for q ≥ 0.9, regardless of the value of τ s . However, now the dynamics is non-ergodic already at q = 0.8 for τ s = 0.02s, and remains so by further increasing τ s . By contrast, the dynamics remain ergodic at q = 0.7 even for τ s = 1s. We conjecture that this is because there are no stable fixed points, in the limit τ s → ∞, in the network at q = 0.7. However, it is also conceivable that, instead, stable fixed points will appear for τ s > 1s, and hence ergodicity will be broken for larger values of τ s . For practical reasons, i.e., simulations run times, we have systematically explored only τ s ≤1s. In any case, we do not expect the dynamics to become non-ergodic at q = 0 because there are no stable fixed points in the network, even in the limit τ s → ∞. Download figure Open in new tab FIG. 6: D as a function of for exponentially decaying synapses with different values of q and of the synaptic time constant: namely, τ s = 0.02 s (solid lines); 0.1 s (dashed lines); and 1 s (dotted lines). The colors refer to q = 0 (blue), q = 0.5 (green), q = 0.7 (yellow), q = 0.8 (orange), q = 0.9 (red), and q = 1.0 (black). These results were obtained using five replicas; other parameters are fixed as in Fig. 1 It is very insightful to visualize ergodic and non-ergodic dynamics in the network. For this purpose, we use the “projection” of the network dynamics into the threedimensional space ( h 1 , h 2 , h 3 ), where h i denotes the total synaptic input to neuron i . In Fig. 7 , we show these projections for five different initial conditions in the same network, with different colors denoting different initial conditions. In the top row, we illustrate the case q = 0.8 and τ s = 0 (in this case, ). As can be seen, there is no evident separation between the representative points at any time; the orbits for the different initial conditions fill “densely” the same volume in the phase space. That this is not an artifact of the low-dimensional projection is, indeed, confirmed by the ergodic distance, , that goes to 0 for . For comparison, in the bottom row, we illustrate the case q = 0 and τ s = 1 s. As is apparent, the qualitative behavior is the same as before. In fact, 0 for also in the case q = 0 and τ s = 1 s. By contrast, the behavior is qualitatively different in the case q = 0.8 and τ s = 1 s, which we illustrate in the middle row. In this case, while the orbits seem, as before, to fill the same volume at intermediate times, they asymptotically reach clearly separated regions of the phase space. Though it should be evident from the figure (see also Fig. 5 ), it is perhaps worth emphasizing the fact that such “asymptotic confinement” of the orbits— that is, reaching equilibrium—takes very long times, i.e., ∼ 10 3 τ s . Download figure Open in new tab FIG. 7: Trajectories in a 3D neural space defined by the total synaptic input of three random neurons; each point corresponds to at the same time. The different colors represent the five different replicas displayed. The upper row corresponds to instantaneous synapses with q = 0.8, while the other two refer to τ s = 1 s and q = 0.8 (middle row) and q = 0 (bottom row). The same parameters as in Fig.1 were used. III. DISCUSSION We have shown that, in a spiking network operating in the dynamically-balanced regime, the over-expression of reciprocal synaptic connections, naturally and robustly leads to long-lived fluctuations in neural activity. In turn, these fluctuations lead to long relaxation times as estimated by the decay of the autocorrelation function, consistent with experimental observations in the active cortex [ 1 – 3 ]. The emergence of this slow dynamics is independent of both the size of the network and of the synaptic time constant, unlike in previous accounts (as discussed in the Introduction). Rather, it depends on the effective excitatory self-coupling generated by the symmetry in the synaptic connectivity, as already pointed out in [ 30 ] and further elaborated upon in the present study. Accordingly, we conclude that long relaxation times are a generic feature of the balanced regime, as soon as one includes in the model a more realistic description of the local synaptic connectivity [ 32 – 36 ]. Besides confirming previous results, here we significantly expand upon them. Specifically, we have shown that the network dynamics becomes non-ergodic when the excitatory self-coupling is suitably strong. Notably, the breaking of ergodicity also occurs in presence of significant levels of fast noise (i.e., for instantaneous synaptic currents) and in partially symmetric networks. By taking long synaptic time constants, we have shown that the network dynamics becomes multi-stable when ergodicity is broken. This suggests that the mechanism underlying the ergodicity breaking in our model network could be a rather classical one, that is, the coexistence of multiple stable fixed points of the (rate) dynamics. It is important to stress that, in our model network, the multi-stability does not result from some low-rank structure embedded in the synaptic connectivity, as is the case, for instance, in model networks of associative memory [ 11 , 45 , 46 ] or, more specifically, in clustered networks [22– By contrast, the connectivity matrix in our network has no “structure” beyond a pre-assigned level of symmetry, and it is, indeed, full rank. Therefore, to the best of our knowledge, this is the first report of non-ergodic dynamics in an “unstructured” spiking network (see [ 31 ] for a related report in a network of rate neurons). The usual caveats attached to any result based on numerical simulations apply to our results as well. We have probed the robustness of our results against finitesize effects by systematically varying, within “reasonable” bounds, both the number of neurons in the network (1000 ≤ N ≤5000) and the average number of connections they make (100 ≤ K ≤3000). In all cases, we observed the same phenomenology as described in Results; in particular, we found breaking of ergodicity for q < 1 and τ s = 0, for suitably large values of g . We have also verified that the basic phenomenology does not critically depend on the specific spiking-neuron model. Concretely, we have also studied networks of quadratic integrate-and-fire (QIF) neurons. The QIF neuron is, in some respects, more realistic than the LIF neuron, since it can approximate reasonably well an en-tire class of neural models—Type I excitability according to Hodgkin classification—in proximity to the tonic firing threshold [ 47 , 48 ]. We found, essentially the same phenomenology as for the LIF model. The results for a balanced network of QIF neurons will be presented in a forthcoming publication. The best way to assess the robustness of our results would be, obviously, to develop an accompanying analytical theory. This is a formidable task for networks of LIF neurons [ 13 ]. On the other hand, using dynamical mean-field theory to achieve an analytical description of the steady states of networks, either of rate elements in presence of white noise or, even better, of spike-response models as in [ 19 ], appears feasible although very far from trivial. Taken together, our results reveal the existence of a novel phase of the collective activity in dynamically-balanced spiking networks. In this phase, the dynamics reaches, asymptotically, different equilibrium states when starting from random initial conditions. In each of these states, the network operates in the balanced regime and, accordingly, the spiking activity of the neurons is unsaturated, temporally-irregular and asynchronous. In this respect, then, it resembles the pattern of activity that one observes in random, ergodic networks. However, unlike in random networks, the ergodicity is broken in this new phase; the memory of the initial conditions never fades away and reaching equilibrium is a very slow process. An obvious question, then, is whether there is any evidence that the cortex operates in a non-ergodic phase. We are not aware of any experimental report definitively ruling out this possibility. As just discussed, many features of neural activity will be qualitatively the same in the ergodic and non-ergodic phases. Our work suggests a simple observable, easy to measure in an experiment (at least in principle), that could be help address this question. That is the ergodic distance, D ( T ), originally introduced in [ 31 ], that we have also used in our spiking network. Note that D ( T ) can be computed using only single-cell recordings. Deviations from 1 /T in its decay at long times (e.g., ∼ 10s) would be supporting evidence of ergodicity breaking, at least in a weak sense (i.e., very long transients but only one equilibrium state). At present, we are actively pursuing this line of research. IV. METHODS A. The model network We consider a network of N inhibitory LIF neurons sparsely connected so that each neuron makes and receives, K synaptic contacts on average [ 18 , 19 , 38 , 39 ]. The neurons also receive a constant excitatory drive from outside the network. The membrane potential of neuron i, v i , measured relative to the rest potential and in units of the threshold, evolves in time according to where τ m = 10 ms is the membrane time constant and h i ( t ) is the total synaptic input (external and recurrent) to neuron i . Note that the synaptic input is measured in the same units as the voltage and is thus dimensionless. Whenever v i ( t ) ≥ 1, the voltage is reset to the rest potential v i = 0, a spike is generated, and, as a result, postsynaptic currents (PSCs) are initiated at all the postsynaptic targets of neuron i . Then, the dynamics Eq. (3) resumes immediately after the spike, starting from v i = 0. The total synaptic input, h i ( t ), is given by where is the external excitatory input and E i ( t ) is given by where τ s is the synaptic time constant, g > 0 is a measure of the synaptic inhibitory strength, A ij = 1 if there exists a synaptic connection from neuron j to neuron i ( A ij = 0 otherwise), and the sum over n is over all the spike times, , of neuron j . In words, E ( t ) is the linear superposition of all the PSCs generated by the spiking activity of the (presynaptic) neurons connected to neuron i . The unitary PSC, generated by a single spike (at time t = 0), is where Θ → ( t ) is the Heaviside function. Thus, the synaptic time constant, τ s , controls the duration of the PSCs. In the limit τ s → 0, instantaneous synaptic currents, E i ( t ) becomes (see Eq. (5)) Note that synaptic currents are measured in units of 1/time (e.g., in Hz). The unitaryPSC, Eq. (6), is proportional to , so that scales with K as the external drive, see Eq. (4). To have h i ( t ) ∼ 1 when K is large (1 ≪ K ≪ N ), the recurrent synaptic input, τ m E i ( t ), should cancel (i.e., to leading order) on average the external drive, . Thus, the network operates in the dynamically-balanced regime when K is large [ 10 , 12 ]. All the simulations described in Results are with N = 4000, K = 1200 and h ext = 0.1. 1. Synaptic connectivity The synaptic connectivity matrix, A ≡ ( A ij ), is generated in the following way [ 30 ]. For each ordered pair of neurons ( i, j ), the probability, p B , that they are reciprocally connected, i.e., ( A ij = 1, A ji = 1), is given by while the probability, p U , that they are unidirectionally connected, i.e., ( A ij = 1, A ji = 0), is given by where the parameter q , 0 ≤ q ≤ 1, controls the overrepresentation of reciprocal connections, i.e., the level of symmetry. In experiments, this over-representation is often quantified in by the ratio, η , between the number of observed reciprocally-connected pairs of neurons divided by the expected number of reciprocally-connected pairs in a directed Erdős-Rényi network with the same probability of connection [ 32 , 36 ]. The probability of connection is simply estimated by dividing the number of connected pairs by the number of pairs sampled in the experiment. In our case, Hence, η = 1 for q = 0 and the network is simply a directed Erdős-Rényi network, while η = N/K for q = 1 and all connections are reciprocal, i.e., the network is fully symmetric. We note that, while q is independent of the probability of connection (i.e, K/N ), η is not. In particular, the largest possible value of η is the reciprocal of the probability of connection. Hence, larger probabilities of connection necessarily result in lower values of η , regardless of the actual level of symmetry, measured more properly by q . For instance, in the cortex local excitatory-to-excitatory connections occur with probability ∼0.1 and η ≃ 5 [ 36 ]; this leads to q ≃ 0.5 using Eq. (10) (taking, obviously, the solution with 0 ≤ q ≤1). Inhibitory-to-inhibitory connections, instead, occur with probability ∼0.3 and η ≃ 3 [ 36 ]; this leads to q ≃ 0.9. Hence, though η is larger for the excitatory-to-excitatory network, the inhibitory-to-inhibitory network is more “symmetric”. Indeed, about 90% of its connections are bidirectional against only about 50% in the excitatory-to-excitatory network. 2. Network simulations We have employed event-driven techniques to simulate the network dynamics. This allowed us to perform fast and accurate simulations of large networks. In fact, the differential equations, Eqs.(3)-(5), describing the network dynamics are exactly integrated between consecutive spikes. This can be done by rewriting the evolution of the v i ’s and of the E i ’s as a mapping between their value at a spiking event occurring at time t n and their value at the next spiking event at t n +1 [ 49 ]. The identity of the neuron firing at the next spiking event and the time, t n +1 , at which the event occurs can be determined in the following way. Suppose that neuron k has fired at the current spiking event, at time t n . Then, v k ( t + ) = 0 (as a result of the resetting) and where is the time immediately after/before the spike. After a time interval Δ, if no other spike occurs, we have and The last term on the r.h.s. is Using Eqs. (13) and (14), we can determine the next spiking time, relative to t n , for each neuron i if no other spike occurs before, i.e., Δ i such that for each i . Then, the next actual spiking event occurs at time and the identity of the neuron spiking at t n +1 , m , is simply that is, the neuron i for which Δ i achieves its minimum value. In the case of instantaneous synapses, we use Eq. (7) into Eq. (13), and the last term on the r.h.s. of Eq. (13) is now B. Observables 1. Single-neuron firing activity The single-neuron activity has been characterized in terms of the following observables: The firing rate For each neuron i , we obtained its average firing rate, ν i , by counting all the spikes that the neuron fired during the simulation and then dividing the spike count by the total simulation time, T . The population-averaged firing rate, ⟨ ν ⟩, is given by Estimates of the ν i ’s were obtained by using simulations with T = 1000s or longer. In some cases (see, e.g., Fig. 4 ), we also estimated the instantaneous firing rate of neuron i, r i , that is, the firing rate in a given time window of duration Δ < T . This was obtained by counting only the spikes that the neuron fired within the time window and then dividing this spike count by Δ. The coefficient of variation The coefficient of variation of neuron i, CV i , is the ratio between the standard deviation of its inter-spike intervals and the average of its interspike intervals [ 4 ]. CV i provides a simple quantification of the degree of temporal irregularity of the spiking process. If neuron i spikes periodically, then CV i = 0; if neuron i spikes according to a Poisson process, then CV i = 1. The average and standard deviation of the inter-spike intervals have been estimated with simulations of T = 1000s. 2. Pairwise spike-count correlations We choose the bin size Δ = 100 ms and divide the interval [0, T ], where T is the total duration of the simulation, into T/ Δ non-overlapping bins. Then, for each neuron i we compute where k = 1, …, T/ Δ is the bin number, ν i is the average firing rate of neuron i estimated in [0, T ], and the sum over a is over all the spike times, , of neuron i ; ν Δ is the expected spike count in a time bin Δ. Hence, for each i = 1, …, N . Then, for each pair of neurons ( i, j ), ρ ij is computed as and it lies between −1 and 1. 3. Autocorrelation of the synaptic input For each neuron i , we compute where T is the total duration of the simulation. In the case of instantaneous synapses, we low-pass filter h i ( t ) with a cutoff frequency f c = 50Hz before using Eq. (22). Hence, for each i = 1, …, N . We define the normalized autocorrelation, C i ( τ ), of the input to neuron i at a time lag τ as so that C i (0) = 1 for each i = 1, …, N . Finally, 4. Spike-triggered input For each neuron i and lag τ , we define Z i ( τ ) as the set of spike numbers, n , for which if τ > 0, or the n ’s for which if τ < 0. Thus, if n ∈ Z ( τ ), there are no spikes between and . The spike-triggered average input to neuron i at lag τ , STA i ( τ ), is given by where | Z i ( τ ) | is the number of elements in Z i ( τ ), the sum is restricted to elements in Z i ( τ ) and denotes the average of h i ( t ) over time, i.e., We low-pass filtered h i ( t ) with a cutoff frequency f c = 100 Hz before computing Eqs. (26) and (27). By averaging over neurons, we finally obtain the spike-triggered average input, STA ( τ ), shown in Fig. 2 , i.e., 5. Ergodic “distance” To investigate the ergodicity of the network dynamics, we consider M different physical replicas of the same network. Each replica is the time evolution of the network starting from different initial conditions. Then, for two different replicas, α and β , we compute where is the spike count of the neuron i up to time T in replica α and is the corresponding firing rate. By averaging over all the distinct pairs of replicas, we finally obtain D ( T ) [ 31 ], that is, ACKNOWLEDGEMENT We acknowledge stimulating interactions with Alberto Bacci, Alejandro Carballosa, and Nina La Miciotta. G.M. would like to thank Stefano Fusi for a useful discussion about slow changes in excitability at the singleneuron level. A. T. received financial support by the Labex MME-DII (Grant No. ANR-11-LBX-0023-01) (together with F.T.) and by CY Generations (Grant No ANR-21-EXES-0008) (together with G.M.), all part of the French program Investissements d’Avenir. Funder Information Declared Agence Nationale de la Recherche, https://ror.org/00rbzpz17 , ANR-11-LBX-0023-01 , ANR-21-EXES-0008 References [1]. ↵ T. Ogawa and H. Komatsu , Journal of neurophysiology 103 , 2433 ( 2010 ). OpenUrl CrossRef PubMed Web of Science [2]. J. D. Murray , A. Bernacchia , D. J. Freedman , R. Romo , J. D. Wallis , X. Cai , C. Padoa-Schioppa , T. Pasternak , H. Seo , D. Lee , et al. , Nature neuroscience 17 , 1661 ( 2014 ). OpenUrl CrossRef PubMed [3]. ↵ S. Nishida , T. Tanaka , T. Shibata , K. Ikeda , T. Aso , and T. Ogawa , Cerebral Cortex 24 , 1671 ( 2014 ). OpenUrl CrossRef PubMed Web of Science [4]. ↵ P. Dayan and L. Abott , “Theoretical neuroscience,” ( 2001 ). [5]. ↵ G. La Camera , A. Rauch , D. Thurbon , H.-R. Luscher , W. Senn , and S. Fusi , Journal of neurophysiology 96 , 3448 ( 2006 ). OpenUrl CrossRef PubMed Web of Science [6]. ↵ C. Pozzorini , R. Naud , S. Mensi , and W. Gerstner , Nature neuroscience 16 , 942 ( 2013 ). OpenUrl CrossRef PubMed [7]. ↵ S. A. Fisher , T. M. Fischer , and T. J. Carew , Trends in neurosciences 20 , 170 ( 1997 ). OpenUrl CrossRef PubMed Web of Science [8]. ↵ M. Beiran and S. Ostojic , PLoS computational biology 15 , e1006893 ( 2019 ). OpenUrl [9]. ↵ M. V. Tsodyks and T. Sejnowski , Network: Computation in Neural Systems 6 , 111 ( 1995 ). OpenUrl CrossRef [10]. ↵ C. Van Vreeswijk and H. Sompolinsky , Science 274 , 1724 ( 1996 ). OpenUrl Abstract / FREE Full Text [11]. ↵ D. J. Amit and N. Brunel , Cerebral cortex (New York, NY: 1991) 7 , 237 ( 1997 ). OpenUrl CrossRef PubMed Web of Science [12]. ↵ A. Renart , J. De La Rocha , P. Bartho , L. Hollender , N. Parga , A. Reyes , and K. D. Harris , science 327 , 587 ( 2010 ). OpenUrl Abstract / FREE Full Text [13]. ↵ A. Lerchner , C. Ursta , J. Hertz , M. Ahmadi , P. Ruffiot , and S. Enemark , Neural computation 18 , 634 ( 2006 ). OpenUrl CrossRef PubMed Web of Science [14]. ↵ A. Roxin , N. Brunel , D. Hansel , G. Mongillo , and C. van Vreeswijk , Journal of Neuroscience 31 , 16217 ( 2011 ). OpenUrl Abstract / FREE Full Text [15]. G. Mongillo , D. Hansel , and C. Van Vreeswijk , Physical review letters 108 , 158101 ( 2012 ). OpenUrl CrossRef PubMed [16]. ↵ G. Mongillo , S. Rumpel , and Y. Loewenstein , Nature neuroscience 21 , 1463 ( 2018 ). OpenUrl CrossRef PubMed [17]. ↵ M. Helias , T. Tetzlaff , and M. Diesmann , PLoS computational biology 10 , e1003428 ( 2014 ). OpenUrl PubMed [18]. ↵ O. Harish and D. Hansel , PLoS computational biology 11 , e1004266 ( 2015 ). OpenUrl PubMed [19]. ↵ J. Kadmon and H. Sompolinsky , Physical Review X 5 , 041030 ( 2015 ). OpenUrl [20]. ↵ D. Dahmen , S. Grün , M. Diesmann , and M. Helias , Proceedings of the National Academy of Sciences 116 , 13051 ( 2019 ). OpenUrl Abstract / FREE Full Text [21]. ↵ D. Angulo-Garcia , S. Luccioli , S. Olmi , and A. Torcini , New Journal of Physics 19 , 053011 ( 2017 ). OpenUrl [22]. ↵ G. Deco and E. Hugues , PLoS computational biology 8 , e1002395 ( 2012 ). OpenUrl PubMed [23]. ↵ A. Litwin-Kumar and B. Doiron , Nature neuroscience 15 , 1498 ( 2012 ). OpenUrl CrossRef PubMed [24]. ↵ L. Mazzucato , A. Fontanini , and G. La Camera , Journal of Neuroscience 35 , 8214 ( 2015 ). OpenUrl Abstract / FREE Full Text [25]. ↵ B. A. Brinkman , H. Yan , A. Maffei , I. M. Park , A. Fontanini , J. Wang , and G. La Camera , Applied Physics Reviews 9 ( 2022 ). [26]. ↵ C. Huang and B. Doiron , Current opinion in neurobiology 46 , 31 ( 2017 ). OpenUrl CrossRef PubMed [27]. ↵ P. Hänggi , P. Talkner , and M. Borkovec , Reviews of modern physics 62 , 251 ( 1990 ). OpenUrl CrossRef Web of Science [28]. ↵ X. Yang and G. La Camera , PLOS Computational Biology 20 , e1012220 ( 2024 ). OpenUrl PubMed [29]. ↵ D. Martí , N. Brunel , and S. Ostojic , Physical Review E 97 , 062314 ( 2018 ). OpenUrl CrossRef [30]. ↵ S. Rao , D. Hansel , and C. van Vreeswijk , Scientific Reports 9 , 3334 ( 2019 ). OpenUrl PubMed [31]. ↵ K. Berlemont and G. Mongillo , bioRxiv, 2022 ( 2022 ). [32]. ↵ S. Song , P.J. Sjöström , M. Reigl , S. Nelson , and D. B. Chklovskii , PLoS biology 3 , e68 ( 2005 ). OpenUrl CrossRef PubMed [33]. Y. Yoshimura , J. L. Dantzker , and E. M. Callaway , Nature 433 , 868 ( 2005 ). OpenUrl CrossRef PubMed Web of Science [34]. H. Ko , S. B. Hofer , B. Pichler , K. A. Buchanan , P.J. Sjöström , and T. D. Mrsic-Flogel , Nature 473 , 87 ( 2011 ). OpenUrl CrossRef PubMed Web of Science [35]. N. Brunel , Nature neuroscience 19 , 749 ( 2016 ). OpenUrl CrossRef PubMed [36]. ↵ L. Campagnola , S. C. Seeman , T. Chartrand , L. Kim , A. Hoggarth , C. Gamlin , S. Ito , J. Trinh , P. Davoudian , C. Radaelli , et al. , Science 375 , eabj5861 ( 2022 ). OpenUrl CrossRef PubMed [37]. ↵ A. Crisanti and H. Sompolinsky , Physical Review A 37 , 4865 ( 1988 ). OpenUrl PubMed [38]. ↵ M. Monteforte and F. Wolf , Physical review letters 105 , 268104 ( 2010 ). OpenUrl CrossRef PubMed [39]. ↵ M. Di Volo and A. Torcini , Physical review letters 121 , 128301 ( 2018 ). OpenUrl CrossRef PubMed [40]. ↵ R. Capocelli and L. Ricciardi , Kybernetik 8 , 214 ( 1971 ). OpenUrl CrossRef PubMed Web of Science [41]. ↵ We used the function S(x) = 25/(1 + e-k(x-x0)) to approximate the data. The two free parameters, x0 and k, were determined by least-squares fitting. We obtained: x0 = 0.32 and k = 5.67 for q = 0, x0 = 0.29 and k = 5.40 for q = 0.25, x0 = 0.22 and k = 4.94 for q = 0.5 . [42]. ↵ Clearly, since STA(t) is a rate, the integral in time R is dimensionless . [43]. ↵ C. W. Gardiner et al. , Handbook of stochastic methods , Vol. 3 ( springer Berlin , 1985 ). [44]. ↵ S. Shoham , D. H. O’Connor , and R. Segev , Journal of Comparative Physiology A 192 , 777 ( 2006 ). OpenUrl CrossRef PubMed Web of Science [45]. ↵ J. J. Hopfield , Proceedings of the national academy of sciences 79 , 2554 ( 1982 ). OpenUrl Abstract / FREE Full Text [46]. ↵ E. Curti , G. Mongillo , G. La Camera , and D. J. Amit , Neural computation 16 , 2597 ( 2004 ). OpenUrl CrossRef PubMed Web of Science [47]. ↵ W. Gerstner , W. M. Kistler , R. Naud , and L. Paninski , Neuronal dynamics: From single neurons to networks and models of cognition ( Cambridge University Press , 2014 ). [48]. ↵ B. Gutkin , in Encyclopedia of computational neuroscience ( Springer , 2022 ) pp. 3412 – 3419 . [49]. ↵ R. Zillmer , R. Livi , A. Politi , and A. Torcini , Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 76 , 046102 ( 2007 ). OpenUrl View the discussion thread. Back to top Previous Next Posted July 31, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Is the cortical dynamics ergodic? A numerical study in partially symmetric networks of spiking neurons Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Is the cortical dynamics ergodic? A numerical study in partially symmetric networks of spiking neurons Ferdinand Tixidre , Gianluigi Mongillo , Alessandro Torcini bioRxiv 2025.07.28.666831; doi: https://doi.org/10.1101/2025.07.28.666831 Share This Article: Copy Citation Tools Is the cortical dynamics ergodic? A numerical study in partially symmetric networks of spiking neurons Ferdinand Tixidre , Gianluigi Mongillo , Alessandro Torcini bioRxiv 2025.07.28.666831; doi: https://doi.org/10.1101/2025.07.28.666831 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Neuroscience Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17691) Bioengineering (13892) Bioinformatics (41936) Biophysics (21452) Cancer Biology (18588) Cell Biology (25504) Clinical Trials (138) Developmental Biology (13378) Ecology (19899) Epidemiology (2067) Evolutionary Biology (24320) Genetics (15609) Genomics (22506) Immunology (17736) Microbiology (40394) Molecular Biology (17181) Neuroscience (88605) Paleontology (666) Pathology (2832) Pharmacology and Toxicology (4824) Physiology (7641) Plant Biology (15153) Scientific Communication and Education (2045) Synthetic Biology (4294) Systems Biology (9825) Zoology (2271)