Predicting brain evoked response to external stimuli from temporal correlations of spontaneous activity

The relation between spontaneous and stimulated global brain activity is a fundamental problem in the understanding of brain functions. This question is investigated both theoretically and experimentally within the context of nonequilibrium fluctuation-dissipation relations. We consider the stochastic coarse-grained Wilson-Cowan model in the linear noise approximation and compare analytical results to experimental data from magnetoencephalography (MEG) of human brain. The short time behavior of the autocorrelation function for spontaneous activity is characterized by a double-exponential decay, with two characteristic times, differing by two orders of magnitude. Conversely, the response function exhibits a single exponential decay in agreement with experimental data for evoked activity under visual stimulation. Results suggest that the brain response to weak external stimuli can be predicted from the observation of spontaneous activity and pave the way to controlled experiments on the brain response under different external perturbations.


I. INTRODUCTION
The brain represents one of the most fascinating systems where several mechanisms at different scales are deeply intertwined, resulting in a complex behavior.One of the main open issues in the understanding of brain functioning is the relation between spontaneous and evoked activity, namely the response of the system to external stimuli.In particular, it has been found that the large variability in the response to repeated presentations of the same stimulus can be attributed to the ongoing spontaneous activity [1][2][3].However, a theoretical framework to formalize this question is still lacking and the quantitative connection between the spontaneous and evoked activity remains unclear [4,5].
Experimental results for temporal correlations of spontaneous brain activity have been reported in a number of studies.For instance, the seminal article [6] focused on correlations for spontaneous alpha oscillations in the healthy human brain, and found a decay characterized by power-law behavior at long times.More recent analyses [7] measured autocorrelations of spiking activity fluctuations in several cortical areas of the macaque monkey.For each area an intrinsic timescale is defined from the exponential fit of the autocorrelation, which decays to a nonzero offset value taking into account longer timescales.A hierarchical ordering of intrinsic characteristic times across areas was evidenced in [8].In particular, a hierarchy of time scales coupling the dynamics of sensory regions at different ordering levels, was found in the auditory functional area [9].The observed relation between time-scales suggests a temporal organization across the cerebral cortex.Other studies on the cortex dynamics, focusing on the effects of sleep deprivation on time correlations, observed exponential decays: Sustained wakefulness appears to decrease correlation timescales in humans [10], and cortical timescales depend on time awake and sleep stage in rats [11].Even if spontaneous activity has been found to be highly structured and to participate in high cognitive functions, its functional role remains poorly understood.In particular, it is well accepted that variability is reduced during tasks at different scales, yet there is no general consensus whether task activation would inflate or reduce functional connectivity and correlations [12,13].The question we address here is to provide a theoretical formulation for the relation between the response and the correlation function of rest activity by means of a simple neuronal model suitable for an analytical approach.
The relation between ongoing and stimulated activity can be addressed theoretically within the general framework of statistical physics, by means of the fluctuationdissipation relations, connecting the spontaneous fluctuations of a system with the response function to external perturbations.Recently, these relations have been extended beyond the context of standard thermodynamics, to nonequilibrium systems, such as active and biological matter [14][15][16][17][18][19].These relations hold for a wide variety of physical systems, not necessarily at the critical point.Their deep predictive power relies on the possibility to estimate the response of a system from the observation of its spontaneous fluctuations.However, in order to write explicit relations, a theoretical model describing the system dynamics at the scale of interest is needed.To this extent, a useful approach to describe the activity of a large number of neurons is provided by the coarse-grained Wilson-Cowan model (WCM) [20,21].The stochastic version of this model can describe the fluctuations of active excitatory and inhibitory neuron populations via two Langevin equations, whose dynamics is coupled through a feed-forward term, and can reproduce the statistics of burst activity [22].Recent studies focused on dynamical stability properties of this model revealing a rich dynamical behavior, see for instance [23][24][25].
In this work, we address the fundamental problem of unveiling the quantitative relation between spontaneous and stimulated brain activity.We compare the analytical expression for correlations and response function in the WCM with spontaneous and stimulated brain activity measured via magnetoencephalography (MEG).We measure the time correlation functions of spontaneous activity in several healthy subjects, finding a temporal behavior mainly characterized by a double-exponential decay.This two-timescale decay observed in experiments is in good agreement with the analytical prediction of the WCM.We also calculate, in the linear response regime, the response function to an external stimulation and compare the result to MEG data of evoked activity, obtained by applying a visual stimulation (pictures of faces) to the same subjects.Our study enlightens how some properties of the induced brain dynamics can be predicted from the observation of spontaneous activity in the absence of stimuli.
The paper is organized as follows.In Sec.II we briefly recall the definition of the WCM, focusing on its linearized version, which is expected to hold in the large size limit.Then, in Sec.III, we report the experimental results for the spontaneous activity autocorrelation functions from MEG data.In Sec.IV we derive the fluctuation-dissipation relation for the WCM and compare its predictions with experimental data for evoked activity.Finally some conclusions are drawn in the last section.Appendices A, B, and C report details on the experimental procedure and on analytical calculations.

II. LINEARIZED WILSON-COWAN MODEL
In order to provide a theoretical framework to understand the relation between spontaneous and evoked brain activity, we consider the linearized version of the stochastic coarse-grained WCM, which allows us to describe the systems dynamics at meso-and macroscales [21].This model describes the excitatory and inhibitory neuron population dynamics, which evolve according to a Master Equation for the probability p k,l (t) to have k excitatory and l inhibitory neurons active at time t [22].The total number of neurons is N E + N I , where N E and N I are the populations of excitatory and inhibitory neurons, respectively.Assuming N E = N I = N , for large N the number of active neurons is written as the sum of a deterministic component (E, I) = (k/N, l/N ), and a stochastic fluctuation (ξ E , ξ I ), scaled by √ N , i.e. k = N E + √ N ξ E and l = N I + √ N ξ I .The Master Equation can be rewritten in terms of these new variables and expanded in power of N obtaining a Fokker-Planck equation.This leads to the dynamic equations for the deterministic and stochastic components.Next, introducing the global variables [22] in the large N limit and close to the fixed point (Σ 0 , ∆ 0 = 0) the deterministic components satisfy the equations where α is the decay rate of the active state, f (s) = tanh(s) for s > 0 and zero otherwise, τ 0 = 1 ms is a fixed microscopic time-scale and s = w 0 Σ + (w E + w I )∆ + h.
Here w E (w I ) is the synaptic strength of the excitatory (inhibitory) population, w 0 = w E − w I and h a small external input.We notice that Σ represents the global neuronal activity, whereas ∆ measures the unbalance in activity between the excitatory and the inhibitory population.The deterministic equations have the unique stable solution (Σ 0 , 0).Note that the fixed point ∆ 0 = 0 represents the condition of balance between excitation and inhibition, as found for spontaneous activity of neuronal systems in healthy state [28].
Concerning the fluctuations of Σ and ∆, ξ Σ and ξ ∆ , in the large N limit (the linear noise approximation) and close to the fixed point (Σ 0 , ∆ 0 = 0), they satisfy [22] where η Σ and η ∆ are zero average, delta-correlated white noises with unitary variance, and 1/τ τ 0 , with s 0 = w 0 Σ 0 + h.The offdiagonal term w f f is called hidden feed-forward term and plays a central role because it rules the coupling between the global activity and the unbalance between the activities of excitatory and inhibitory populations [29].In particular, a fluctuation in ∆ affects the temporal evolution of Σ but not vice-versa.
Eqs. (3) are two coupled linear Langevin equations, the dynamics of which can be easily solved analytically [30].In particular, the solution in vectorial form reads where M is the coupling matrix appearing in Eq. (3).Therefore, we evaluate the correlation matrix at stationarity where ξ Σ (t)ξ Σ (0) , takes the form Therefore, the behavior of the correlations of the active neuron population is characterized by a doubleexponential decay, with two characteristic times τ 1 and τ 2 , which are a function of the model parameters as specified above.Similar double exponential decay is also found for the C Σ∆ (t), whereas a single exponential function with characteristic time τ 2 is obtained for C ∆Σ (t) and C ∆∆ (t) (see Appendix C).

III. TIME CORRELATION OF SPONTANEOUS ACTIVITY FROM MEG DATA
The variable Σ describes the brain spontaneous activity and therefore allows for a comparison with experimental data obtained from MEG measurements.We analyse the global signal, by summing the signals at all sensors, the signal from sensors monitoring the visual area in the occipital lobe and the signal from sensors monitoring the temporal lobes, which include the area specialized in face recognition, e.g., the Fusiform Face Area (FFA).Spontaneous activity was recorded from 7 healthy human subjects with a whole-head 248-channel magnetometer array (4-D Neuroimaging, Magnes 3600WH) in the MEG facility at the EMBI Unit, Bar-Ilan University, Israel.The MEG signal was recorded at a sampling rate of 1017.25 Hz and offline band-pass filtered between 0.8 and 80 Hz as well as underwent a cleaning procedure of potential artefacts.The absolute amplitudes were summed across the 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 from time series (of total length about 4 minutes for each subject) of the absolute value {|x(t)|} of rest activity signal, with mean µ and standard deviation σ 0 .
We focus here on the short time behavior (up to one second), where the experimental data show a decay mainly characterized by two typical times (see Fig. 1).In this study we neglect oscillations at very short time scales, which probably correspond to α brain rhythm, and consider the global decay behavior.Analyzing 7 subjects (2 data sets for each subject), we observe that this complex behavior is quite stable (see Table in Appendix B) and can be fitted by Eq. ( 6).Moreover, the same functional behavior is found for the global signal and for the signal averaged only on sensors from temporal and occipital lobes, corresponding to areas involved in processing of the stimuli.The double exponential behavior is also detected for data at single sensors in these two areas, where the autocorrelation function, as expected, exhibits larger statistical fluctuations (see Appendix B).The average values of the characteristic times are found to be τ 1 = 8.8 ± 1.5 ms and τ 2 = 515 ± 200 ms, differing by almost two orders of magnitude.The short time scale, causing an abrupt decay in the correlation function, is of the order of magnitude of synaptic time scales as well as single neuron time scales, such as the refractory period [31].Conversely, the long characteristic time is compatible with the time scale of low-frequency rhythms, as the θ rhythm.Interestingly, the long characteristic time τ 2 controls the single exponential decay of the correlations C ∆Σ (t) and C ∆∆ (t) (see Appendix C).This suggests that a fluctuation in Σ or ∆ affects the temporal behavior of ∆ itself over a long timescale, implying a slow recovery of the balance condition if the system is moved away from the fixed point ∆ 0 = 0. From the fitting procedure we also obtain an estimation for the feed-forward coefficient w f f = 0.008 ± 0.0035 ms −1 , which characterizes the coupling between excitatory and inhibitory neuron populations.Finally, we note that, at longer times (t 2 s), slower decay (power-law) of the autocorrelation function can take place.This functional decay cannot be obtained by a linearized model.

IV. PREDICTING RESPONSE FROM SPONTANEOUS ACTIVITY VIA THE FLUCTUATION-DISSIPATION THEOREM
The present theoretical framework allows us to address the fundamental question of the relation between spontaneous and evoked activity.We can evaluate the linear response function to a weak external perturbation, de-fined as where i, j = (Σ, ∆).This represents the average response of the variable ξ i (t) at time t to an impulsive perturbation applied to the variable ξ j at time 0, and • denotes nonstationary averages over many realizations.From Eq.( 4) one immediately obtains the response matrix From Eq. ( 5), we obtain a fluctuation-dissipation relation connecting the linear response with spontaneous fluctuations [14].In particular, since we are interested in the response of the total activity Σ to a small perturbation of Σ itself, the response function reads where the explicit form of the inverse matrix σ −1 is reported in Appendix C. The relevance of relation (11) relies on the fact that it allows us to get information on the response of the system to an external perturbation by simply looking at its unperturbed dynamics.In general, as evident from Eq. ( 11), both the autocorrelation and the cross-correlation are required to reconstruct the response to a weak external stimulus.This is due to the presence of nonequilibrium conditions breaking detailed balance, as noted in several systems (see for instance [14,32,33]).
Let us note however, that, in our particular case, due to the upper triangular form of the matrix M in Eq. ( 3), the response function of the model takes the simple exponential form where the characteristic time τ 1 is the same ruling the short-time behavior of the correlation function.This is consistent with Eq. ( 11) because one can easily check that the cross-correlation appearing in Eq. ( 11) exactly cancels the term ∝ exp(−t/τ 2 ) in Eq. ( 6), giving the single-exponential decay for the response function (see Appendix C).Therefore, in this approach, measurements of the spontaneous fluctuations in the global brain activity Σ alone could provide a prediction for the system response.Concerning the other response functions, we notice that R Σ∆ (t) shows a double exponential behavior, whereas, as expected from the triangular form of M, R ∆Σ (t) = 0 and R ∆∆ (t) = exp(−t/τ 2 ) (see Appendix C).The response function derived by the analytical calculation characterizes the behavior close-in-time to the stimulus application.In order to fairly compare the analytical prediction to experimental data, the experimental perturbation has to be sufficiently weak to remain in the linear regime.Face processing is one of the most studied cognitive abilities.Humans are considered experts in face processing, which accordingly does not involve long-lasting cognitive engagement, potentially mimicking a small perturbation to the system.The evoked response is fast and of relatively short (typically characterized by 100, 170 and 250 ms evoked response fields (ERFs), and some late components at 400 and 600 ms) and the return to baseline in both power and spectral content occurs prior to the 1 sec limit [26,27].We measured the brain activity of the same 7 participants evoked by visual stimuli, by showing a series of pictures of human faces.Each stimulus is presented for 1 sec, with varying inter stimulus intervals (larger than 1 sec).Further details are reported in Appendix A.
In Fig. 2 we show the response function averaged over the 7 subjects for the global signal and for the signal averaged only on sensors from temporal and occipital lobes, corresponding to areas involved in processing of the stimuli (see Appendix B for the response function for each subject and at single sensor).Since each subject shows a different latency activity associated to each stimulus which cannot be accounted for by analytical predictions, experimental data for the response have been shifted in time and rescaled in amplitude in order to have that the maximum in each dataset is reached at time zero and its value is normalized to 1.We have to stress that the applied perturbation cannot be considered impulsive, as in the analytical definition, and moreover the system takes a certain time to develop the response signal, as shown in Appendix B. However, experimental data on evoked activity qualitatively confirm the single exponential decay predicted by the WCM for all datasets.Indeed, imposing a fit with a double exponential leads in all cases to an amplitude close to zero for the second exponential decay.The exponential behavior is also detected for data at single sensors, where the response function, as expected, exhibits larger statistical fluctuations (see Appendix B).The averaged fitted value for the characteristic time τ R = 55 ± 25ms, which is of order of magnitude comparable to τ 1 (see Table in Appendix B for the characteristic time values in each dataset).The quantitative difference between these two timescales could be due to a number of factors, as the properties of the stimulation, which in the theoretical approach is supposed to be instantaneous and small.The real stimulation requires a time delay before the decay sets in, raising the need to identify the initial time to fit the response function.Fluctuations in experimental conditions then provide a wider range of characteristic times that, averaged over trials, lead to a larger τ R .Random fluctuations also hide a possible correlation between τ 1 and τ R across subjects.Response experiments under different stimulation protocols and controlled conditions are required to test more quantitatively the timescale agreement.The robust functional behavior found for global signals, functional areas signals and data at single sensor, is an indication of scale-invariance, namely the correlation and the response function do not depend on the sample size, therefore evidencing the self-similar properties of the system.This observation could be the outcome of the coarse-grained measure of the activity by MEG, since sensors are spaced at ∼ 2cm and therefore the signal at a single sensor is already representative of the activity of a large population of neurons.We are aware that at a finer scale the response to a stimulus has a complex spatio-temporal organization, however this information cannot be accounted for by our analytical approach.
Let us finally observe that within this approach one can derive another relevant quantity characterizing the nonequilibrium behavior of the neuronal system, namely the entropy production rate, which characterizes how much the system is far from equilibrium [33] This makes clear that the nonequilibrium source in this model is the presence of the feed-forward term w f f , which couples the fluctuations in the two variables Σ and ∆.

V. CONCLUSIONS
We have proposed a theoretical approach to predict the relation between global spontaneous and evoked brain activity, based on the linear noise approximation of the WCM, which allows for analytical calculations.The theoretical approach relies on two main assumptions: The model is two-dimensional and parameters are tuned to achieve balanced excitation and inhibition, which leads to the fixed point ∆ = 0.The comparison with MEG data confirms that temporal correlations of spontaneous fluctuations are mainly characterized by a double-exponential decay, with two well separated typical times τ 1 ≃ 10 ms and τ 2 ≃ 500 ms.The short time scale is related to synaptic timescales, while the long one is compatible with slow brain rhythms.This analytical approach is therefore able to rationalize experimental results for the correlation and response functions providing a coherent framework for the dual functional behavior found experimentally.An important aspect of our results is that the functional behavior found in MEG data is scale-invariant, with stable value of the characteristic times across scales, from the cm scale (single sensor data) to the entire brain (global signal).This observation suggests that at first approximation the linear model well accounts for the brain behavior at different scales, even if the connectome has a complex modular structure.
The presence of a single exponential decay with the short characteristic time for the response function is fully coherent with the efficient performance of the human brain in visual tasks.Indeed, it is well known that the human eye can appreciate about 10 images per second, resulting in a temporal resolution of about 100 ms [36].A larger relaxation time, of the order of τ 2 ≃ 500 ms, would imply an overlap in the response to close-in-time stimulations, affecting the performance.We deem that future experiments on evoked activity, properly designed to the application of the fluctuation-dissipation theorem, could shed new light on the brain functionality at large scale, with strong impact in neurobiology and neuroscience.Moreover, this result opens the way to numerical studies implementing microscopic models, as integrate and fire neuronal models, on complex networks able to investigate in detail the role of the network structure and of different temporal scales in neurocognitive behavior.
head movements throughout the recordings, head localization measurements were performed before and after each experiment.Head position and shape were determined by Pollhemus FASTTRAK digitizer and five coils attached to the participant's head, measuring position relative to the MEG sensors.The MEG was recorded at a sampling rate of 1017.25 Hz and analog band-pass filtered online at 0.1-400Hz.Reference coils were used to remove environmental noise.Accelerometers (Bruel and Kjaer) attached to the gantry were used to remove vibration noise.The 50-Hz signal from the power outlet was recorded by an additional channel and the average power-line response to a power cycle was subtracted from every MEG sensor [34].
E-prime 2.0 (Psychology Software Tools Inc.) was used for experimental control.During rest, participants were instructed to fixate their eyes on a fixation cross at the center of a black screen.During evoked, the stimuli were grey-scale pictures of human faces.Each stimulus was presented for 1000 ms with inter-stimulus intervals varying between 1,300 and 1,700 ms.Stimuli were back projected on a screen placed in front of the subjects, by a video projector situated outside the room.To each subject, 540 pictures of faces were presented (details: photographs of 5 different human male models in 9 head posters and emotional expressions, giving 45 pictures and photographs of 3 different human female models with 3 emotional expressions, giving 9 pictures; All pictures were repeated randomly 10 times, once in each experimental block, that is 10 blocks).During the MEG scan, participants completed an oddball gender-detection task, pressing a response button only when a female face was presented (16.67%).This ensured that all 450 face presentation of male models were task irrelevant.Only these presentations underwent analysis.Subsequently, among these 450 picture presentations, between 3 and 30 trials per subject were removed during the cleaning procedure due to artefacts.
Data processing was performed using Fieldtrip opensource toolbox for Advanced MEG Analysis [35].MEG recordings were first cleaned for line frequency, building vibration and heartbeats artefacts with an in-house opensource software based on external cues [34].Stimulusevoked data were segmented to include the 1 sec trials as well as an additional 0.2 pre-trial interval and head and tail of 0.4 sec that were later cut from analysis.All data were band-pass filtered offline between 0.8 and 80 Hz (zero-phase two-pass Butterworth IIR filter of order 4 and 53 dB stopband attenuation, upper band limit was chosen to minimize the effect of muscle artifacts).Epochs containing a false-positive response or contaminated by jump in the MEG sensors or muscle artefacts, displaying variance higher than 3 SD in power above 60 Hz, were discarded.One mal-functioning MEG sensor was discarded from all datasets.Independent component analysis (ICA) was performed on the remaining data, to ensure the removal of eye-movements, blinks and leftover heartbeats artefacts.ICA components reflecting such -0,2 0 0,2 0,4 0,6 0,8 1 artefacts, as determined by visual inspection of the 2D scalp maps and time course of that ICA component, were rejected and remaining components were used to reconstruct the data.For additional information, see Ref. [26].
The resultant absolute amplitude were summed across the MEG sensors and formed a global signal of the associated brain activity.are also the regions that demonstrated the strongest activity.For both data sets, we analysed the temporal autocorrelations for rest activity and the response in the total signal from each area and for individual sensors in each area.We also present data analysis for the different subjects.
In Tab.I we report the values of the parameters obtained from the fit of the theoretical predictions of the linearized Wilson-Cowan to the experimental data, for rest and evoked activity, for global activity and for occipital and temporal area.For each subject we have two datasets of spontaneous activity and one data set for evoked activity.Note that the values of τ 1 are quite stable within the set of subjects, while the values of τ 2 show a larger variability.
In figures 5 and 6, we show the results for the autocorrelation and the response function compared with results for global data, for one subject (similar behavior was verified also for the other six subjects).In Figs. 5 thin lines show data for each single sensor, whereas the thick black line is the average over the eight sensors.Data show that the autocorrelation function (in linear and log-linear scale in the inset) calculated at a single sensor both in the occipital and temporal area and averaged over eight sensors of the same area, confirm the double exponential decay with characteristic times similar to the ones for global activity.Moreover, in Fig. 6 (left panel), the autocorrelation function calculated for the global signal for each area displays the same functional behavior with the same characteristic times found for the global signal.
Fig. 6 (bottom panel) shows the response function for the average signal from the temporal and occipital areas (thick lines, averaged over the number of trials), as well as for signals at two representative sensors in each area.
In this case, data for each sensor are, as expected, very noisy in all cases, however the fit of the thick lines is again optimal with a single exponential decay where the characteristic time (43ms for the temporal area) is in good agreement with τ R found for the global signal.where M T denotes the transpose matrix.Solving Eq. (C4) one obtains (C5) The inverse matrix σ −1 then reads The elements of the time correlation matrix C(t) are obtained from the equations [30] C ij (t) = (e Mt σ) ij . (C7) The matrix M has eigenvalues (−1/τ 1 , −1/τ 2 ) and eigenvectors (1, 0) T and (−w f f τ 1 τ 2 /(τ 1 − τ 2 ), 1) T .Diagonal- We notice that in the case of equal characteristic times τ 1 = τ 2 the upper right element in the above matrix takes the value te −t/τ1 .Next, computing the matrix product in Eq. (C7) one obtains the explicit expressions We notice that C ∆Σ (t) is different than zero and describes the decay of the initial correlation ξ ∆ (0)ξ Σ (0), controlled by the characteristic time τ 2 .

1 FIG. 2 .
FIG.2.(Rescaled and shifted) average response function for the global signal, the signal from the occipital and the temporal lobes.The typical time for the response function from an exponential fit decay, according to Eq. (12), is τR = 55 ± 25 ms.In the inset the same data are plotted in semi-log scale.

FIG. 5 .
FIG.5.Autocorrelation function for the spontaneous activity at eight sensors in the occipital (a) and temporal (b) areas of a single subject: Thin lines show data for each single sensor, whereas the thick black line is the average over the eight sensors.

1 =
FIG. 6.(a) The autocorrelation function calculated for the global signal for each area.(b) The response function for the average signal from the temporal and occipital areas (thick lines, averaged over the number of trials), as well as for signals at two representative sensors in each area.
• • • denote average over noise, i, j = (Σ, ∆) and σ is the covariance matrix which is a known function of the model parameters (see Appendix C).The specific case of the autocorrelation for the total activity Σ, C ΣΣ ≡

TABLE I .
Parameters obtained from the double-exponential fit with the function f (t) = Ae −t/τ 1 +(1−A)e −t/τ 2 to experimental data of rest activity (two data set for each subject) and from a single exponential fit with the function g(t) = e −t/τ R for evoked activity, for all the subjects, for global activity, occipital area and temporal area.Times are measured in ms.