PHYSICAL REVIEW E 91, 012702 (2015)

Existence of anticorrelations for local field potentials recorded from mice reared in standard condition and environmental enrichment F. Vallone,1,2,3 A. Cintio,1 M. Mainardi,3,* M. Caleo,3 and A. Di Garbo1,† 1

Institute of Biophysics, CNR—National Research Council, 56124 Pisa, Italy 2 The BioRobotics Institute, Scuola Superiore S. Anna, 56026 Pisa, Italy 3 Neuroscience Institute, CNR—National Research Council, 56124 Pisa, Italy (Received 7 February 2014; revised manuscript received 28 July 2014; published 9 January 2015) In the present paper, we analyze local field potentials (LFPs) recorded from the secondary motor cortex (M2) and primary visual cortex (V1) of freely moving mice reared in environmental enrichment (EE) and standard condition (SC). We focus on the scaling properties of the signals by using an integrated approach combining three different techniques: the Higuchi method, detrended fluctuation analysis, and power spectrum. Each technique provides direct or indirect estimations of the Hurst exponent H and this prevents spurious identification of scaling properties in time-series analysis. It is well known that the power spectrum of an LFP signal scales as 1/f β with β > 0. Our results indicate the existence of a particular power spectrum scaling law 1/f β with β < 0 for low frequencies (f < 4 Hz) for both SC and EE rearing conditions. This type of scaling behavior is associated to the presence of anticorrelation in the corresponding LFP signals. Moreover, since EE is an experimental protocol based on the enhancement of sensorimotor stimulation, we study the possible effects of EE on the scaling properties of secondary motor cortex (M2) and primary visual cortex (V1). Notably, the difference between Hurst’s exponents in EE and SC for individual cortical regions (M2) and (V1) is not statistically significant. On the other hand, using the detrended cross-correlation coefficient, we find that EE significantly reduces the functional coupling between secondary motor cortex (M2) and visual cortex (V1). DOI: 10.1103/PhysRevE.91.012702

PACS number(s): 87.10.−e, 87.19.ll, 84.35.+i, 05.40.Jc

I. INTRODUCTION

Physiological recordings in a variety of biological systems are characterized by complex fluctuations exhibiting powerlaw scaling behavior [1–3]. The scaling properties of neural signals of humans and animals can be used to get information on the mechanisms modulating neural activities [4–9]. The existence of scaling implies the absence of characteristic time scales (scale-free) for the physiological signal of interest. Scale-free dynamics is a hallmark of long-range correlation (LRC) that indicates the presence of memory properties in the system. It is well known that the index called Hurst’s exponent (0 < H < 1) characterizes the nature of LRC [10]: if H > 1/2 the process presents positive LRC, otherwise if H < 1/2 the process exhibits anticorrelations. Usually, detrended fluctuation analysis (DFA) is employed to assess the existence of LRC in signals [11]. There are many examples of the application of DFA to neural recordings [4,5,12]. For instance, the correlation properties of filtered (10 and 20 Hz) magnetoencephalography and elecroencephalography (EEG) signals (amplitude envelopes) from normal human brain were studied and found to exhibit positive LRC [4]. Although less frequently, also anticorrelations have been observed in physical and biological complex systems [13–16]. However, as pointed out in several research works, spurious identification of LRC can occur [12,17]. For instance, it has been shown that the presence of trends or periodic components in a signal impacts the performance of the DFA method [12,18]. In particular, the existence of anticorrelations in a time

* Current address: Institute of Human Physiology, Catholic University of the Sacred Heart, 00168 Rome, Italy. † Corresponding author: [email protected]

1539-3755/2015/91(1)/012702(15)

series could be strongly biased and to assess their presence it is convenient to use an integrated approach [17]. In this work, the LRC properties of a signal will be studied by using three different methods simultaneously. Specifically, the Higuchi method [19], the DFA [11], and spectral analysis [20] will be used to estimate the fractal dimension, the Hurst exponent, and the scaling properties of the power spectrum, respectively. Since the values of these parameters are related to each other [17,21], we obtain a consistent scheme for the quantification of the LRC. Here we will analyze recordings of local field potentials (LFPs) from the secondary motor cortex (M2) and visual cortex (V1) of freely moving mice, reared in standard (SC) and environmental enrichment (EE) conditions. EE is an experimental protocol based on a complex sensorimotor stimulation that affects brain development, including corticocortical coupling [22]. It is worth remarking that our work fits well in the general framework concerning the dynamics of complex systems subject to external inputs [23]. Indeed, each neuron can be modelled as a nonlinear oscillator [24] and LFP recordings as a superposition of the overall synaptic activities generated by such networks of coupled nonlinear oscillators subject to external perturbations. Our attention will focus on the study of LRC properties of LFPs recorded from mice reared in SC and EE. Usually, electrophysiological recordings of LFPs are characterized by a power spectrum scaling law that “inversely” relates the magnitude of LFP power to frequency f , i.e., 1/f β , where β > 0 [25–27]. We found that LFPs recorded from the secondary motor cortex (M2) and visual cortex (V1) of freely moving mice, in SC and EE conditions, exhibit a power spectrum scaling law that “proportionally” relates the magnitude of LFP power to frequency f , i.e., 1/f β , where β < 0, for low frequencies 0.3 Hz < f < 4 Hz. This

012702-1

©2015 American Physical Society

VALLONE, CINTIO, MAINARDI, CALEO, AND DI GARBO

PHYSICAL REVIEW E 91, 012702 (2015)

particular type of scaling behavior is an evidence of the existence of anticorrelations in the LFP signals. On the other hand, we achieved the usual power spectrum scaling law for frequencies ranging in the interval 13 Hz < f < 30 Hz. To the best of our knowledge, this is the first time that this dual kind of power spectrum scaling law has been found in LFPs. Moreover, EE affects several electrophysiological and neuroanatomical parameters in the adult brain. In particular, EE modulates a series of factors involving synaptic plasticity [28]. It is known that the main source of LFPs are extracellular currents arising from the synaptic activity generated by thousands of coupled neurons [25]. Thus, the increased somatosensory stimulation experienced by an animal raised in EE could impact the dynamical properties of LFPs. In fact, in filtered (10 and 20 Hz) EEG signals (amplitude envelopes) from normal human brain, it was shown that LRC are attenuated by somatosensory stimulations [5]. However, to our knowledge, this is the first attempt aiming at investigating the LRC properties of LFP signals from two different cortical areas of freely moving mice previously reared in EE. We found no statistically significant difference between Hurst’s exponents in EE and SC for individual cortical regions (M2, V1). On the other hand, using the detrended cross-correlation coefficient (ρDCCA ) [29], we found that EE reduces the functional coupling between the secondary motor cortex (M2) and the visual cortex (V1) with respect to mice reared in SC. The paper is organized as follows. In Sec. II A we will introduce the basic properties of LRC and we will continue in Sec. II B with the discussion of the properties of the fractional Gaussian noise (FGn) and fractional Gaussian motion (FBm). Then, in Sec. II C we will describe the data and how they have been acquired and we will give a concise account of the methods that will be used to study the LFP signals. In Sec. III A we will refer to FGn models, with known Hurst’s exponents, to test the validity of the integrated approach. In Sec. III B we will discuss the main results obtained by studying the LRC properties of LFP recordings. Finally, conclusions will be presented in Sec. IV. II. THEORY AND ANALYSIS METHODS A. Long-range correlation (LRC)

The notion of long-range correlation (long memory) turned out to be helpful in different applications of time-series analysis. In particular, this concept was successfully employed to describe dynamical properties of signals from different sources: for instance, electrical activities recorded on the scalp [23,30], heart-rate variability [31,32], the sequence of base pairs in DNA [11,33], gait fluctuations [21,34], and financial volatility [35]. Intuitively, a stationary time series is characterized by long memory if observations that are far apart (in time or space) maintain some degree of correlation. The above definition can be stated formally as follows. Let X = (xi : i = 1,2, . . . ,∞) be a one-dimensional, stationary stochastic process in discrete time with zero-mean, finite variance σ 2 = E(xk2 ) < ∞, covariances γn = E(xk xn+k ), and correlations ρn = γn /σ 2 , n = 0, 1, 2, . . . ,∞ (where E means the  expected value). Now, let us define the random walk Rn = ni=1 xi and Var(Rn ) be the corresponding  variance. Then it can be shown that Var(Rn ) = σ 2 n[1 + 2 nk=1 (1 − k/n)ρn ] (by assuming that

the correlations depend only on the lag) [36]. Moreover, it can be shown [36] that if one assumes the summability of correlations, i.e., ∞ 

ρn < ∞,

(1)

n=0

then one concludes that the variance of the partial sums of the random walk grows linearly with the number of steps. On the other hand, if the correlations are not summable, i.e., (1) does not hold, then the variance of partial sum grows faster than linearly with n and it follows from Lamperti’s theorem (see, e.g., Theorem 2.1.1 in Ref. [37]) that convergence to the Brownian motion is impossible. On the basis of the above considerations a stochastic process is said to exhibit long-range dependence when the series in the left-hand side of (1) is divergent. That is the case if and only if there exists δ ∈ (0, 1) such that lim ρn ≈ n−δ .

n→∞

(2)

Thus, through the analysis of the summability properties of the correlations function  one can establish the existence, or not, of a time scale τ = ∞ s=0 ρs [38]. Therefore, when a long-range dependence is present, no finite time scale is provided and the processes generating the signal are scale free. B. FBm and FGn processes

A stochastic process Y = (yi : i = 1,2, . . . ,∞) is called self-similar if exists H such that for all c > 0 one has yk = c−H yc k , d

(3)

d

where by = one means the equality of distributions and H is called the Hurst exponent. Let Y be a zero-mean [E(yn ) = 0] Gaussian process with y0 = 0, self-similar with exponent H and with stationary increments. Then it can be shown that 0 < H < 1 [36,39]. In this case, the stochastic process Y = (yi : i = 1,2, . . . ,∞) is called (discrete) fractional Brownian motion (FBm). The process Y has stationary increments when, for all k  0, the distribution probability of (yk + h − yk ) is independent of h. Let Y be a FBm and xn = yn+1 − yn (n = 0,1, . . . ,∞) be the corresponding increment process. Then xn is a zero-mean Gaussian process with covariance given by γn = E(x0 xn ) =

σ2 (|n + 1|2 H − 2|n|2 H + |n − 1|2 H ), 2 (4)

where σ 2 = E(xn2 ) [36]. This stochastic process is called (discrete) fractional Gaussian noise (FGn) and is fully specified by the value of its variance (σ 2 ) and Hurst’s exponent (0 < H < 1) [40]. We will now discuss briefly some points that are relevant when the above theory is applied to analyze real-world signals. (i) When searching for long-range correlations in a stationary signal, the basic assumption is that it can be modelled with a FGn. Then the consistence of this hypothesis can be assessed by estimating the self-similarity parameter (H ),

012702-2

EXISTENCE OF ANTICORRELATIONS FOR LOCAL FIELD . . .

PHYSICAL REVIEW E 91, 012702 (2015)

which completely specifies the statistical properties of the FGn and of the corresponding FBm. (ii) From (4) it follows that the Hurst exponent determines the asymptotic behavior of the correlation ρn = γn /σ 2 for large lags:

with a sampling rate of 100 Hz as the differential between two adjacent electrode sites placed in the same cortical area, 50 000× amplified and 0.3- to 30-Hz band-passed. Before the analysis, all time series were visually inspected to confirm the absence of recordings artifacts. As discussed above, the correlation properties of LFP signals will be studied by using an integrated approach. In particular, we will implement specific algorithms to estimate the fractal dimension (H = 2 − D), the DFA exponent α (H = α), and the scaling property ) of the power spectrum and then we will check for (H = β+1 2 the consistence of the estimations of the Hurst exponents. Moreover, the DCCA will be used to study the coupling between the LFPs recorded from visual cortex (V1) and secondary motor cortex (M2). DCCA will be described in Sec. II F.

lim ρn ≈ [H (2 H − 1) n2 H − 2 ],

n→+∞

H = 1/2.

(5)

Hence, the FGn has positive correlations (or is persistent) whenever 1/2 < H < 1 and negative correlations (or is antipersistent) if 0 < H < 1/2. In the last case it can be shown  ρ = 0 [36]. Strictly speaking, a FGn with negthat ∞ n n=−∞ ative correlation cannot be thought as a long-range correlated process, because it does not satisfy Eq. (2). However, this process has a value of H = 1/2 and we will still refer to it as a long-range correlated process. It is important to stress that it  is very difficult to test the condition ∞ n=−∞ ρn = 0 on real data [17,36]. (iii) Let D be the Hausdorff dimension of the sample path of a FBm with Hurst’s exponent H . The parameters D and H describe different features of the FBm: D characterizes the local properties of the walk and H the global ones (correlation properties). As shown in Ref. [41], these two parameters are related by the following well-known relationship: D =2−H

(6)

that can be used in applications to get an estimation of the Hurst exponent starting from the fractal dimension of the corresponding FBm path. (iv) Similarly, the quantitative characterization of the scaling property of the power spectrum of a FGn can be used to estimate the corresponding Hurst’s exponent. Let X = (x(i) : i = 1, . . . ,∞) represent a FGn with Hurst’s exponent H and ρn be the corresponding autocorrelation function. Then the corresponding power spectrum is given by (Wiener-Khinchin’s theorem) SH (f ) = ρ0 + 2

∞ 

ρn cos(2π f n).

(7)

1. Fractal dimension

Fractal analysis offers a way of describing the extent of self-similarity in biological signals. Monofractal signals are characterized by a uniformly consistent scaling that is quantified by the fractal dimension D (or Hurst’s exponent H = 2 − D). In order to estimate the fractal dimension, the Higuchi method will be employed [19]. Let x(i),i = 1, . . . ,N be a given time series, and then the implementation of the above method is performed as follows: (i) Construct a new time series,     N −m k k , (10) xm : x(m) , x(m + k), . . . , x m + k m = 1,2, . . . ,k, where [ ] denotes the Gauss’ notation and both k and m are integers. (ii) Define the length of the curve xmk as follows: ⎛ N−m ⎞ [ k ]  N −1 |x(m + ik) − x[m + (i − 1)]k|⎠ N−m 2 Lm (k) = ⎝ k k i=1 (11)

n=1

Then it can be shown that the scaling behavior of SH (f ) for small frequencies is given by [17] SH (f ) ∼

1 , fβ

f → 0,

(8)

and consider the average length of the curve L(k) for a fixed time interval k [obtained by averaging over m = 1,2, . . . ,k of Lm (k)]. (iii) If x(i) is a signal with fractal dimension D, then L(k) ∼ k −D

where β = 2H − 1.

(9)

Consequently, the scaling properties of the power spectrum can be used to estimate the Hurst exponent of the corresponding FGn. It is worth remarking that the relation in Eq. (9) between Hurst’s exponent and the scaling coefficient β is valid only for small frequencies. C. Data description and analysis methods

LFP recordings were performed in awake, freely moving mice using the protocol described in Ref. [22]. Each data set consists of a bivariate time series representing the LFPs simultaneously recorded from visual cortex (V1) and secondary motor cortex (M2). Cortical LFP signals were acquired

(12)

and the curve has fractal dimension D. 2. Detrended fluctuation analysis

The DFA method was introduced by Peng et al. [11] and has been widely applied to characterize the correlation properties of time series [4,5,11,12,14,15,31–33]. The DFA method is used to estimate the Hurst exponent and the corresponding algorithm consists of the following steps: (i) The data series x(i),i = 1, . . . ,N is shifted by the mean x and integrated (cumulatively summed),  y(k) = ki=1 (x(i) − x ) and then segmented into windows of various size n. (ii) In each window the integrated data are locally fitted by a linear polynomial yn (k) and the mean-squared residual F (n)

012702-3

VALLONE, CINTIO, MAINARDI, CALEO, AND DI GARBO

(fluctuations) is calculated:   N 1  [y(k) − yn (k)]2 . F (n) =  N k=1

PHYSICAL REVIEW E 91, 012702 (2015)

(13)

(iii) Whenever a power-law scaling is present, then F (n) scales as F (n) ∼ nα .

(14)

If 0 < α < 1, then the signal x(i) can be modeled by means of a FGn process with H = α. On the other hand, when 1 < α < 2 the process is nonstationary and can be described by a FBm with H = α − 1.

us assume to know the power spectrum SH (f ) ∝ 1/f 2H −1 of the FGn. Then, we consider a sample path in the frequency domain, i.e., a sequence of complex numbers {sk } whose squared moduli are an ordered sequence of values of SH (f ). Then, by using the inverse discrete Fourier transform, we can obtain the sequence {xi }, i.e., the time-domain counterpart to {sk }. Because, by construction, the power spectrum of {xi } is SH (f ), and because autocorrelation and power spectrum form a Fourier pair, then {xi } is guaranteed to have the same autocorrelation function ρi as the FGn process. Consequently, the time series {xi } is a realization of such a process with a given H . F. Detrended cross-correlation analysis

D. Determination of scaling intervals and scaling exponents

Scaling intervals and scaling exponents of the quantities described by the Eqs. (8), (12), and (14) are estimated by looking at the local slopes θ ’s of both SH (f ), L(k) , and F (n) in double logarithmic plots [42,43]. The local slopes are calculated using a numerical derivative. For the DFA this reads θ (n)dfa =

ln F (n + 1) − ln F (n) . ln(n + 1) − ln(n)

(15)

To determine the scaling interval, we define the mean local slope and its standard error over an interval Iln(n),δn = [ln nl ,ln nr ] containing δn points as follows:

δn θ¯dfa

nr 1  δn ¯θdfa := θ (n)dfa , δn n=n l   nr 

1 δn 2  θ (n)dfa − θ¯dfa := √ . (δn − 1)δn n=nl

(16)

(17)

We choose as scaling interval the biggest region Iln(n),N among a set of intervals {Iln(n),δn }, where the ratio of the mean δn δn local slope θ¯dfa and its standard error θ¯dfa satisfies a requested precision , i.e., Iln(n),N := max

(18)

δn θ¯dfa  δn θ¯dfa

δn θ¯dfa , {Iln(n),δn } θ¯ δn dfa

∀ Iln(n),δn .

(19)

The DFA exponent and its error are associated with the mean local slope and its standard error over the scaling interval, i.e., N , α ≡ θ¯dfa N α ≡ θ¯dfa .

(20) (21)

The use of the local slope is useful to avoid bias due to the logarithmic transformation on the interested variables. The nonlinearity of the transformation hides significant curvatures in plots like ln(F (n)) against ln(n) [42,43]. E. Algorithm to generate self-similar processes

To synthesize a self-similar process with Hurst’s exponent H , the Fourier filtering method was used [44,45]. This is an approximated method to generate FGn. However, it works well and is largely employed in several applications [17,44–47]. Let

So far we have considered methods that allow to study the scaling properties of single signals. Signals from different cortical areas are univariate components of a complex system that are probably interacting. A natural step forward is to establish whether the time series are power-law cross-correlated, that is, the time lagged cross-correlation function between the pair of signals asymptotically takes the form of a power law. The DCCA [48] provides a tool for investigating the presence of a correlation across various time scales between cortical signals by extending the DFA technique. We assume to have two time series (xi , i = 1, . . . ,N) and (yi , i = 1, . . . ,N ) that are power-law autocorrelated with Hurst’s exponents H and G, respectively. We are interested in testing the presence of long-range cross-correlations, i.e., whether there exist some 0 < βxy < 1 and c > 0 such that, as k → ∞, c (22) E(xi yi+k ) ∼ βxy . k In the same way as for the long-range dependence, the information conveyed, in the limit of large time lags, by a cross-covariance function is connected to its Fourier transform (cross-power spectrum) as f → 0. For that reason an empirical time-lagged cross-covariance is unreliable as a guide to the analysis of power-law properties. Actually, the lowest-frequency part of its spectrum is often contaminated by deterministic trends. In precise analogy to DFA estimation of the Hurst parameter of a covariance function, it is established that there is a given notion of cross-correlation which is not influenced by deterministic nonstationarities. The correspond2 ing cross-covariance Fx,y function between the data series (xi ) and (yi ) is defined by taking into account that (i) the fluctuation analysis is carried out on the aggregate  series Rx (k) = ki=1 (x(i) − x ), k = 1, . . . , N and the analogous series Ry (k) for the y data; (ii) the aggregate signals are split into N − n overlapping boxes of size n; in the box j = 1, . . . , N − n the least-squares j j fit R˜ x (R˜ y ) of (xi ) ((yi )) is computed; (iii) for the detrending of the nonstationarities in each box j we compute the covariance of the residuals: 1  ˜j ˜j fx2y (n, j ) = (n−1) “box j  (Rx − Rx ) (Ry − Ry ) The cross-covariance function

012702-4

Fx2y (n)

N−n 1  2 = f (n, j ), N − n j xy

(23)

EXISTENCE OF ANTICORRELATIONS FOR LOCAL FIELD . . .

PHYSICAL REVIEW E 91, 012702 (2015)

TABLE I. Hurst’s exponents HD,α,β , presented in the form of mean value over 10 realizations and standard error of the mean, of artificial FGn signals with positive correlations (H = 0.7) and negative correlations (H = 0.3). Positive correlation

Negative correlation

FGn model

HD

Hα

Hβ

HD

Hα

Hβ

H ± SEH

0.68 ± 0.01

0.69 ± 0.01

0.69 ± 0.01

0.31 ± 0.01

0.31 ± 0.01

0.29 ± 0.01

defined using the aggregate signals, is the detrended crosscovariance between (xi ) and (yi ). It is a function of the scale n which is used for the box-splitting and for averaging the fluctuations around the time trends. In Ref. [48] it was established that whenever (xi ) and (yi ) are long-range crosscorrelated [see Eq. (22)], then the covariance between their aggregate signals is a scaling

with given Hurst’s exponent, and the corresponding results will be discussed before presenting the results of the analysis of LFP data. Finally, we will present the results obtained from the implementation of the DCCA coefficient ρDCCA (see Sec. II F). The DCCA coefficient ρDCCA quantifies the coupling between LFPs recorded from V1 and M2.

Fx2y (n) ∼ n2 λx y

A. Artificial signals

(24)

as n → ∞. The parameter λx y is the bivariate Hurst exponent. The covariance analysis associated with Fx2y is called DCCA. Let us consider cross-correlated series (xi ) and (yi ) with Hurst’s exponents H and G, i.e., Fx1 (n) ∼ nH and Fy1 (n) ∼ nG . In DCCA method one only focuses on the scaling of Fx2y . In Ref. [29], in analogy to Pearson’s correlation coefficient, a new parameter is introduced, ρDCCA (n) =

Fx2y (n) Fx1 (n) Fy1 (n)

(25)

as a starting point in testing for long-range interdependence. It is called DCCA coefficient. One can show that ρDCCA lies between −1 and 1 as for the standard variance-covariance approach [49]. The DCCA coefficient can be interpreted as a standard correlation coefficient with ρDCCA (n) = −1 for perfectly anticorrelated series, ρDCCA (n) = 1 for perfectly correlated series and ρDCCA (n) = 0 for uncorrelated processes. In Ref. [48] it is indicated that, in general, the exponent λx y in the scaling law of Fx2y [see Eq. (23)] approaches H +G for 2 large-enough n. Hence, the coefficient ρDCCA asymptotically becomes independent of n. This aspect makes apparent that F 2 provides less information than ρDCCA . In fact, the exponent λx y (i.e., F 2 ) only determines the long-range power-law correlations, while ρDCCA provides the level of cross-correlations [29]. III. RESULTS

In this section we will present and discuss the results of the analysis of LFP signals recorded from different cortical areas (M2 and V1) of awake, freely moving mice raised in either SC or EE conditions. In particular, it will be shown that the scaling properties of the LFP recordings are consistent with those of an anticorrelated FGn in both SC and EE mice. As discussed in Sec. II C, the study of the correlation properties of the LFP recordings was performed by using an integrated approach that was implemented by employing three different methods (fractal dimension, DFA, power spectrum), each one giving direct or indirect informations on the scaling properties of a signal. This approach prevents false identification of LRC and allows for three independent estimations of Hurst’s exponent H . All the employed methods were tested on artificial signals,

Artificial FGn signals will be used to test the algorithm for the estimation of the fractal dimension (see Sec. II C 1). To this aim, FGn signals with positive and negative correlation properties were generated by using the Fourier filtering method described in Sec. II E. In particular, 10 realizations (each containing N = 218 data points) of FGn with positive (H = 0.7) and negative (H = 0.3) correlations were considered. Then the fractal dimension of the corresponding FBm was calculated and Hurst’s exponent HD was estimated by means of Eq. (6). Each FBm is generated by xFBm (i) =

3 10 +i 

x(j ),

(26)

j =1

where x(j ) is one of the FGn realizations with zero mean and standard deviation 1. The value 103 is chosen to avoid statistical error on the first FBm value. The results are presented in the form of the mean value over 10 realizations and standard error of the mean. The Hurst exponents turns out to be HD = 0.68 ± 0.01 and HD = 0.31 ± 0.01 for the positively and negatively correlated signals, respectively (see Table I). As an example, in top-left panel of Fig. 1, the results corresponding to the estimate of the fractal dimensions D obtained with a single realization of the FGn are plotted. The local slopes give D = 1.32 and D = 1.69 for the positively and negatively correlated FBm, respectively. Second, the Hurst exponent was estimated by using the DFA method (see Sec. II C 2). The time window (n) ranges from nmin = 3 to nmax = 700, with n incremented by a unit step each time. The corresponding results are Hα = 0.69 ± 0.01 and Hα = 0.31 ± 0.01 (see Table I). In the top-right panel of Fig. 1 an example of application of the DFA technique on a single realization of FGn is shown. The corresponding values for the DFA exponent are α = 0.69 and α = 0.31 for the positively correlated and the anticorrelated FGn signals, respectively. Last, the estimation of Hurst’s exponents was obtained by using the scaling properties of the power spectrum [SH (f ) ∼ f1β ]. We estimated the power spectrum by averaging modified periodograms. We adopted the Welch procedure of average on nonoverlapped windows each one containing N = 1024 data point, using Hanning windowing. The results

012702-5

VALLONE, CINTIO, MAINARDI, CALEO, AND DI GARBO

PHYSICAL REVIEW E 91, 012702 (2015)

4

α=0.31 -8

0

ln[F(n)]

ln[〈L(k)〉]

α=0.7

-10

D=1.69

-4

-12

-8

D=1.32 -14

-12

0

2

4

ln[k]

6

8

1

2

3

4

5

ln[n]

6

7

7

9

8

β=-0.39

β=0.41

6

ln[S(f k )]

ln[S(f k )]

7

5

6

5

4 4

3

-8

-6

-4

ln[f k ]

-2

0

3 -8

-6

-4

ln[f k ]

-2

0

2

FIG. 1. Integrated approach applied to signals generated using the Fourier filtering method, with given Hurst’s exponents. In the top-left panel L(k) against k is plotted in a double-logarithm scale. In the top-right panel ln(F (n)), obtained through the DFA, is plotted against ln(n) in a double-logarithm scale. The corresponding DFA exponent α is in keeping with Hurst’s exponent estimated using the fractal dimension, i.e., α = 0.69 and α = 0.31. In the bottom panels the corresponding power spectra in a double logarithm scale are shown. For the case H = 0.7 and H = 0.3 the linear fitting gives β = 0.41 (bottom left) and β = −0.39 (bottom right), respectively.

are Hβ = 0.70 ± 0.01 and Hβ = 0.29 ± 0.01 (see Table I). In the bottom panels of Fig. 1 the opposite behavior, for low frequencies, between the power spectra of positively (β = 0.41) and negatively (β = −0.39) correlated FGn signals is clearly shown. The comparison of the mean values of Hurst’s exponents HD ,Hα ,Hβ was done by using one-way ANOVA (p = 0.05 as significance threshold). The corresponding results show that the mean values of Hurst’s exponents calculated with the three methods are not statistically different (p > 0.05 both for positively and negatively correlated FGn signals). These findings guarantee the consistence of the integrated approach and establish a safe ground for the analysis of LFP recordings. In several papers the scaling properties of neural data are investigated by analyzing the scaling properties of the envelope of the raw signal (or of a certain frequency band) [4,5,12]. However, if we are interested in assessing the scaling properties of the raw recording, how does this filtering change the scaling properties of the original signal? In other words, let xi be a signal with Hurst’s exponent H and yi be the corresponding envelope signal. Does the Hurst exponent of yi differ from H ? This topic is investigated below by using

artificial FGn signals and by extracting the corresponding envelope signal through the Hilbert transform. Let x be a signal in continuous time, and then the Hilbert transform of x is the limit  1 x(τ ) ˜ = lim dτ (27) x(t) δ→0 π |t−τ |>δ (t − τ ) defined almost everywhere [50]. It is a well-known theorem that x can be continu˜ = ously extended to the complex function z(t) = x(t) + i x(t) a(t)eiφ(t) that is analytic in the upper half of the complex plane. The function a(t) is called the envelope signal of x and φ(t) is its instantaneous phase. For a narrow-band signal the envelope is defined as ˜ 2 )1/2 a(t) = |z(t)| = (x(t)2 + x(t)

(28)

and the instantaneous phase is determined by φ(t) = ˜ arctan(x(t)/x(t)) [51]. Let now x(t) be an artificial FGn with a fixed value of H and a(t) be the corresponding envelope signal. Let us consider 10 realizations of x(t) for H = 0.3, and then the mean value of Hurst’s exponents of the corresponding envelope signals is Hαe = 0.51 ± 0.01. Similarly, for the case

012702-6

EXISTENCE OF ANTICORRELATIONS FOR LOCAL FIELD . . . -6

PHYSICAL REVIEW E 91, 012702 (2015)

power spectrum is shown. The power spectra of the signals used in the present work, in a double-logarithm diagram, show a common feature: increasing β < 0 and decreasing β > 0 behaviors, depending on the frequency intervals. The intervals are identified by the shaded areas and represent two well-defined physiological frequency bands, a low-frequency interval called a delta band, f < 4 Hz, and a beta band, 13–30 Hz (see right panel of Fig. 3). Regarding the scaling of the power spectrum in the Beta band, it is important to stress out that we cannot relate the scaling property quantified by Hurst’s exponent with the scaling coefficient β. The Eq. (9) that links the Hurst exponent with the scaling coefficient β makes sense only for small frequencies, that in our analysis are comprised within in the delta band, f < 4 Hz. However, we can estimate the scaling coefficient β without relating it to Hurst’s exponent. The results for secondary motor cortex (M2) + + = 1.12 ± 0.2 and βSC = 1.49 ± 0.2 for the EE and are βEE SC conditions, respectively. For V1 the corresponding scaling + + coefficients are βEE = 1.18 ± 0.2 and βSC = 0.73 ± 0.2. Concerning the study of the LRC properties, here it is assumed that each LFP recording can be thought of as a realization of a FGn. Consequently, the LRC properties of these signals will be characterized by applying the integrated approach discussed above. As in the case of FGn models, we first estimate the fractal dimension of the signal and then the corresponding Hurst’s exponent by using the relation H = 2 − D. The results for the Hurst exponents are represented as the mean value over the number of animals (eight for the SC group and six for the EE group) together with the standard error. Examples of fractal analysis are shown in the top panels of Figs. 4 and 6. The determinations of the scaling intervals and fractal dimensions have been done using a numerical derivative (see Sec. II D). Examples are shown in the top panels of Figs. 5 and 7, where the scaling regions fall in Iln(k) = [1.7, 3.7] and Iln(k) = [2.0, 4.2], respectively. The mean local slopes and their standard errors are equal to θfd ± θfd = −1.77 ± 0.01 and

α=0.51

ln[F(n)]

-8

-10

α=0.56

-12

-14 1

2

3

4

ln[n]

5

6

7

FIG. 2. DFA applied to envelope a(t) [see Eq. (28)] of the artificial FGn signals, generated using the Fourier filtering method with given Hurst’s exponents: specifically H = 0.3 (H = 0.7) for the top (bottom) curve. The fluctuation, obtained through the DFA, against n in a double-logarithm scale is shown. The corresponding DFA exponents α are in disagreement with Hurst’s exponents of the whole signals, i.e., α = 0.51 (α = 0.56) for the top (bottom) curve.

H = 0.7, the result is Hαe = 0.57 ± 0.01. In Fig. 2 the results for a typical realization of the FGn are reported. In conclusion, the above findings indicate that, if we are interested in the scaling properties of LFP recordings, then we have to consider the complete signal instead of a partial component such as the envelope. Last, the results discussed in this section show the importance of using consistent methodological approaches when studying the scaling properties of a signal. B. Scaling properties of LFP recordings

In the left panel of Fig. 3 the time course of a typical LFP recording is plotted, and in the right panel the corresponding

0.1

8 6

ln[S(f k )]

x(i)

0.05

0

4 2

-0.05

β0

-2 -0.1

0

10

5

2×10

5

3×10

5

i

-2

-1

0

1

2

3

4

ln[f k ]

FIG. 3. (Left) LFP signal from visual cortex (V1) of an animal in SC. (Right) In the double-logarithm diagram the spectrum of the LFP (left figure) has typical features. The plot qualitatively shows an increasing β < 0 and decreasing β > 0 behavior of the spectrum depending on frequency intervals (the frequency fk is expressed in Hz). Specifically, two frequency intervals where the spectrum can be well approximated by a linear fit can be identified. The intervals are identified by the shaded areas and represent two well-defined physiological frequency bands, the delta band f < 4 Hz and beta band 13–30 Hz. 012702-7

VALLONE, CINTIO, MAINARDI, CALEO, AND DI GARBO

PHYSICAL REVIEW E 91, 012702 (2015)

6

0

θ=-1.77

4

-2

θ(k) fd

ln[〈L(k)〉]

5

D=1.77 3 2

-4

1 0 1.5

2.5

-6

3.5

ln[k]

2

3

4

5

ln[k]

6

7

6

-3.5

4

α=0.22

2

ln[F(n)]

θ(n) dfa

-4

θ=0.22 0

-2

-4.5

-4

3

3.5

4

4.5

ln[n]

-6

5

7

10

6

5

θ(f k ) ps

ln[S(f k )]

-5

5

β=-0.7 4

3

4

5

ln[n]

6

θ=0.7

0

-5

-0.4

0

ln[f k ]

0.4

-10 -2

0.8

FIG. 4. Integrated approach applied to a LFP signal from visual cortex (V1) of an animal in SC (see Fig. 3). In the top panel is shown the length of the fractal curve, calculated through the Higuchi method, against k in a double-logarithm scale. The fractal dimension D turns out to be D ± D = 1.77 ± 0.01. In the middle panel the results obtained through DFA are shown. The corresponding value of the DFA exponent is α ± α = 0.22 ± 0.01 and is consistent with the value of Hurst’s exponent calculated through the fractal dimension HD = 2 − D = 0.23 ± 0.01. In the bottom panel, the power spectrum S(fk ) of the signal against fk is shown in a double-logarithm scale (the frequency fk is expressed in Hz). The corresponding value of β is β ± β = −0.7 ± 0.1 and the = 0.15 ± 0.05. corresponding value of Hurst’s exponent is Hβ = β+1 2

-1

ln[f k ]

0

1

FIG. 5. Determination of the scaling intervals and scaling exponents of fractal dimension, DFA and power spectrum techniques applied to a LFP signal from visual cortex (V1) of an animal reared in SC (see Fig. 4). In the top panel are shown the local slopes θ(k)fd of the length of the fractal curve L(k) . The scaling interval falls in the region Iln(k) = [1.7, 3.7], where the mean local slope is equal to θfd ± θfd = −1.77 ± 0.01. In the middle panel the results obtained for the DFA are shown. The corresponding scaling region is determined by the interval Iln(n) = [3.3, 4.3] in which the mean local slope assumes the value θdfa ± θdfa = 0.22 ± 0.01. In the bottom panel, the scaling region is Iln(fk ) = [−0.3, 0.65] and is characterized by values of the mean local slope θps ± θps = 0.7 ± 0.1 (the frequency fk is expressed in Hz).

012702-8

EXISTENCE OF ANTICORRELATIONS FOR LOCAL FIELD . . .

PHYSICAL REVIEW E 91, 012702 (2015)

5

0

θ=-1.83

4 -2

ln[〈L(k)〉]

3

θ(k) fd

D=1.83

2

-4

1 0.5 2

2.5

3

3.5

-6

4

2

4

ln[k] -3.0

6

ln[k]

6 4

α=0.17

2

ln[F(n)]

θ(n) dfa

-3.5

-4.0

θ=0.17 0

-2 -4

-4.5 3.0

3.5

4.0

4.5

-6

5.0

4

ln[n]

4.5

2

4

θ(f k ) ps

4

ln[S(f k )]

5

β=-0.3

3.5

3 -1.2

ln[n]

5

6

θ=0.3 0

-2

-0.8

-0.4

-4

0

ln[f k ]

-1.0

-0.5

0.0

0.5

ln[f k ]

FIG. 6. Integrated approach applied to LFP signal from secondary motor cortex (M2) of a mouse reared in EE. In the top panel is shown the length of the fractal curve, calculated through the Higuchi method, against k in a double-logarithm scale. The local slopes give the fractal dimension D ± D = 1.83 ± 0.01. In the middle panel the results obtained through DFA are shown. The value of the DFA exponent is α ± α = 0.17 ± 0.01 and is consistent with the value of Hurst’s exponent calculated through the fractal dimension HD = 2 − D = 0.17 ± 0.01. In the bottom panel, the power spectrum S(fk ) of the signal against fk is shown in a double-logarithm scale (the frequency fk is expressed in Hz). The corresponding value of β is β ± β = −0.3 ± 0.1 = and the corresponding value of Hurst’s exponent is Hβ = β+1 2 0.35 ± 0.05.

FIG. 7. Determination of the scaling intervals and scaling exponents of fractal dimension, DFA, and power spectrum techniques applied to a LFP signal from secondary motor cortex (M2) of a mouse reared in EE (see Fig. 6). In the top panel are shown the local slopes θ(k)fd of the length of the fractal curve L(k) . The scaling region is located in the interval Iln(k) = [2.0, 4.2] where the mean local slope assumes a value θfd ± θfd = −1.83 ± 0.01. In the middle panel the results obtained for the DFA are shown. The corresponding scaling interval falls within the region Iln(n) = [3.4, 4.5] determined by value of the mean local slope equal to θdfa ± θdfa = 0.17 ± 0.01. In the bottom panel the scaling region for the power spectrum S(fk ) of the signal is determined. The scaling interval turns out to be Iln(fk ) = [−1, 0] and is characterized a mean local slope θps ± θps = 0.3 ± 0.1 (the frequency fk is expressed in Hz).

012702-9

VALLONE, CINTIO, MAINARDI, CALEO, AND DI GARBO

PHYSICAL REVIEW E 91, 012702 (2015)

5

ln[〈L(k)〉]

4

D=1.5

3

2

2

2.5

3

3.5

ln[k]

4

-3

α=0.5

ln[F(n)]

-3.5

-4

-4.5

3

3.5

4

4.5

ln[n]

5

8

7

ln[S(f k )]

θfd ± θfd = −1.83 ± 0.01, respectively. For the M2 Hurst exponents are HD = 0.22 ± 0.03 and HD = 0.15 ± 0.03, for the EE and SC rearing conditions, respectively; on the other hand, for the visual cortex (V1) the corresponding Hurst exponents are HD = 0.25 ± 0.03 and HD = 0.21 ± 0.03 (see Tables II and III). These results indicate that LFP recordings from M2 and V1 cortices can be described as anticorrelated FGn signals. Each LFP recording has a proper scaling interval, reflecting the diversity of the neural signals, and from data analysis it follows that the typical time window where scaling occurs is around τ D ≈ 500 ms. The presence of scale-free dynamics suggests that LFP recordings could be the output of a collective activity of a large number of interacting components with different time scales. A possible biological interpretation for the presence of anticorrelated LRC in these neural signals will be discussed below. Let us consider now the estimation of the Hurst exponent by means of the DFA method (see Sec. II C 2). The time window n varies between nmin = 3 and nmax = 700, incrementing n by a unit step each time. Examples of the application of the DFA method to LFP signals are shown on the middle panels of Figs. 4 and 6. Furthermore, in the middle panels of Figs. 5 and 7 the local slopes are shown. For M2 we have Hα = 0.31 ± 0.04 and Hα = 0.24 ± 0.01 for the EE and SC rearing conditions, respectively. For V1 the corresponding Hurst exponents are Hα = 0.28 ± 0.02 and Hα = 0.29 ± 0.03 for the EE and SC rearing conditions, respectively (see Table III). These results indicate that LFP signals are characterized by anticorrelated LRC, in agreement to the previous results. The typical time window where scaling occurs is around τ α ≈ 535 ms and has the same order of magnitude found with fractal analysis. To further check the above results, the scaling properties of LFP signals were assessed by using the power spectrum analysis [see Eq. (9)]. Some examples are shown in the bottom panels of Figs. 4 and 6. The results for all animals are the following. M2: Hβ = 0.25 ± 0.05 and Hβ = 0.23 ± 0.03 for animals in EE and SC rearing conditions, respectively. V1: Hβ = 0.27 ± 0.03 and Hβ = 0.22 ± 0.03 for the EE and SC rearing conditions, respectively (see Tables II and III ). The frequency scaling intervals are within the physiological delta band (0.3 Hz < ν < 4 Hz). The width of a typical frequency scaling interval is around ν ≈ 1.5 Hz (see, for example, the bottom panel of Fig. 5). The amplitude of this frequency window is qualitatively in agreement with the mean value of the temporal scales arising from the fractal analysis and the D α DFA (1/ν ≈ τ +τ ). 2 An additional test on the LRC properties of the LFP signals was performed by implementing the method of surrogate data [20]. After a random shuffling of raw signals, the corresponding values of the Hurst exponents were consistent with H = 0.5 (see Table IV), i.e., the expected value for a memoryless signal. In Fig. 8 the corresponding results using Higuchi’s method, DFA, and power spectrum are shown. In conclusion, the above results suggest that all analyzed LFP recordings are characterized by anticorrelated LRC properties. To assess the consistence of the values of Hurst’s exponents, one-way ANOVA followed by post hoc Tukey’s test was employed and the corresponding results are reported in Tables V and VI. Using p = 0.05 as threshold of significance we found that Hurst’s exponents are consistent in EE for

β=0 6

5

4

0

0.4

ln[f k ]

0.8

FIG. 8. Integrated approach applied to a surrogate LFP signal from secondary motor cortex (M2) of an animal reared in EE. In the top panel the length of the fractal curve, calculated through the Higuchi method, against k in a double-logarithm scale, is shown. The fractal dimension D turns out to be D = 1.5. In the middle panel the results obtained through DFA are shown. In the bottom panel it is shown the power spectrum S(fk ) of the signal against fk in double-logarithm scale (the frequency fk is expressed in Hz).

M2 and V1 (see Table V). On the other hand, we found inconsistency in SC mice and for M2 for the pair (Hα ,HD ) (see Table VI). This inconsistency could be due to several factors, such as the finite length of the data, noisy fluctuations superposed on the power-law spectrum, some bias in the

012702-10

EXISTENCE OF ANTICORRELATIONS FOR LOCAL FIELD . . .

PHYSICAL REVIEW E 91, 012702 (2015)

TABLE II. Values of the fractal dimension D, DFA exponent α and power spectrum exponent β for all LFP recordings. Secondary motor cortex (M2) D

α

−β

D

α

−β

1 2 3 4 5 6

1.83 1.75 1.79 1.84 1.83 1.67

0.17 0.31 0.37 0.33 0.43 0.26

0.3 0.3 0.7 0.8 0.6 0.7

1.69 1.67 1.78 1.74 1.79 1.85

0.27 0.22 0.35 0.24 0.26 0.33

0.3 0.4 0.7 0.5 0.6 0.7

1 2 3 4 5 6 7 8

1.74 1.95 1.89 1.73 1.87 1.89 1.91 1.83

0.26 0.19 0.22 0.24 0.30 0.22 0.23 0.24

0.5 0.6 0.8 0.7 0.5 0.9 0.7 0.5

1.81 1.84 1.91 1.83 1.7 1.69 1.79 1.77

0.27 0.41 0.26 0.44 0.26 0.25 0.22 0.22

0.5 0.8 0.8 0.5 0.5 0.6 0.7 0.7

Animal

EE

SC

Visual cortex (V1)

DFA, or the presence of non-Gaussian noise. Nevertheless, our findings show that the dynamics of LFP signals is characterized by anticorrelated LRC. Notably, the impact of the EE protocol on the scaling properties of the signals is found to be not statistically significant. Indeed, using the t test and the significance threshold p = 0.05 the difference between Hurst’s exponents in EE and SC is not statistically significant (see Table VII). This result suggests that the scaling properties of LFP recordings are not affected by the stimulation caused by the EE rearing condition. Moreover, it is interesting to note that the EE rearing condition does not affect the scaling coefficient of the power spectrum in the Beta band. Indeed, using a t test we found no statistically significant difference between the scaling + + coefficients βEE and βSC for M2 and V1. On the other hand, using the DCCA coefficient (see Sec. II F), we found a significant effect of the EE on the functional coupling between M2 and V1. In principle, there is no cross-correlation whenever ρDCCA = 0. However, in practice that is true only for an infinitely long time series. For a finite one, even if cross-correlations are not present, ρDCCA would be probably some small nonzero value due to finite-size effects. In Fig. 9 the plots of the ρDCCA coefficients corresponding to the EE and SC conditions are shown. Furthermore, for all the individual data series (i.e., for LFP recordings obtained from each mouse), we generated a surrogate of its LFP, i.e., a completely decorrelated series produced by random and independent shuffling of the LFP data. We created two sets of decorrelated data series obtained

from the EE and SC conditions. In Fig. 9 the mean values of ρDCCA computed using the EE surrogates (EES ) and the SC ones (SCS ) and the corresponding standard errors are shown. These data show that, in the region n > 3000, the values of ρDCCA for the raw data are larger than those corresponding to surrogate data. Consequently, we can judge with a reasonable degree of confidence that the ρDCCA coefficients indicate actual nonzero correlations for both the EE and SC groups. The SC EE mean values ρDCCA (ρDCCA ) for all EE (SC) signals were estimated by averaging the corresponding ρDCCA (n) values EE = 0.17 ± 0.01 in the region n > 3000. The results are ρDCCA SC and ρDCCA = 0.24 ± 0.01 for the EE and SC conditions, respectively. The two mean values were compared by using the t test and the corresponding result indicates statistically significant difference with p < 0.05. The result is in agreement with Ref. [22], where, using the mutual information method, it was found that the coupling between the visual and motor areas is larger in SC-reared than in EE-reared mice. C. Physiological interpretation

Here we will provide a qualitative explanation for the anticorrelation properties of LFP signals. Furthermore, we will attempt to propose possible physiological mechanisms causing the lack of substantial effects on the scaling properties of individual cortical regions (M2, V1) of EE mice, along with the significant reduction in the functional coupling between M2 and V1. The anticorrelation shown by LFPs from mice reared in both SC and EE rearing conditions is a property related to the

TABLE III. Mean values and standard errors of Hurst’s exponents HD,α,β for LFP recordings from secondary motor cortex (M2) and visual . cortex (V1). The values of HD and Hβ where calculated by using HD = 2 − D and Hβ = 1+β 2 Secondary motor cortex (M2) Animal-LFP signals HEE ± SEHEE HSC ± SEHSC

Visual cortex (V1)

HD

Hα

Hβ

HD

Hα

Hβ

0.22 ± 0.03 0.15 ± 0.03

0.31 ± 0.04 0.24 ± 0.01

0.25 ± 0.05 0.23 ± 0.03

0.25 ± 0.03 0.21 ± 0.03

0.28 ± 0.02 0.29 ± 0.03

0.27 ± 0.03 0.22 ± 0.03

012702-11

VALLONE, CINTIO, MAINARDI, CALEO, AND DI GARBO

PHYSICAL REVIEW E 91, 012702 (2015)

s for surrogate data from secondary motor cortex (M2) and visual TABLE IV. Mean values and standard errors of Hurst’s exponents HD,α,β s s s cortex (V1). The values of HD , Hβ , and Hα are in agreement with Hurst’s exponent for a memoryless signal H = 0.5. The surrogate method is an additional validation for the presence of LRC in LFP signals.

Secondary motor cortex (M2) Animal-surrogate data s s HEE ± SEHEE s s HSC ± SEHSC

Visual cortex (V1)

HDs

Hαs

Hβs

HDs

Hαs

Hβs

0.51 ± 0.01 0.50 ± 0.01

0.50 ± 0.01 0.50 ± 0.01

0.50 ± 0.01 0.49 ± 0.01

0.51 ± 0.01 0.50 ± 0.01

0.50 ± 0.01 0.50 ± 0.01

0.51 ± 0.01 0.50 ± 0.01

information contents carried by the low-frequency components of the LFP signals. The information content of the high-frequency component of the LFP should be related to “local” properties of the dynamics of the cortical network [52]. On the other hand, low-frequency components are less local than high-frequency ones, since the size of signal-generation region surrounding an electrode is larger and the LFP generated by a synaptically activated population spreads further outside the population edge due to volume conduction [53,54]. Since the content of LFP signals in low frequencies carries nonlocal information, the anticorrelated scaling properties of LFP recordings found at low frequencies (delta band) could be a consequence of the existence of a negative feedback between spatially separated cortical regions. As pointed out recently, the brain is thought to be organized into distinct dynamically anticorrelated networks [55]. This implies that tasks or stimuli determine an increase of the activity in regions supporting the processing of such inputs and a decreased activity in areas supporting unrelated or irrelevant processes

[55]. In addition, in Ref. [56] it is stressed that in the brain, even at rest, large regions activate in bulk at the same time as other regions deactivate. This balance is the only possibility to avoid both total quiescence, in which the brain is shutdown, and massive excitation, in which the entire cortex is fired up. However, the knowledge of the mechanisms responsible for the maintenance of stability it still incomplete. Moreover, it is strongly conceivable that the presence of anticorrelations in biological systems is a hallmark of adaptivity, a property of enormous importance for biological systems (see Refs. [14,57]). Qualitatively, if the system that generates an anticorrelated signal is affected by an external perturbation, that signal will continue to maintain smaller fluctuations, reflecting an internal stability of the system. Intuitively, a negative LRC related to the low-frequency component of a neural signal is expected, instead of a positive one, which could not preserve the internal stability of a neural system subject to an external perturbation. A similar concept can be derived from the experimental observations of homeostatic regulation

TABLE V. Comparison between Hurst’s exponents of LFP recordings of animals reared in EE estimated through fractal dimension, DFA, and power spectrum. At the entry i,j = D,α,β, the table displays the values pi,j that are obtained by one-way ANOVA followed by post hoc Tukey’s test comparing Hi ,Hj . Above the diagonal, Hurst’s exponents for visual cortex signals are compared, whereas, below the diagonal, the Hurst exponents for motor cortex signals are compared. The statistical significance threshold was p = 0.05.

TABLE VI. Comparison between Hurst’s exponents of LFP recordings of animals reared in SC estimated through fractal dimension, DFA, and power spectrum. At the entry i,j = D,α,β, the table displays the values pi,j that are obtained by one-way ANOVA followed by post hoc Tukey’s test comparing Hi ,Hj . Above the diagonal, Hurst’s exponents for visual cortex signals are compared, whereas, below the diagonal, the Hurst exponents for motor cortex signals are compared. The statistical significance threshold was p = 0.05.

EE

SC

Tukey

HD

HD

Hβ

0.804

Tukey

0.703

0.867

HD

0.952

0.518 M2

HD

Hα

0.03

Hβ

0.067

Hα

Hβ

0.127

0.995

0.131

0.92 M2

012702-12

V1

0.22

Hβ

V1

Hα

Hα

EXISTENCE OF ANTICORRELATIONS FOR LOCAL FIELD . . .

PHYSICAL REVIEW E 91, 012702 (2015)

TABLE VII. Statistical comparison between Hurst’s exponents in EE and SC. In the table are reported the values of pD,α,β obtained from in the statistical comparison of HD,α,β for EE and SC, calculated from fractal dimension, DFA, and power spectrum, respectively. The statistical significance threshold was p = 0.05. Secondary motor cortex (M2)

HD,α,β (EE) − HD,α,β (SC)

pD

pα

pβ

pD

pα

pβ

0.123

0.053

0.638

0.325

0.748

0.249

of neuronal networks, which allow the maintenance of nearly constant firing rates despite variations in afferent activity [58]. In particular, stabilization of neuronal activity may be achieved through a dynamic adaptation of the modification threshold for long-term potentiation (LTP) and long-term depression (LTD), depending on the preexisting excitability state [59]. For instance, a high level of previous activity slides to the right the modification threshold for LTP, favoring the induction of LTD, while a sustained reduction in postsynaptic activity has opposite effects, promoting the induction of LTP. These mechanisms may underlie the ability of cortical networks to keep neuronal activity within a functional dynamic range, thus preserving their stability regardless of external perturbations. It is conceivable that these properties may leave a mark at the level of LFPs and translate into a negative LRC. Moreover, in a model proposed by Levina et al. [60] to explain neuronal avalanches, electrical activity patterns of cortical networks which have been observed both in vitro [61] and in vivo [62], space and time distributions of avalanches display scale invariance. In Ref. [63], where avalanches are assumed to be a self-organized and quasicritical process [64], the Levina model is used to compute synaptic currents that are integrated into a LFP signal. The corresponding power spectrum possesses a maximum at a frequency fM and increases in the interval [0,fM ] [63]. Therefore, we found a possible model providing a dynamical mechanism which yields scale invariance (and the concomitant power-law distributions) and anticorrelations in the low-frequency band. In addition, evidences of the increase in the power spectrum of LFP recordings for increasing frequencies (theta and delta bands) come from experimental results from the left medial prefrontal cortex of awake rats [65]. We now discuss the possible physiological underpinnings of the effects of EE that result in a significantly lower functional coupling between motor and visual cortices, while leaving the LRC of these individual regions unaffected. EE is known to affect several electrophysiological and neuroanatomical parameters in the adult brain. Specifically, EE modulates a series of factors involved in the regulation of synaptic plasticity, including neurotrophic factors, neuromodulators, NMDA receptors, and inhibitory circuitry [28]. In particular, experiments in the visual system have clearly shown that a reduction of intracortical inhibition is a major mediator of the effects of EE [66,67]. Indeed, administration of the GABA-A receptor antagonist diazepam blocks the plasticizing effects of EE, indicating that EE acts by preserving low, juvenile levels of GABA-A-mediated inhibition into adulthood [67]. These data may be used to explain the lower functional coupling between V1 and M2 found in our analysis and the lack of effect of EE on the scaling properties of LFPs recorded in individual

cortical areas. The GABA-A receptor is a ionotropic receptor which mediates fast inhibitory neurotransmission. Thus, a reduction of GABA-A-mediated inhibition following EE is expected to reduce the correlation level of V1 and M2, as the LFPs of these two individual regions may oscillate with greater independence, thus lowering their coupling. Indeed, lowering the GABA-A-mediated inhibition impacts the firing activities of the pyramidal cells by decreasing their synchrony, thus making them more independent from each other [68–70]. At a network level this mechanism implies that the corresponding LFP signals fluctuate with lower correlation levels. At the same time, given the time scale of GABA-A-mediated inhibition (about 10 ms) one should not expect an effect on the LRC (time scale in the order of 100 ms) of signals recorded in individual regions. Thus, the impact of EE on GABA-Amediated neurotransmission predicts an effect on the coupling of interconnected networks but not on the degree of LRC in single cortical areas, in complete agreement with our experimental findings. Moreover, recent experimental results [71,72] suggest that an enhancement of visual stimulations (like in the EE condition) should not affect the scaling of low-frequency component of the LFP spectrum. In Refs. [71] the low-frequency component (0–4 Hz) of LFP recordings from the V1 cortex of a monkey during visual stimulations was not affected. Furthermore, in Refs. [72] LFP 0.7 SC EE SCS EES

0.6 0.5 0.4

ρDMCA

t test

Visual cortex (V1)

0.3 0.2 0.1 0 -0.1 0

1000

2000

3000

4000

5000

n

FIG. 9. Plots of the DCCA coefficients of representative SC and EE animals. In both cases ρDCCA assumes nearly constant values for n > 3000. The plots of the mean ρDCCA ’s computed from the surrogate data of the LFP signals are also shown. For each mean ρDCCA the corresponding standard errors (shadowed area) are displayed. As expected, the surrogate DCCA coefficients oscillate near zero for large-enought n and are much smaller than the asymptotic value for the corresponding LFP signals.

012702-13

VALLONE, CINTIO, MAINARDI, CALEO, AND DI GARBO

PHYSICAL REVIEW E 91, 012702 (2015)

recordings were performed in the V1 area of a macaque during presentation of grating patterns of increasing contrast. Also in this case, the low-frequency part of the spectrum did not change with contrast. In conclusion, these results on the impact of external stimulations on different frequency bands of LFP signals are in agreement with our findings.

The results discussed in this paper suggest the existence of anticorrelation properties of the LFP signals recorded from secondary motor cortex (M2) and visual cortex (V1) of freely moving mice reared in either environmental enrichment (EE) or standard condition (SC). The presence of anticorrelations was assessed by using three different techniques (fractal dimension, DFA, and power spectrum) that implement an integrated approach. We found no difference in the scaling

properties of LFP recordings when comparing EE and SC and this suggests that the scaling properties of the signal are not affected by the rearing condition. On the other hand, we found a significant effect of EE on the functional coupling between M2 and V1. Our results show that sensory stimulation could affect the electrical activities of neural networks in different ways. Therefore, our study could provide additional insight into basic fundamental questions about the mechanism of functional segregation and integration in the nervous system [73]. Moreover, the increasing β < 0 and decreasing β > 0 behaviours of power spectra, depending on frequency intervals, could have implications in modeling LFPs as well as self-organized criticality for brain models. Finally, our findings could be potentially relevant in clinical protocols (e.g., poststroke neurorehabilitation, Alzheimer’s disease) because they contribute to the understanding of the mechanisms driving the regulation of the neural activity evoked by sensorimotor stimulation.

[1] J. B. Bassingthwaighte, L. S. Liebovitch, and B. J. West, Fractal Physiology (Oxford University Press, New York, 1994). [2] P. Bak, How Nature Works: The Science of Self Organized Criticality (Springer, Berlin, 1996). [3] B. J. West, Front. Physiol. 1, 12 (2010). [4] K. Linkenkaer-Hansen, V. V. Nikouline, J. M. Palva, and R. J. Ilmoniemi, J. Neurosci. 21, 1370 (2001). [5] K. Linkenkaer-Hansen, V. V. Nikulin, J. M. Palva, K. Kaila, and R. J. Ilmoniemi, Eur. J. Neurosci. 19, 203 (2004). [6] B. He, J. Zempel, A. Z. Snyder, and M. E. Raichle, Neuron 66, 353 (2010). [7] D. Fraiman and D. R. Chialvo, Front. Physiol. 3, 307 (2012). [8] S.-S. Poil, R. Hardstone, H. D. Mansvelder, and K. LinkenkaerHansen, J. Neurosci. 32, 9817 (2012). [9] J. M. Palva, A. Zhigalov, J. Hirvonen, O. Korhonen, K. Linkenkaer-Hansen, and S. Palva, Proc. Natl. Acad. Sci. USA 110, 3585 (2013). [10] B. Mandelbrot, The Fractal Geometry of Nature (Times Book, New York, 1982). [11] C.-K. Peng, S. V. Buldyrev, S. Havlin, M. Simons, H. E. Stanley, and A. L. Goldberger, Phys. Rev. E 49, 1685 (1994). [12] R. Hardstone, S. Poil, G. Schiavone, R. Jansen, V. V. Nikulin, H. D. Mansvelder, and K. Linkenkaer-Hansen, Front. Physiol. 3, 450 (2012). [13] P. M. Richards, Phys. Rev. B 16, 1393 (1977). [14] C.-K. Peng, J. Mietus, J. M. Hausdorff, S. Havlin, H. E. Stanley, and A. L. Goldberger, Phys. Rev. Lett. 70, 1343 (1993). [15] S. Bahar, J. W. Kantelhardt, A. Neiman, H. H. A. Rego, D. F. Russell, L. Wilkens, A. Bunde, and F. Moss, Europhys. Lett. 56, 454 (2001). [16] C. Pan, B. Zheng, Y. Wu, Y. Wang, and X. Tang, Phys. Lett. A 329, 130 (2004). [17] G. Rangarajan and M. Ding, Phys. Rev. E 61, 4991 (2000). [18] K. Hu, P. C. Ivanov, Z. Chen, P. Carpena, and H. E. Stanley, Phys. Rev. E 64, 011114 (2001). [19] T. Higuchi, Physica D 31, 277 (1988). [20] H. Kantz and T. Schreiber, Nonlinear Time Series and Analysis (Cambridge University Press, New York, 1997). [21] C. Heneghan and G. McDarby, Phys. Rev. E 62, 6103 (2000).

[22] A. Di Garbo, M. Mainardi, S. Chillemi, L. Maffei, and M. Caleo, PLoS ONE 6, e25285 (2011). [23] D. R. Chialvo, Nat. Phys. 6, 744 (2010). [24] A. L. Hodgkin and A. F. Huxley, J. Physiol. 117, 500 (1952). [25] G. Buzs´aki, C. Anastassiou, and C. Koch, Nat. Rev. Neurosci. 13, 407 (2012). [26] J. Milstein, F. Mormann, I. Fried, and C. Koch, PLoS ONE 4, e4338 (2009). [27] K. Miller, L. Sorensen, J. Ojemann, and M. den Nijs, PLoS Computat. Biol. 5, 1000609 (2009). [28] A. Sale, N. Berardi, and L. Maffei, Physiol. Rev. 94, 189 (2014). [29] G. Zebende, Physica A 390, 614 (2011). [30] R. Hwa and T. Ferree, Physica A 338, 246 (2004). [31] C.-K. Peng, S. Havlin, H. E. Stanley, and A. L. Goldberger, Chaos 5, 82 (1995). [32] M. Barbi, S. Chillemi, A. Di Garbo, R. Balocchi, C. Carpeggiani, M. Emdin, C. Michelassi, and E. Santarcangelo, Chaos Soliton. Fract. 9, 507 (1998). [33] S. V. Buldyrev, A. L. Goldberger, S. Havlin, R. N. Mantegna, M. E. Matsa, C.-K. Peng, M. Simons, and H. E. Stanley, Phys. Rev. E 51, 5084 (1995). [34] J. M. Hausdorff, S. L. Mitchell, R. Firtion, C. K. Peng, M. E. Cudkowicz, J. Y. Wei, and A. L. Goldberger, J. Appl. Physiol. 82, 262 (1997). [35] N. Vandewalle and M. Ausloos, Physica A 246, 454 (1997). [36] J. Beran, Statistics for Long-Memory Processes (Chapman & Hall, New York, 1994). [37] P. Embrechts and M. Maejima, Self-Similar Processes (Princeton University Press, Princeton, NJ, 2002). [38] A. Montagnini, P. Allegrini, S. Chillemi, A. Di Garbo, and P. Grigolini, Phys. Lett. A 244, 237 (1998). [39] B. B. Mandelbrot and J. W. V. Ness, SIAM Rev. 10, 422 (1968). [40] G. Samorodnitsky and M. Taqqu, Stable Non-Gaussian Random Processes: Stochastic Models with Infinite Variance (Chapman & Hall, New York, 1994). [41] K. Falconer, Fractal Geonetry, Mathematical Foundations and Applications (Wiley, Chichester, 1990). [42] R. Bryce and K. Sprague, Sci. Rep. 2, 315 (2012).

IV. CONCLUSIONS

012702-14

EXISTENCE OF ANTICORRELATIONS FOR LOCAL FIELD . . .

PHYSICAL REVIEW E 91, 012702 (2015)

[43] Y.-H. Shao, G.-F. Gu, Z.-Q. Jiang, W.-X. Zhou, and D. Sornette, Sci. Rep. 2, 835 (2012). [44] H. Peitgen and D. Saupe, The Science of Fractal Images (Springer, New York, 1988). [45] J. Feder, Fractals (Plenum Press, New York, 1988). [46] C.-K. Peng, S. Havlin, M. Schwartz, and H. E. Stanley, Phys. Rev. A 44, R2239 (1991). [47] S. Prakash, S. Havlin, M. Schwartz, and H. E. Stanley, Phys. Rev. A 46, R1724 (1992). [48] B. Podobnik and H. Stanley, Phys. Rev. Lett. 100, 084102 (2008). [49] B. Podobnik, Z. Q. Jiang, W. X. Zhou, and H. E. Stanley, Phys. Rev. E 84, 066118 (2011). [50] E. Titchmarsh, Introduction to the Theory of Fourier Integrals (Oxford University Press, New York, 1948). [51] B. Boashash, Proc. IEEE 80, 520 (1992). [52] C. B´edard, H. Kroger, and A. Destexhe, Phys. Rev. E 73, 051911 (2006). [53] A. Sirota and G. Buzs´aki, Thalamus. Relat. Syst. 3, 245 (2005). [54] S. Leski, H. Linden, T. Tetzlaff, K. Pettersen, and G. Einevoll, PLoS Comput. Biol. 9, e1003137 (2013). [55] M. D. Fox, A. Z. Snyder, J. L. Vincent, M. Corbetta, D. C. Van Essen, and M. E. Raichle, Proc. Natl. Acad. Sci. USA 102, 9673 (2005). [56] D. Fraiman, P. Balenzuela, J. Foss, and D. R. Chialvo, Phys. Rev. E 79, 061922 (2009). [57] M. I. Rabinovich, P. Varona, A. I. Selverston, and H. D. I. Abarbanel, Rev. Mod. Phys. 78, 1213 (2006). [58] G. Turrigiano, Annu. Rev. Neurosci. 34, 89 (2011).

[59] W. Abraham, Nat. Rev. Neurosci. 9, 387 (2008). [60] A. Levina, J. M. Herrmann, and T. Geisel, Nat. Phys. 3, 857 (2007). [61] J. M. Beggs and D. Plenz, J. Neurosci. 23, 11167 (2003). [62] E. Gireesh and D. Plenz, Proc. Natl. Acad. Sci. USA 105, 7576 (2008). [63] J. A. Bonachela, S. de Franciscis, J. J. Torres, and M. A. Munoz, J. Stat. Mech. (2010) P02015. [64] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59, 381 (1987). [65] S. A. Neymotin, H. Lee, E. Park, A. A. Fenton, and W. W. Lytton, Front. Computat. Neurosci. 5, 19 (2011). [66] A. Sale, J. Vetencourt, P. Medini, M. Cenni, L. Baroncelli, R. De Pasquale, and L. Maffei, Nat. Neurosci. 10, 679 (2007). [67] F. Greifzu, J. Pielecka-Fortuna, E. Kalogeraki, K. Krempler, P. Favaro, O. Schl¨uter, and S. L¨owel, Proc. Natl. Acad. Sci. USA 111, 1150 (2014). [68] S. Cobb, E. Buhl, K. Halasy, O. Paulsen, and P. Somogyi, Nature 378, 75 (1995). [69] C. McBain and A. Fisahn, Nat. Rev. Neurosci. 2, 11 (2001). [70] J. Hjorth, K. T. Blackwell, and J. H. Kotaleski, J. Neurosci. 29, 5276 (2009). [71] A. Gail, H. Brinksmeyer, and R. Eckhorn, Cereb. Cortex 14, 300 (2004). [72] A. Henrie and R. Shapley, J. Neurophysiol. 94, 479 (2005). [73] G. Tononi, O. Sporns, and G. M. Edelman, Proc. Natl. Acad. Sci. USA 91, 5033 (1994).

012702-15