Predicting brain evoked response to external stimuli from temporal correlations of spontaneous activity
aa r X i v : . [ q - b i o . N C ] S e p Predicting brain evoked response to external stimuli from temporal correlations ofspontaneous activity
A. Sarracino, ∗ O. Arviv,
2, 3, † O. Shriki,
2, 3, 4, ‡ and L. de Arcangelis § Engineering Department, University of Campania “Luigi Vanvitelli”, 81031 Aversa (Caserta), Italy Department of Cognitive and Brain Sciences, Ben-Gurion University of the Negev, Beer-Sheva, Israel The Inter-Faculty School for Brain Sciences, Zlotowski Center for Neuroscience,Ben-Gurion University of the Negev, Beer-Sheva, Israel Department of Computer Science, Ben-Gurion University of the Negev, Beer-Sheva, Israel
The relation between spontaneous and stimulated global brain activity is a fundamental prob-lem in the understanding of brain functions. This question is investigated both theoretically andexperimentally within the context of nonequilibrium fluctuation-dissipation relations. We considerthe stochastic coarse-grained Wilson-Cowan model in the linear noise approximation and compareanalytical results to experimental data from magnetoencephalography (MEG) of human brain. Theshort time behavior of the autocorrelation function for spontaneous activity is characterized by adouble-exponential decay, with two characteristic times, differing by two orders of magnitude. Con-versely, the response function exhibits a single exponential decay in agreement with experimentaldata for evoked activity under visual stimulation. Results suggest that the brain response to weakexternal stimuli can be predicted from the observation of spontaneous activity and pave the way tocontrolled experiments on the brain response under different external perturbations.
I. INTRODUCTION
The brain represents one of the most fascinating sys-tems where several mechanisms at different scales aredeeply intertwined, resulting in a complex behavior.One of the main open issues in the understanding ofbrain functioning is the relation between spontaneousand evoked activity, namely the response of the systemto external stimuli. In particular, it has been found thatthe large variability in the response to repeated presen-tations of the same stimulus can be attributed to theongoing spontaneous activity [1–3]. However, a theoret-ical framework to formalize this question is still lackingand the quantitative connection between the spontaneousand evoked activity remains unclear [4, 5].Experimental results for temporal correlations of spon-taneous brain activity have been reported in a num-ber of studies. For instance, the seminal article [6] fo-cused on correlations for spontaneous alpha oscillationsin the healthy human brain, and found a decay charac-terized by power-law behavior at long times. More recentanalyses [7] measured autocorrelations of spiking activ-ity fluctuations in several cortical areas of the macaquemonkey. For each area an intrinsic timescale is definedfrom the exponential fit of the autocorrelation, which de-cays to a nonzero offset value taking into account longertimescales. A hierarchical ordering of intrinsic charac-teristic times across areas was evidenced in [8]. In par-ticular, a hierarchy of time scales coupling the dynamicsof sensory regions at different ordering levels, was found ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] in the auditory functional area [9]. The observed rela-tion between time-scales suggests a temporal organiza-tion across the cerebral cortex. Other studies on thecortex dynamics, focusing on the effects of sleep depriva-tion on time correlations, observed exponential decays:Sustained wakefulness appears to decrease correlationtimescales in humans [10], and cortical timescales dependon time awake and sleep stage in rats [11]. Even if spon-taneous activity has been found to be highly structuredand to participate in high cognitive functions, its func-tional role remains poorly understood. In particular, it iswell accepted that variability is reduced during tasks atdifferent scales, yet there is no general consensus whethertask activation would inflate or reduce functional connec-tivity and correlations [12, 13]. The question we addresshere is to provide a theoretical formulation for the rela-tion between the response and the correlation function ofrest activity by means of a simple neuronal model suit-able for an analytical approach.The relation between ongoing and stimulated activitycan be addressed theoretically within the general frame-work of statistical physics, by means of the fluctuation-dissipation relations, connecting the spontaneous fluctu-ations of a system with the response function to externalperturbations. Recently, these relations have been ex-tended beyond the context of standard thermodynamics,to nonequilibrium systems, such as active and biologicalmatter [14–19]. These relations hold for a wide varietyof physical systems, not necessarily at the critical point.Their deep predictive power relies on the possibility to es-timate the response of a system from the observation ofits spontaneous fluctuations. However, in order to writeexplicit relations, a theoretical model describing the sys-tem dynamics at the scale of interest is needed. To thisextent, a useful approach to describe the activity of alarge number of neurons is provided by the coarse-grainedWilson-Cowan model (WCM) [20, 21]. The stochasticversion of this model can describe the fluctuations of ac-tive excitatory and inhibitory neuron populations via twoLangevin equations, whose dynamics is coupled througha feed-forward term, and can reproduce the statistics ofburst activity [22]. Recent studies focused on dynamicalstability properties of this model revealing a rich dynam-ical behavior, see for instance [23–25].In this work, we address the fundamental problem ofunveiling the quantitative relation between spontaneousand stimulated brain activity. We compare the analyt-ical expression for correlations and response function inthe WCM with spontaneous and stimulated brain activ-ity measured via magnetoencephalography (MEG). Wemeasure the time correlation functions of spontaneousactivity in several healthy subjects, finding a temporalbehavior mainly characterized by a double-exponentialdecay. This two-timescale decay observed in experimentsis in good agreement with the analytical prediction of theWCM. We also calculate, in the linear response regime,the response function to an external stimulation and com-pare the result to MEG data of evoked activity, obtainedby applying a visual stimulation (pictures of faces) to thesame subjects. Our study enlightens how some proper-ties of the induced brain dynamics can be predicted fromthe observation of spontaneous activity in the absence ofstimuli.The paper is organized as follows. In Sec. II we brieflyrecall the definition of the WCM, focusing on its lin-earized version, which is expected to hold in the largesize limit. Then, in Sec. III, we report the experimen-tal results for the spontaneous activity autocorrelationfunctions from MEG data. In Sec. IV we derive thefluctuation-dissipation relation for the WCM and com-pare its predictions with experimental data for evokedactivity. Finally some conclusions are drawn in the lastsection. Appendices A, B, and C report details on theexperimental procedure and on analytical calculations. II. LINEARIZED WILSON-COWAN MODEL
In order to provide a theoretical framework to under-stand the relation between spontaneous and evoked brainactivity, we consider the linearized version of the stochas-tic coarse-grained WCM, which allows us to describe thesystems dynamics at meso- and macroscales [21]. Thismodel describes the excitatory and inhibitory neuronpopulation dynamics, which evolve according to a MasterEquation for the probability p k,l ( t ) to have k excitatoryand l inhibitory neurons active at time t [22]. The totalnumber of neurons is N E + N I , where N E and N I are thepopulations of excitatory and inhibitory neurons, respec-tively. Assuming N E = N I = N , for large N the numberof active neurons is written as the sum of a deterministiccomponent ( E, I ) = ( k/N, l/N ), and a stochastic fluctu-ation ( ξ E , ξ I ), scaled by √ N , i.e. k = N E + √ N ξ E and l = N I + √ N ξ I . The Master Equation can be rewrittenin terms of these new variables and expanded in power of N obtaining a Fokker-Planck equation. This leads to thedynamic equations for the deterministic and stochasticcomponents. Next, introducing the global variables [22]Σ = ( E + I ) / , ∆ = ( E − I ) / , (1)in the large N limit and close to the fixed point (Σ , ∆ =0) the deterministic components satisfy the equations d Σ dt = − α Σ+(1 − Σ) f ( s ) /τ , d ∆ dt = − ∆( α + f ( s ) /τ ) , (2)where α is the decay rate of the active state, f ( s ) =tanh( s ) for s > τ = 1 ms is a fixedmicroscopic time-scale and s = w Σ + ( w E + w I )∆ + h .Here w E ( w I ) is the synaptic strength of the excitatory(inhibitory) population, w = w E − w I and h a smallexternal input. We notice that Σ represents the globalneuronal activity, whereas ∆ measures the unbalance inactivity between the excitatory and the inhibitory popu-lation. The deterministic equations have the unique sta-ble solution (Σ , = 0represents the condition of balance between excitationand inhibition, as found for spontaneous activity of neu-ronal systems in healthy state [28].Concerning the fluctuations of Σ and ∆, ξ Σ and ξ ∆ , inthe large N limit (the linear noise approximation) andclose to the fixed point (Σ , ∆ = 0), they satisfy [22] ddt (cid:18) ξ Σ ξ ∆ (cid:19) = (cid:18) − /τ w ff − /τ (cid:19) (cid:18) ξ Σ ξ ∆ (cid:19) + p α Σ (cid:18) η Σ η ∆ (cid:19) , (3)where η Σ and η ∆ are zero average, delta-correlated whitenoises with unitary variance, and 1 /τ = α + f ( s ) /τ − (1 − Σ ) w f ′ ( s ) /τ , w ff = (1 − Σ )( w E + w I ) f ′ ( s ) /τ ,1 /τ = α + f ( s ) /τ , with s = w Σ + h . The off-diagonal term w ff is called hidden feed-forward term andplays a central role because it rules the coupling betweenthe global activity and the unbalance between the ac-tivities of excitatory and inhibitory populations [29]. Inparticular, a fluctuation in ∆ affects the temporal evolu-tion of Σ but not vice-versa.Eqs. (3) are two coupled linear Langevin equations, thedynamics of which can be easily solved analytically [30].In particular, the solution in vectorial form reads ξ ( t ) = e M t ξ (0) + p α Σ Z t e M ( t − t ′ ) η ( t ′ ) dt ′ , (4)where M is the coupling matrix appearing in Eq. (3).Therefore, we evaluate the correlation matrix at station-arity C ij ( t ) ≡ h ξ i ( t ) ξ j (0) i = ( e M t σ ) ij , (5)where h· · · i denote average over noise, i, j = (Σ , ∆) and σ is the covariance matrix which is a known functionof the model parameters (see Appendix C). The specificcase of the autocorrelation for the total activity Σ, C ΣΣ ≡ time (ms) C ( t ) / C ( ) Average correlation (global activity)Average correlation (occipital area)Average correlation (temporal area)fit: A exp(-t/ τ ) + (1-A) exp(-t/ τ ) FIG. 1. Normalized time autocorrelation of the spontaneousactivity. in linear (main) and semi-log (inset) scale for globalactivity and for the two functional areas. Lines representthe average over 14 samples (2 per subject). The blue thickdashed line represents the best fit curve f ( t ) = A exp( − t/τ )+(1 − A ) exp( − t/τ ), with parameters A = 0 . ± . τ =8 . ± . τ = 515 ±
200 ms. h ξ Σ ( t ) ξ Σ (0) i , takes the form C ΣΣ ( t ) = α Σ τ − − τ − ) (6) × τ − − τ − − w ff τ − e − t/τ + w ff τ − e − t/τ ! . Therefore, the behavior of the correlations of the ac-tive neuron population is characterized by a double-exponential decay, with two characteristic times τ and τ , which are a function of the model parameters as speci-fied above. Similar double exponential decay is also foundfor the C Σ∆ ( t ), whereas a single exponential function withcharacteristic time τ is obtained for C ∆Σ ( t ) and C ∆∆ ( t )(see Appendix C). III. TIME CORRELATION OF SPONTANEOUSACTIVITY FROM MEG DATA
The variable Σ describes the brain spontaneous activ-ity and therefore allows for a comparison with experimen-tal data obtained from MEG measurements. We analysethe global signal, by summing the signals at all sensors,the signal from sensors monitoring the visual area in theoccipital lobe and the signal from sensors monitoring thetemporal lobes, which include the area specialized in facerecognition, e.g., the Fusiform Face Area (FFA). Sponta-neous activity was recorded from 7 healthy human sub-jects with a whole-head 248-channel magnetometer array(4-D Neuroimaging, Magnes 3600WH) in the MEG facil-ity at the EMBI Unit, Bar-Ilan University, Israel. TheMEG signal was recorded at a sampling rate of 1017.25Hz and offline band-pass filtered between 0.8 and 80 Hzas well as underwent a cleaning procedure of potentialartefacts. The absolute amplitudes were summed acrossthe recorded channels to give a global rest activity signal. More details can be found in Appendix A and in [26, 27].We compute the normalized autocorrelation function C ( t ) C (0) = P N − ts =1 ( | x ( s ) | − µ )( | x ( s + t ) | − µ ) σ (7)from time series (of total length about 4 minutes foreach subject) of the absolute value {| x ( t ) |} of rest ac-tivity signal, with mean µ and standard deviation σ .We focus here on the short time behavior (up to one sec-ond), where the experimental data show a decay mainlycharacterized by two typical times (see Fig. 1). In thisstudy we neglect oscillations at very short time scales,which probably correspond to α brain rhythm, and con-sider the global decay behavior. Analyzing 7 subjects (2data sets for each subject), we observe that this com-plex behavior is quite stable (see Table in Appendix B)and can be fitted by Eq. (6). Moreover, the same func-tional behavior is found for the global signal and for thesignal averaged only on sensors from temporal and occip-ital lobes, corresponding to areas involved in processingof the stimuli. The double exponential behavior is alsodetected for data at single sensors in these two areas,where the autocorrelation function, as expected, exhibitslarger statistical fluctuations (see Appendix B). The av-erage values of the characteristic times are found to be τ = 8 . ± . τ = 515 ±
200 ms, differing byalmost two orders of magnitude. The short time scale,causing an abrupt decay in the correlation function, is ofthe order of magnitude of synaptic time scales as well assingle neuron time scales, such as the refractory period[31]. Conversely, the long characteristic time is compati-ble with the time scale of low-frequency rhythms, as the θ rhythm. Interestingly, the long characteristic time τ controls the single exponential decay of the correlations C ∆Σ ( t ) and C ∆∆ ( t ) (see Appendix C). This suggests thata fluctuation in Σ or ∆ affects the temporal behaviorof ∆ itself over a long timescale, implying a slow recov-ery of the balance condition if the system is moved awayfrom the fixed point ∆ = 0. From the fitting procedurewe also obtain an estimation for the feed-forward coeffi-cient w ff = 0 . ± . − , which characterizes thecoupling between excitatory and inhibitory neuron pop-ulations. Finally, we note that, at longer times ( t & IV. PREDICTING RESPONSE FROMSPONTANEOUS ACTIVITY VIA THEFLUCTUATION-DISSIPATION THEOREM
The present theoretical framework allows us to addressthe fundamental question of the relation between spon-taneous and evoked activity. We can evaluate the linearresponse function to a weak external perturbation, de-fined as R ij ( t ) ≡ δξ i ( t ) δξ j (0) , (8)where i, j = (Σ , ∆). This represents the average responseof the variable ξ i ( t ) at time t to an impulsive perturbationapplied to the variable ξ j at time 0, and · denotes non-stationary averages over many realizations. From Eq.(4)one immediately obtains the response matrix R ( t ) = e M t for t > . (9)From Eq. (5), we obtain a fluctuation-dissipation relation R ( t ) = C ( t ) σ − , (10)connecting the linear response with spontaneous fluctu-ations [14]. In particular, since we are interested in theresponse of the total activity Σ to a small perturbationof Σ itself, the response function reads R ΣΣ ( t ) = ( σ − ) h ξ Σ ( t ) ξ Σ (0) i + ( σ − ) h ξ Σ ( t ) ξ ∆ (0) i , (11)where the explicit form of the inverse matrix σ − is re-ported in Appendix C. The relevance of relation (11) re-lies on the fact that it allows us to get information onthe response of the system to an external perturbationby simply looking at its unperturbed dynamics. In gen-eral, as evident from Eq. (11), both the autocorrelationand the cross-correlation are required to reconstruct theresponse to a weak external stimulus. This is due tothe presence of nonequilibrium conditions breaking de-tailed balance, as noted in several systems (see for in-stance [14, 32, 33]).Let us note however, that, in our particular case, dueto the upper triangular form of the matrix M in Eq. (3),the response function of the model takes the simple ex-ponential form R ΣΣ ( t ) = exp( − t/τ ) , (12)where the characteristic time τ is the same ruling theshort-time behavior of the correlation function. This isconsistent with Eq. (11) because one can easily checkthat the cross-correlation appearing in Eq. (11) exactlycancels the term ∝ exp( − t/τ ) in Eq. (6), giving thesingle-exponential decay for the response function (seeAppendix C). Therefore, in this approach, measurementsof the spontaneous fluctuations in the global brain activ-ity Σ alone could provide a prediction for the systemresponse. Concerning the other response functions, wenotice that R Σ∆ ( t ) shows a double exponential behav-ior, whereas, as expected from the triangular form of M , R ∆Σ ( t ) = 0 and R ∆∆ ( t ) = exp( − t/τ ) (see AppendixC).The response function derived by the analytical cal-culation characterizes the behavior close-in-time to thestimulus application. In order to fairly compare the ana-lytical prediction to experimental data, the experimental time (ms) R ( t ) Average response (global activity)Average response (occipital area)Average response (temporal area)fit: exp(-t/ τ R ) FIG. 2. (Rescaled and shifted) average response function forthe global signal, the signal from the occipital and the tempo-ral lobes. The typical time for the response function from anexponential fit decay, according to Eq. (12), is τ R = 55 ± perturbation has to be sufficiently weak to remain in thelinear regime. Face processing is one of the most stud-ied cognitive abilities. Humans are considered expertsin face processing, which accordingly does not involvelong-lasting cognitive engagement, potentially mimickinga small perturbation to the system. The evoked responseis fast and of relatively short (typically characterized by100, 170 and 250 ms evoked response fields (ERFs), andsome late components at 400 and 600 ms) and the re-turn to baseline in both power and spectral content oc-curs prior to the 1 sec limit [26, 27]. We measured thebrain activity of the same 7 participants evoked by visualstimuli, by showing a series of pictures of human faces.Each stimulus is presented for 1 sec, with varying interstimulus intervals (larger than 1 sec). Further details arereported in Appendix A.In Fig. 2 we show the response function averaged overthe 7 subjects for the global signal and for the signal av-eraged only on sensors from temporal and occipital lobes,corresponding to areas involved in processing of the stim-uli (see Appendix B for the response function for eachsubject and at single sensor). Since each subject showsa different latency activity associated to each stimuluswhich cannot be accounted for by analytical predictions,experimental data for the response have been shifted intime and rescaled in amplitude in order to have that themaximum in each dataset is reached at time zero and itsvalue is normalized to 1. We have to stress that the ap-plied perturbation cannot be considered impulsive, as inthe analytical definition, and moreover the system takesa certain time to develop the response signal, as shownin Appendix B. However, experimental data on evokedactivity qualitatively confirm the single exponential de-cay predicted by the WCM for all datasets. Indeed, im-posing a fit with a double exponential leads in all casesto an amplitude close to zero for the second exponen-tial decay. The exponential behavior is also detected fordata at single sensors, where the response function, asexpected, exhibits larger statistical fluctuations (see Ap-pendix B). The averaged fitted value for the characteris-tic time τ R = 55 ± τ (see Table in Appendix B for the char-acteristic time values in each dataset). The quantitativedifference between these two timescales could be due toa number of factors, as the properties of the stimula-tion, which in the theoretical approach is supposed to beinstantaneous and small. The real stimulation requiresa time delay before the decay sets in, raising the needto identify the initial time to fit the response function.Fluctuations in experimental conditions then provide awider range of characteristic times that, averaged overtrials, lead to a larger τ R . Random fluctuations also hidea possible correlation between τ and τ R across subjects.Response experiments under different stimulation proto-cols and controlled conditions are required to test morequantitatively the timescale agreement. The robust func-tional behavior found for global signals, functional ar-eas signals and data at single sensor, is an indication ofscale-invariance, namely the correlation and the responsefunction do not depend on the sample size, therefore ev-idencing the self-similar properties of the system. Thisobservation could be the outcome of the coarse-grainedmeasure of the activity by MEG, since sensors are spacedat ∼ t h W t i = 2 w ff α Σ ( τ − + τ − ) . (13)This makes clear that the nonequilibrium source in thismodel is the presence of the feed-forward term w ff , whichcouples the fluctuations in the two variables Σ and ∆. V. CONCLUSIONS
We have proposed a theoretical approach to predictthe relation between global spontaneous and evokedbrain activity, based on the linear noise approximationof the WCM, which allows for analytical calculations.The theoretical approach relies on two main assump-tions: The model is two-dimensional and parametersare tuned to achieve balanced excitation and inhibition,which leads to the fixed point ∆ = 0. The comparisonwith MEG data confirms that temporal correlations ofspontaneous fluctuations are mainly characterized by adouble-exponential decay, with two well separated typi-cal times τ ≃
10 ms and τ ≃
500 ms. The short time scale is related to synaptic timescales, while the long oneis compatible with slow brain rhythms. This analyticalapproach is therefore able to rationalize experimental re-sults for the correlation and response functions providinga coherent framework for the dual functional behaviorfound experimentally. An important aspect of our re-sults is that the functional behavior found in MEG datais scale-invariant, with stable value of the characteristictimes across scales, from the cm scale (single sensor data)to the entire brain (global signal). This observation sug-gests that at first approximation the linear model wellaccounts for the brain behavior at different scales, evenif the connectome has a complex modular structure.The presence of a single exponential decay with theshort characteristic time for the response function is fullycoherent with the efficient performance of the humanbrain in visual tasks. Indeed, it is well known that thehuman eye can appreciate about 10 images per second,resulting in a temporal resolution of about 100 ms [36]. Alarger relaxation time, of the order of τ ≃
500 ms, wouldimply an overlap in the response to close-in-time stimu-lations, affecting the performance. We deem that futureexperiments on evoked activity, properly designed to theapplication of the fluctuation-dissipation theorem, couldshed new light on the brain functionality at large scale,with strong impact in neurobiology and neuroscience.Moreover, this result opens the way to numerical studiesimplementing microscopic models, as integrate and fireneuronal models, on complex networks able to investi-gate in detail the role of the network structure and ofdifferent temporal scales in neurocognitive behavior.
ACKNOWLEDGMENTS
LdA would like to thank MIUR projectPRIN2017WZFTZP for financial support. AS acknowl-edges support from MIUR project PRIN201798CZLJ.LdA and AS acknowledge support from Program (VAn-viteLli pEr la RicErca: VALERE) 2019 financed by theUniversity of Campania “L. Vanvitelli”.
Appendix A: Experimental procedure, Dataacquisition and Preprocessing
Spontaneous resting state activity and stimulus-evokedresponse were recorded from healthy human subjects (n= 7, age= 22 ± C ( t ) / C ( ) C ( t ) / C ( ) time (ms) C ( t ) / C ( ) time (ms) FIG. 3. Autocorrelation function of the global spontaneousactivity for subjects 1-6 (a-f). -0,2 0 0,2 0,4 0,6 0,8 1 s i gn a l ( T ) -0,2 0 0,2 0,4 0,6 0,8 1 -0,2 0 0,2 0,4 0,6 0,8 1 s i gn a l ( T ) -0,2 0 0,2 0,4 0,6 0,8 1 -0,2 0 0,2 0,4 0,6 0,8 1 time (s) s i gn a l ( T ) -0,2 0 0,2 0,4 0,6 0,8 1 time (s) x10 -11 x10 -11 (a) (b)(c)(e) (d)(f) FIG. 4. Signal of evoked activity for subjects 1-6 (a-f). Thevisual stimulation is applied at time 0. artefacts, as determined by visual inspection of the 2Dscalp maps and time course of that ICA component, wererejected and remaining components were used to recon-struct the data. For additional information, see Ref.[26].The resultant absolute amplitude were summed acrossthe MEG sensors and formed a global signal of the asso-ciated brain activity.
Appendix B: Autocorrelation and response functionsfor different functional areas and different subjects
In Fig. 3 and Fig. 4 we report the (global) time auto-correlations and response functions, respectively, for theseveral observed subjects.We also show further data analyses considering, ratherthan the global signal from all sensors, 1) the signal fromsensors corresponding to visual areas (in the occipitallobe) and 2) the signal from sensors corresponding tothe temporal lobes, which include the area specialized inface recognition, e.g., the fusiform face area. The later
Subject τ global occ. temp. τ global occ. temp. A global occ. temp. τ R global occ. temp.1a 7.78 7.99 8.09 259 267 249 0.902 0.906 0.930 83.3 45.5 96.441b 7.84 7.87 7.93 367 394 366 0.866 0.864 0.9052a 10.82 9.45 10.1 515 662 598 0.789 0.756 0.84 52.6 58.7 42.22b 8.51 8.70 8.71 400 425 402 0.775 0.722 0.8093a 7.75 7.60 8.22 537 568 740 0.757 0.849 0.890 37.0 35.7 43.53b 8.17 7.28 8.66 698 498 881 0.725 0.754 0.7724a 7.65 6.42 7.20 437 661 591 0.832 0.874 0.921 90.9 80.4 1084b 7.37 6.37 6.68 582 689 544 0.752 0.747 0.7945a 10.65 8.80 10.6 442 455 439 0.742 0.897 0.878 22.7 27.2 23.55b 9.79 7.34 8.43 637 633 674 0.877 0.937 0.9346a 8.29 8.28 8.08 1049 722 686 0.793 0.805 0.741 52.6 37.8 54.16b 7.04 7.52 7.04 258 409 244 0.610 0.751 0.5867a 10.84 9.32 9.34 691 887 767 0.894 0.922 0.921 41.7 128.7 30.27b 10.77 8.93 9.76 815 890 1314 0.855 0.872 0.888Average 8.8 ± ± ± ±
200 583 ±
190 607 ±
280 0.8 ± ± ± ±
25 60 ±
35 57 ± f ( t ) = Ae − t/τ +(1 − A ) e − t/τ to experimentaldata of rest activity (two data set for each subject) and from a single exponential fit with the function g ( t ) = e − t/τ R for evokedactivity, for all the subjects, for global activity, occipital area and temporal area. Times are measured in ms. are also the regions that demonstrated the strongest ac-tivity. For both data sets, we analysed the temporal au-tocorrelations for rest activity and the response in thetotal signal from each area and for individual sensors ineach area. We also present data analysis for the differentsubjects.In Tab. I we report the values of the parameters ob-tained from the fit of the theoretical predictions of the lin-earized Wilson-Cowan to the experimental data, for restand evoked activity, for global activity and for occipitaland temporal area. For each subject we have two datasetsof spontaneous activity and one data set for evoked ac-tivity. Note that the values of τ are quite stable withinthe set of subjects, while the values of τ show a largervariability.In figures 5 and 6, we show the results for the autocor-relation and the response function compared with resultsfor global data, for one subject (similar behavior was ver-ified also for the other six subjects). In Figs. 5 thin linesshow data for each single sensor, whereas the thick blackline is the average over the eight sensors. Data show thatthe autocorrelation function (in linear and log-linear scalein the inset) calculated at a single sensor both in the oc-cipital and temporal area and averaged over eight sensorsof the same area, confirm the double exponential decaywith characteristic times similar to the ones for globalactivity. Moreover, in Fig. 6 (left panel), the autocorre-lation function calculated for the global signal for eacharea displays the same functional behavior with the samecharacteristic times found for the global signal.Fig. 6 (bottom panel) shows the response function forthe average signal from the temporal and occipital areas(thick lines, averaged over the number of trials), as wellas for signals at two representative sensors in each area. In this case, data for each sensor are, as expected, verynoisy in all cases, however the fit of the thick lines isagain optimal with a single exponential decay where thecharacteristic time (43ms for the temporal area) is ingood agreement with τ R found for the global signal. Appendix C: Details on analytical computations
The Wilson-Cowan model in the linear noise approxi-mation consists of the linear system˙ ξ = M ξ + p α Σ η , (C1)where ξ = ( ξ Σ , ξ ∆ ) T , M = − /τ w ff − /τ ! , (C2)and η = ( η Σ , η ∆ ) T are independent delta-correlatednoises, with zero average and unit variance. The gen-eral solution of Eq. (C1) is ξ ( t ) = e M t ξ (0) + p α Σ Z t e M ( t − t ′ ) η ( t ′ ) dt ′ . (C3)The covariance matrix at stationarity σ , whose matrix el-ements are the equal-time correlations σ ij = h ξ i (0) ξ j (0) i with i, j = (Σ , ∆), where h· · · i denotes an average in thestationary state, satisfies the matrix relation [30] − α Σ α Σ ! = M σ + σ M T , (C4) time (ms) C ( t ) / C ( ) averagedouble exponential fit (A=0.97, τ =6.9 ms, τ =376 ms)0 100 200 300 400 5000.010.11 (a) time (ms) C ( t ) / C ( ) averagedouble exponential fit (A=0.97, τ =7.1 ms, τ =309 ms)0 100 200 300 400 5000.010.11 (b) FIG. 5. Autocorrelation function for the spontaneous activ-ity at eight sensors in the occipital (a) and temporal (b) areasof a single subject: Thin lines show data for each single sen-sor, whereas the thick black line is the average over the eightsensors. where M T denotes the transpose matrix. SolvingEq. (C4) one obtains σ = α Σ τ (cid:16) w ff τ τ τ + τ (cid:17) w ff τ τ τ + τ w ff τ τ τ + τ τ . (C5)The inverse matrix σ − then reads σ − = 2 α Σ τ + τ τ τ + τ + τ (1 + w ff τ ) × τ + τ τ − w ff τ − w ff τ τ + τ ++ w ff τ τ τ ! . (C6)The elements of the time correlation matrix C ( t ) are ob-tained from the equations [30] C ij ( t ) = ( e M t σ ) ij . (C7)The matrix M has eigenvalues ( − /τ , − /τ ) and eigen-vectors (1 , T and ( − w ff τ τ / ( τ − τ ) , T . Diagonal- time (ms) C ( t ) / C ( ) signal from all sensorssignal from sensors in the occipital areasignal from sensors in the temporal areadouble exponential fit (A=0.9, τ =7.8, τ =257)0 100 200 300 400 5000.010.11 (a) -0.2 0 0.2 0.4 0.6 0.8 1 time (s) -0.500.51 R ( t ) occipital area (average over 8 sensors)temporal area (average over 8 sensors)occipital area (sensor A105, averaged over trials)temporal area (sensor A210, averaged over trials)single exponential fit ( τ R =43 ms) (b) FIG. 6. (a) The autocorrelation function calculated for theglobal signal for each area. (b) The response function for theaverage signal from the temporal and occipital areas (thicklines, averaged over the number of trials), as well as for signalsat two representative sensors in each area. izing M , one obtains the matrix exponential e M t = − w ff τ τ τ − τ ! e − t/τ e − t/τ ! × − w ff τ τ τ − τ ! − = e − t/τ τ τ w ff ( e − t/τ − e − t/τ ) τ − τ e − t/τ ! . (C8)We notice that in the case of equal characteristic times τ = τ the upper right element in the above matrix takesthe value te − t/τ . Next, computing the matrix productin Eq. (C7) one obtains the explicit expressions C ΣΣ ( t ) = α Σ τ τ τ − τ ) × h ( τ − − τ τ − − τ w ff ) e − t/τ + τ w ff e − t/τ i (C9) C Σ∆ ( t ) = α Σ τ τ w ff τ − τ ) h τ e − t/τ − ( τ + τ ) e − t/τ i (C10) C ∆Σ ( t ) = α Σ τ τ w ff τ + τ ) e − t/τ (C11) C ∆∆ ( t ) = α Σ τ e − t/τ . (C12) We notice that C ∆Σ ( t ) is different than zero and describesthe decay of the initial correlation ξ ∆ (0) ξ Σ (0), controlledby the characteristic time τ .The matrix exponential e M t coincides with the re-sponse function matrix R ij = δξ i ( t ) /δξ j (0), so that R ΣΣ ( t ) = e − t/τ (C13) R Σ∆ ( t ) = τ τ w ff ( e − t/τ − e − t/τ ) τ − τ (C14) R ∆Σ ( t ) = 0 (C15) R ∆∆ ( t ) = e − t/τ . (C16) [1] A. Arieli, A. Sterkin, A. Grinvald, A. Aertsen, Science , 1868 (1996).[2] M. D. Fox, A. Z. Snyder, J. M. Zacks and M. E. Raichle,Nature Neuroscience , 23 (2006).[3] M. D. Fox and M. E. Raichle, Nature Reviews Neuro-science , 700 (2007).[4] B. J. He, J. Neurosci. , 4672 (2013).[5] Z. Huang, J. Zhang, A. Longtin, G. Dumont, N. W. Dun-can, J. Pokorny, P. Qin, R. Dai, F. Ferri, X. Weng, andG. Northoff, Cerebral Cortex , 1037 (2017).[6] K. Linkenkaer-Hansen, V. V. Nikouline, J. M. Palva, andR. J. Ilmoniemi, J. Neurosci. , 1370 (2001).[7] J. D. Murray, A. Bernacchia, D. J. Freedman, R. Romo,J. D. Wallis, X. Cai, C. Padoa-Schioppa, T. Pasternak,H. Seo, D. Lee, and X.-J. Wang, Nature Neuroscience ,1661 (2014).[8] C. J. Honey, T. Thesen, T. H. Donner, L. J. Silbert, C.E. Carlson, O. Devinsky, W. K. Doyle, N. Rubin, D. J.Heeger, U. Hasson, Neuron , 668 (2012).[9] G. J. Stephens, C. J. Honey, U. Hasson, J Neurophysiol , 2019 (2013).[10] C. Meisel, K. Bailey, P. Achermann, and D. Plenz, Sci.Rep. , 11825 (2017).[11] C. Meisel, A. Klaus, V. V. Vyazovskiy, and D. Plenz, J.Neurosci. , 10114 (2017).[12] I. Ferezou, T. Deneux, Neurophoton. , 031221 (2017).[13] M.W. Cole, T. Ito, D. Schultz, R. Mill, R. Chen, C.Cocuzza, Neuroimage , 111 (2007).[15] L. F. Cugliandolo, J. Phys. A: Math. Theor. , 483001(2011).[16] U. Seifert, Rep. Prog. Phys. , 126001 (2012).[17] M. Baiesi and C. Maes, New J. Phys. , 013004 (2013).[18] A. Puglisi, A. Sarracino, and A. Vulpiani, Phys. Rep. , 1 (2017).[19] A, Sarracino and A. Vulpiani, Chaos , 083132 (2019).[20] H. R. Wilson and J. D. Cowan, Biophys. J. , 1 (1972).[21] T. Ohira, J.D. Cowan, in Mathematics of neural net-works: models, algorithms, and applications, Ellacort S,Anderson I, eds. Springer, p. 290 (1997).[22] M. Benayoun, J. D. Cowan, W. van Drongelen, E. Wal-lace, PLoS Comput Biol. , 1 (2010).[23] E. Negahbani, D.A. Steyn-Ross, M.N. Steyn-Ross, M.T.Wilson, J.W. Sleigh, J. Math. Neurosci. 5, 1 (2015)[24] C. Zankoc, T. Biancalani, D. Fanelli, and R. Livi, Chaos,Solitons and Fractals , 504 (2017).[25] S. di Santo, P. Villegas, R. Burioni, and M. A. Mu˜noz,J. Stat. Mech. , 073402 (2018).[26] O. Arviv, A. Goldstein, and O. Shriki, J. Neurosci. ,13927 (2015).[27] O. Arviv, A. Goldstein, O. Shriki, Sci. Rep. 9, 13319(2019).[28] F. Lombardi, H.J. Herrmann, C. Perrone-Capano, D.Plenz, L. de Arcangelis, Phys. Rev. Lett. , 228703(2012)[29] B. K. Murphy and K. D. Miller, Neuron , 635 (2009).[30] Risken, H. (1996). The Fokker-Planck Equation.Springer, Berlin, Heidelberg.[31] O. Shriki, D. Hansel, H. Sompolinski, Neural Computa-tion , 1809 (2003).[32] A. Sarracino, D. Villamaina, G. Gradenigo, and A.Puglisi, EPL , 34001 (2010).[33] A. Crisanti, A. Puglisi, and D. Villamaina, Phys. Rev. E , 1 (2012).[34] I. Tal and M. Abeles, J. Neurosci. Methods , 31(2013).[35] R. Oostenveld, P Fries, E. Maris, and J.-M. Schoffelen,Comput. Intell. Neurosci. , 156869 (2011).[36] A. Jain, R. Bansal, A. Kumar, K.D. Singh, Int. J. Appl.Basic Med. Res.5