## Abstract

The phenomenon of synchronization between two or more areas of the brain coupled asymmetrically is a relevant issue for understanding mechanisms and functions within the cerebral cortex. Anticipated synchronization (AS) refers to the situation in which the receiver system synchronizes to the future dynamics of the sender system while the intuitively expected delayed synchronization (DS) represents exactly the opposite case. AS and DS are investigated in the context of causal information formalism. More specifically, we use a multi-scale symbolic information-theory approach for discriminating the time delay displayed between two areas of the brain when they exchange information.

## 1. Introduction

Brain function relies on the ability of neurons to communicate with each other. Interneuronal communication primarily takes place at synapses, where information from one neuron is conveyed to another neuron. In mammals, electrical transmission mediates reciprocal synchronizing interaction between neurons as well as ‘feed-forward’ excitation. The mammalian neocortex and hippocampus often generate synchronized, rhythmic patterns of activity that vary with behavioural state [1]. Synchronization between cortical neurons is therefore an astonishing collective phenomenon in the brain. Indeed, one of the main experimental challenges in neuroscience during the last 30 years was to demonstrate that distributed neural populations in the visual cortex process information in a synchronized way. The visual cortex is composed of a large number of areas, which contain neurons that are tuned to different visual features. Temporally synchronized activity of individual neuronal pairs within the visual cortex has been investigated in many laboratories since the early 1980s [2–6], most often with the motivation to reveal functional coupling between cells. This idea is in general agreement with the hypothesis that neighbouring cells with similar functional properties are tightly coupled to form a neural group. The concept that whole groups of neighbouring neurons could discharge synchronously in response to the same visual object has been attracting neuroscientists’ attention for many years [7]. Evidence has in fact been obtained that many neurons within a column of cat visual cortex can engage in a state of highly synchronous activity in response to an optimally oriented moving light bar [8,9].

Interestingly, local field potential (LFP) oscillations can only be observed when many neurons fire in synchrony, since otherwise the individual neurons’ electric fields would simply cancel out. Thus, the occurrence of high-frequency LFP oscillations demonstrates that local synchrony is generated with high temporal precision, although the LFP arises from extracellular dendritic fields rather than spike firing. Synchronization across neurons has been extensively investigated in the brain, where it has been hypothesized to underlie neurocognitive phenomena such as binding [10], temporal coding [11], spatial attention [12] and other higher cognitive functions [13]. As neural integration is usually understood as the algebraic representation and summation of excitatory or inhibitory postsynaptic potentials, which govern the potential for firing in the postsynaptic neuron, synchronization of presynaptic inputs can be thought of as a mechanism of neural integration. The neuronal activity of a single neuron or a group of neurons depends on intrinsic biophysical properties, and the interactions between different neuronal ensembles [14]. Importantly, the parietal and the motor cortex hold similar organizational principles to the visual cortex, and also consist of numerous areas. Any cerebral activity involves large numbers of areas and coordinated activity between neurons can be present in these areas. This coordination has been investigated by studying the synchronization between field potentials, which reflects the average activity of large groups of neurons in the vicinity of a recording electrode [15]. The strength of coupling between transcortically recorded field potentials in different cortical areas changes dynamically during the performance of a behavioural task. Synchronization has been found between areas of the visual and parietal cortex, and between areas of the parietal and motor cortex, in awake cats and monkeys [16]. Therefore, synchronization on a fine temporal scale appears to be a natural mechanism for the integration and coordination of neuronal activity between different brain regions.

When the connectivity between two regions of the brain (or in general two dynamical systems) is such that one of the regions (the sender) strongly influences the other one (the receiver), a positive time lag is often expected. In this situation, if the two regions synchronize the state is called delayed synchronization (DS); the sender’s dynamic determines the receiver’s one, in agreement with the common intuitive understanding. In the reverse situation, when the influence is transmitted from a sender to a receiver but the receiver’s dynamics leads the sender’s one in time, anticipated synchronization (AS) can occur [17,18]. This type of counterintuitive synchronization was originally found for two identical, autonomous dynamical systems coupled in a sender–receiver configuration. The only difference between the two systems is that the receiver was also subject to a negative delayed feedback [17,19]. The presence of this feedback, a type of ‘memory term’, enables the existence of a solution in which the activity of the receiver system at a time *t* is the same as that of the sender system at a time *t*+*τ* in the future, where *τ* is the feedback time. In other words, the receiver predicts the sender’s future dynamics. Interestingly, this can happen even for systems operating in a chaotic regime. AS has been shown to be stable in a wide variety of systems, ranging from purely academic theoretical studies to experimental works [20–27]. Ciszak *et al.* [24] were the first to find AS in neuronal models. In their study, the authors considered two unidirectional coupled dynamical systems, described by the FitzHugh–Nagumo and Hodgkin–Huxley neuronal models. Both systems were driven by a common white noise, and the receiver contained a negative delay feedback term. Interestingly, it was later shown in neuronal models that the delayed feedback term in the receiver system could be replaced by a dynamical inhibitory loop mediated by chemical synapses [28], transferring the concept to a more realistic situation. Under this condition, AS was found in small neuronal circuits in the presence of an interneuron [28], as well as in neuronal populations [18]. In the latter, it was shown that the model reproduced delay times, as well as coherence and Granger causality spectra, obtained from cortical data of monkeys performing a visual discrimination task. More recently, AS was also observed for other neuronal models and found to depend on the synaptic delays [29] and depolarization parameters [30].

In this paper, we use data from the model describing the dynamics between two brain regions presented in [18] together with experimental data from the sensorimotor cortex recorded between different cortical areas [14], which exhibit both DS and AS, in order to quantify the ‘relative synchronization phase’ between them. We estimate the time delay between different brain areas using subtle measures accounting for the nonlinear dynamic effects of the temporal signal: Shannon entropy [31,32] and the Martín–Platino–Rosso (MPR) statistical complexity [31,32] within the multi-scale entropy–complexity causality plane [33].

## 2. Information-theory quantifiers

### (a) Correlational structure of the time series

Sequences of measurements constitute the basic elements for the study of complex systems. These sequences are commonly called time series, and one should judiciously extract information on the dynamical systems under study. An information-theory quantifier can be defined as a measure that is able to characterize some property of the probability distribution associated with an observable or measurable quantity (i.e. membrane potential). Given a time series , a set of *M* measures of the observable and the associated probability distribution function (PDF), given by *P*≡{*p*_{j}; *j*=1,…,*N*} with and *N* the number of possible states of the system under study, Shannon’s logarithmic information measure [34] is defined by
2.1This functional is equal to zero when we are able to predict with full certainty which of the possible outcomes *j*, whose probabilities are given by *p*_{j}, will actually take place. Our knowledge of the underlying process, described by the probability distribution, is maximal in this instance. By contrast, this knowledge is commonly minimal for a uniform distribution *P*_{e}={*p*_{j}=1/*N*, ∀*j*=1,…,*N*}.

The Shannon entropy *S* is a measure of ‘global character’ that is not too sensitive to strong changes in the PDF taking place in a small region. The degree of structure present in a process is not quantified by randomness measures and, thus, measures of statistical complexity are needed to gain a better understanding of time series. The opposite extremes of perfect order and maximal randomness are too simple to describe as they do not have any structure, and complexity should be zero in both cases. It is important to point out that, at any given distance from these extremes, a wide range of possible degrees of physical structure can exist. The complexity measure allows us to quantify the different types of structures [35]. Here, we consider the MPR statistical complexity [36] as it is able to quantify critical details of dynamical processes underlying the dataset.

Based on the seminal notion advanced by López-Ruiz *et al.* [37], the MPR statistical complexity measure is defined through the product
2.2of the normalized Shannon entropy
2.3with , () and the disequilibrium defined in terms of the Jensen–Shannon divergence. That is,
2.4with
2.5where the above-mentioned Jensen–Shannon divergence and *Q*_{0}, a normalization constant (), are equal to the inverse of the maximum possible value of . This value is obtained when one of the components of *P*, say *p*_{m}, is equal to 1 and the remaining *p*_{j} are equal to zero. The Jensen–Shannon divergence, which quantifies the difference between two (or more) probability distributions, is especially useful to compare the symbolic composition between different sequences [38]. Note that the above-introduced statistical complexity measure depends on two different probability distributions, the one associated with the system under analysis, *P*, and the uniform distribution, *P*_{e}. Furthermore, it was shown that, for a given value of *H*_{S}, the range of possible *C*_{JS} values varies between a minimum and a maximum , restricting the possible values of the statistical complexity measure in a given entropy–complexity plane [39]. Thus, it is clear that important additional information related to the correlational structure between the components of the physical system is provided by evaluating the statistical complexity measure. In order to calculate the two information-theory-derived quantifiers mentioned previously, a probability distribution should be estimated from the time series of the system. Bandt & Pompe (BP) [40] introduced a successful methodology for the evaluation of the PDF associated with scalar time-series data using a symbolization technique.

The pertinent symbolic data are (i) created by ranking the values of the series and (ii) defined by reordering the embedded data in ascending order, which is tantamount to a phase-space reconstruction with embedding dimension (pattern length) *D* and time lag *τ*. In this way, it is possible to quantify the diversity of the ordering symbols (patterns) derived from a scalar time series. Note that the appropriate symbol sequence arises naturally from the time series and no model-based assumptions are needed. In fact, the necessary ‘partitions’ are devised by comparing the order of neighbouring relative values rather than by apportioning amplitudes according to different levels. This technique, as opposed to most of those in current practice, takes into account the temporal structure of the time series generated by the physical process under study. This feature allows us to uncover important details concerning the ordinal structure of the time series [41–43] and can also yield information about temporal correlation [31,32]. It is clear that this type of analysis of time series entails losing some details of the original series’ amplitude information. Nevertheless, by just referring to the series’ intrinsic structure, a meaningful difficulty reduction has indeed been achieved by BP with regard to the description of complex systems. The symbolic representation of time series by recourse to a comparison of consecutive (*τ*=1) or non-consecutive (*τ*>1) values allows for an accurate empirical reconstruction of the underlying phase-space, even in the presence of weak (observational and dynamic) noise [40]. The advantages of the method reside in (i) its simplicity, we need few parameters: the pattern length/embedding dimension *D* and the embedding delay *τ*, and (ii) the extremely fast nature of the pertinent calculation process [44]. The BP methodology can be applied not only to time series representative of low-dimensional dynamical systems, but also to any type of time series (regular, chaotic, noisy or reality based). In fact, the existence of an attractor in the *D*-dimensional phase space is not assumed. The only condition for the applicability of the BP methodology is a very weak stationarity assumption (that is, for *k*≤*D*, the probability for *x*_{t}<*x*_{t+k} should not depend on *t* [40]).

To use the Bandt & Pompe [40] methodology for evaluating the PDF, *P*, associated with the time series (dynamical system) under study, one starts by considering partitions of the pertinent *D*-dimensional space that will hopefully ‘reveal’ relevant details of the ordinal structure of a given one-dimensional time series with embedding dimension *D*>1 () and embedding time delay *τ* (). We are interested in ‘ordinal patterns’ of order (length) *D* generated by (*s*)↦(*x*_{s−(D−1)τ},*x*_{s−(D−2)τ},…,*x*_{s−τ},*x*_{s}), which assigns to each time *s* the *D*-dimensional vector of values at times *s*,*s*−*τ*,…,*s*−(*D*−1)*τ*. Clearly, the greater the *D*-value, the more information on the past is incorporated into our vectors. By ‘ordinal pattern’ related to the time (*s*), we mean the permutation *π*=(*r*_{0},*r*_{1},…,*r*_{D−1}) of [0,1,…,*D*−1] defined by *x*_{s−rD−1τ}≤*x*_{s−rD−2τ}≤⋯≤*x*_{s−r1τ}≤*x*_{s−r0τ}. In order to get a unique result, we set *r*_{i}<*r*_{i−1} if *x*_{s−ri}=*x*_{s−ri−1}. This is justified if the values of *x*_{t} have a continuous distribution so that equal values are very unusual. Thus, for all the *D*! possible permutations *π* of order *D*, their associated relative frequencies can be naturally computed by the number of times this particular order sequence is found in the time series divided by the total number of sequences. The embedding dimension *D* plays an important role in the evaluation of the appropriate probability distribution because *D* determines the number of accessible states *D*! and also conditions the minimum acceptable length *M*≫*D*! of the time series that one needs in order to work with reliable statistics [42].

### (b) Multi-scale entropy–complexity causality plane

The above discussion is based on information theory quantifiers *H*_{S} and *C*_{JS} evaluated using BP’s PDF, which allow us to define the causality information: *H*×*C*. The Shannon entropy–complexity causality plane, *H*×*C*, is based only on global characteristics of the associated time-series BP PDF as both quantities are defined in terms of Shannon entropies. In this case, *H*×*C*, the variation range is (with and the minimum and maximum statistical complexity values, respectively, for a given *H*_{S} value [39]). This causal information plane has been profitably used to separate and differentiate among chaotic and deterministic systems [41,42]; for visualization and characterization of different dynamical regimes when the system parameters vary [41,45,46]; to study time dynamic evolution [47]; to identify periodicities in natural time series [48]; to identify deterministic dynamics contaminated with noise [49,50]; and to estimate intrinsic time scales of delayed systems [51–53], among other applications (see [54] and references therein).

BP suggested working with 4≤*D*≤6 and specifically considered an embedding delay *τ*=1 in their cornerstone paper [40]. However, other values of *τ* can provide additional information since the embedding delay *τ* is the time separation between symbols, and it physically corresponds to multiples of the sampling time of the signal under analysis. More specifically, different time scales are considered by changing the embedding delays of the symbolic reconstruction [33]. The underlying chaotic or stochastic nature of a system may depend on the resolution of the data record [33]. Thus, it is more appropriate to define the concept of deterministic or stochastic behaviour on a certain range of scales. This is to say, a scale-dependent scheme should be considered when dealing with complex multi-scaled neuronal data. The main idea is therefore to generalize the estimation of both symbolic quantifiers, permutation entropy and permutation statistical complexity, accounting for different embedding delays. We refer to the multi-scale entropy–complexity causality plane as the parametric curve described by the permutation quantifiers estimated from a time series with the embedding delay *τ* as a parameter, and by considering a fixed embedding dimension *D*. The importance of selecting an appropriate embedding delay *τ* in the estimation of the permutation quantifiers (*H*_{S} and *C*_{JS}) resides in estimating the intrinsic time scales of delayed systems [31–33].

## 3. Computational modelling of anticipated/delayed synchronization in the brain

In order to investigate the synchronization properties between populations representing cortical regions, we build two populations (figure 1*a*) composed of hundreds of neurons described by the Izhikevich model [55]:
3.1and
3.2where *v* is the membrane potential and *u* is the recovery variable which accounts for activation (inactivation) of K^{+} (Na^{+}) ionic currents. *I*_{x} account for the currents provided by the interaction with other neurons and external inputs. If *v*≥30 mV, then *v* is reset to *c* and *u* to *u*+*d*. For each excitatory neuron, the dimensionless parameters are: (*a*,*b*)=(0.02,0.2) and (*c*,*d*)=(−65,8)+(15,−6)*σ*^{2}. Similarly for each inhibitory neuron: (*a*,*b*)=(0.02,0.25)+(0.08,−0.05)*σ* and (*c*,*d*)=(−65,2). Here, *σ* is a random variable uniformly distributed on the interval [0,1]. These parameters determine the spiking behaviour of each neuron mimicking known types of cortical neurons [55].

Neurons are connected to each other through chemical synapses mediated by AMPA (A) and GABA_{A} (G). The synaptic currents are given by
3.3where *x*=*A*,*G*,*E*_{A}=0 mV, *E*_{G}=−65 mV and *r*_{x} are the fraction of bound synaptic receptors whose dynamics is given by
3.4Without loss of generality, we take *τ*_{A}=5.26 ms and *τ*_{G}=5.6 ms [56]. The summation over *k* stands for presynaptic spikes at times *t*_{k}. Each neuron produces an independent spike train described by a Poisson distribution with rate *R*_{P}. The input reproduces excitatory synapses (with conductances *g*_{E}=0.5 nS) from *n* presynaptic neurons external to the population, each spiking with a Poisson rate *R*_{P}/*n*. Neurons can also receive a constant external current *I*_{c}. We use Euler’s method for numerical integration with a time step of 0.05 ms.

The sender (S) population is composed of 500 neurons (400 excitatory and 100 inhibitory), each one receiving 50 synapses (10% connectivity) from randomly selected neurons of the same population [18]. Assuming a total external rate *R*_{P}=2400 Hz and an external constant current *I*_{c}=0, the mean membrane potential *V* of this population oscillates with a mean period *T*_{S}≈130 ms (see the black curve in figure 1*b*,*c*), corresponding to a frequency *f*≈7.7 Hz, which is related to theta oscillations reported in several experiments [57–59].

The receiver (R) population is also composed of 500 neurons (80% excitatory and 20% inhibitory). Each excitatory neuron receives 40 excitatory synapses from excitatory neurons belonging to R, 10 synapses from the inhibitory neurons in the R population (with conductances *g*_{IR}) and 20 synapses from excitatory neurons from the S population (with conductances *g*_{SR}). Each inhibitory neuron in the R population receives 10 inhibitory synapses from interneurons within the same R population, 40 excitatory synapses from neurons within the R population and 20 synapses from excitatory neurons belonging to the S population. All presynaptic neurons are randomly selected. The conductances *g*_{SR} and *g*_{IR} are key parameters to control the phase difference between the activity of the sender and the receiver populations (figure 1*a*). Unless otherwise stated all other excitatory (inhibitory) synaptic conductance are fixed to the value *g*_{E}=0.5 nS (*g*_{I}=4.0 nS). We analyse long time series in which the S and the R populations oscillate with *f*≈7.7 Hz (see the time series in figure 1*b*,*c*). For this frequency, an inhibitory conductance of *g*_{IR}=8 nS produces DS, while *g*_{IR}=4 nS produces AS.

## 4. The experimental sensorimotor cortex data

Local field potential data were recorded via up to 15 bipolar microelectrodes (51 μm diameter, 2.5 mm tip separation) chronically implanted in the sensorimotor cortex (right hemisphere) of an adult male rhesus macaque monkey, as described by Brovelli *et al.* [14]. Data were acquired during a GO/NO-GO task. The monkey was required to depress a hand lever until two consecutive visual stimuli appeared on the screen. There were two stimulus types. In response to one stimulus type, the monkey was required to release the hand lever (GO trials). In response to the other stimulus type, the monkey was required to keep the hand lever depressed (NO-GO trials). Our analysis focuses on 710 trials of the 90 ms period (18 points, 200 Hz sample rate) before the visual stimulus onset (wait window). Only correct trials (both GO and NO-GO) were analysed. Considering the whole task, each trial lasted for 500 ms.

We analyse experimental LFP data from electrodes at four cortical sites: primary motor, primary somatosensory and two sites at the posterior parietal cortices (figure 1*d*,*e*). These sites are synchronized and the peak frequency in the coherence spectrum is *f*≈24 Hz. Based on previous Granger causality measures reported in [14,18], we know the direction of the information flow between these areas. Considering the two electrodes in the parietal cortex the sender leads the receiver and the relative synchronization is positive, characterizing DS. On the contrary, the primary somatosensory cortex LFP Granger causes the primary motor cortex LFP but lags behind it, yielding an AS regime.

## 5. Estimating the intrinsic time scales of delayed and anticipated cortical areas

Specific cognitive tasks require the activation of different brain regions and patterns. Therefore, neuronal population models should encompass two main aspects. First, to capture the large-scale inter-areal behaviour at multiple temporal scales as well as neuronal scale features. Second, to relate the activity patterns during different situations to the underlying anatomical connectivity of the brain. Here we use the model of two brain regions developed in [18]. The brain regions are coupled with a well-defined directional influence that displays similar features to those observed in the experimental data [14]. The model is inspired by the theoretical framework of AS developed in the field of dynamical systems [18]. That is, we take a dynamical systems model of two cortical regions, coupled with a well-defined directional influence, to generate the LFPs exhibiting DS and AS as in [18]. Depending on the synaptic conductances, the system can manifest DS, with *τ*_{Lag}>0 as in figure 1(*b*,*c*), or AS, with *τ*_{Lag}<0 as in figure 1(*d*,*e*). Figure 1 shows the average membrane potential *V* of sender (black) and receiver (grey) populations in DS (*b*,*d*) and AS (*c*,*e*) regimes. More details of the model can be found in [18].

In the following, we estimate the intrinsic time scales of the DS (AS) between cortical areas by investigating the multi-scale entropy–complexity causality plane. AS occurs when a unidirectional influence from a dynamical system (the sender) to another dynamical system (the receiver) is accompanied by a negative phase difference (or time lag) between sender and receiver [17,18]. Therefore, when considering the AS regime, the receiver’s trajectory can predict the sender’s future behaviour. Figure 2 shows the permutation complexity *C*_{JS} versus the normalized permutation entropy *H*_{S} considering different embedding delay *τ*, with a fixed embedding dimension *D*=6. The dark grey circles and triangles in figure 2*a*,*b* correspond to the sender, for the AS and the DS regimes, respectively. The light grey circles and triangles in figure 2*c*,*d* correspond to the receiver, for the AS and the DS regimes, respectively. *C*_{JS} is maximized when the embedding delay *τ* of the symbolic reconstruction matches the intrinsic time delay *τ* of the system. In the case of AS the maxima are at *τ*_{Sender}=3.4 ms for the sender, and *τ*_{Receiver}=0.8 ms for the receiver. Thus, the ‘relative synchronization time’ between the receiver and the sender is *τ*_{Receiver}−*τ*_{Sender}=−2.6 ms (AS). By contrast, for the case of DS the maxima are at *τ*_{Sender}=0.05 ms for the sender, and *τ*_{Receiver}=0.55 ms for the receiver. That is, the time delay between them is *τ*_{Receiver}−*τ*_{Sender}=0.5 ms (DS).

In the electronic supplementary material, we present the figures corresponding to the variation of normalized permutation Shannon entropy *H*_{S} and the MPR permutation statistical complexity *C*_{JS} as a function of the time lag *τ* in the case of sender and receiver, for the computational model proposed by Matias *et al.* [18], as well as for the experimental data of Brovelli and co-workers [14].

We have analysed the above artificial LFP data obtained from the model of two brain regions developed by Matias *et al.* [18], coupled with a well-defined directional influence. It has been empirically observed that a dominant directional influence between areas of sensorimotor cortex may be accompanied by either a negative or a positive time delay [14]. Indeed Brovelli *et al.* [14] reported that, in monkeys engaged in processing a cognitive task, a dominant directional influence from one area of sensorimotor cortex to another may be accompanied by either a negative or a positive time delay in analogy to the dynamical systems modelling proposed in [18]. In the following, we estimate the intrinsic time scale of the LFP data from the cortical dataset of Brovelli and co-workers [14]. We consider therefore two cases: one with negative and the other with a positive time delay estimated by Brovelli *et al.* [14]. We compute the ‘relative synchronization phase’ between the different areas by means of the multi-scale entropy–complexity, the relative ‘synchronization phase’, between electrode sites. Figure 3 shows the permutation complexity *C*_{JS} versus the normalized permutation Shannon entropy *H*_{S} considering different embedding delay *τ* adjusted to the experimental data, for a fixed embedding dimension *D*=6. The dark grey circles and triangles in figure 3*a*,*b* correspond to the sender, for AS and DS regimes, respectively, while the light grey circles and triangles in figure 3*c*,*d* correspond to the receiver, for AS and DS regimes, respectively. In the case of a negative phase shift between cortical areas the maxima are at *τ*_{Sender}=775 ms and *τ*_{Receiver}=755 ms. Thus, the ‘relative synchronization time’ between the receiver and the sender is *τ*_{Receiver}−*τ*_{Sender}=−20 ms (AS). In contrast, for the case of a positive relative phase between areas the maxima are at *τ*_{Sender}=7.80 ms and *τ*_{Receiver}=8.25 ms. That is, the time delay between both is *τ*_{Receiver}−*τ*_{Sender}=0.45 ms (DS). See the electronic supplementary material for details.

Each of the synthetic datasets has *M*=241 599 points, which guarantees no bias deviation in the estimation of the information quantifiers. Regarding the experimental data we have considered *M*=12 780 points, which corresponds to 710 trials and which might produce some bias deviation from the methodology. Despite this possible limitation, the results we obtained in this section are in quite good agreement with the typical values obtained by other authors for the ‘relative synchronization time’ of AS/DS synchronization [18].

Confidence error intervals cannot be provided within the BP methodology. As we mentioned above, the selection of the embedding dimension, *D*, is relevant for obtaining an appropriate probability distribution because not only does *D* determine the number of accessible states (equal to *D*!) but also the length of the time series, *M*, needs to have a reliable statistics and therefore the requirement is that the condition *M*≫*D*! must be satisfied. BP suggested working with 4≤*D*≤6 in their cornerstone paper [40]. We have currently analysed also the cases of *D*=4 and *D*=5 without finding any sensible difference from the case of *D*=6 when analysing the different values of *τ*.

Thus, the location of the maximum in the multi-scale *H*×*C* causality plane allows us to infer useful information about the underlying dynamics of the LFP’s time series. It provides us with a novel methodology for estimating the ‘relative synchronization time’ between two areas of the monkey brain. In a sender–receiver configuration, the direction of information flow is from the sender to the receiver. Hence, the ‘relative synchronization phase’ is detected by using the permutation entropy quantifier analysing the sender and receiver signals, when the embedding delay matches the intrinsic time delay *τ* of the system within the multi-scale entropy–complexity causality plane.

## 6. Discussion and conclusion

Cortical neurons are mainly coupled via chemical synapses, which can be excitatory or inhibitory. In both cases, the coupling is directional and highly nonlinear, typically requiring a supra-threshold activation (e.g. a spike) of the presynaptic neuron to trigger the release of neurotransmitters. These neurotransmitters need to diffuse through the synaptic cleft and bind to the receptors in the membrane of the postsynaptic neuron. Binding leads to the opening of specific channels, allowing ionic currents to change the postsynaptic membrane potential [60]. This means not only that the membrane potentials are not directly coupled, but also that the synapses themselves are dynamical systems. Matias *et al.* [18,28] proposed to bridge this gap by investigating whether AS can occur in biophysically plausible model neurons coupled via chemical synapses, when replacing the typical self-feedback loop by a dynamical inhibitory loop mediated by an interneuron. Such an inhibitory feedback loop is one of the most canonical neuronal motifs in the brain [61,62]. AS might therefore have important roles in information transmission in the brain: the predictive powers of the brain could emerge from different dynamics of the rhythms generated by the neurons.

In this paper, we used numerical and experimental data to estimate the time delay between the activities of two brain areas using the permutation quantifiers in the multi-scale entropy–complexity causality plane. We showed that the AS/DS dynamics between different areas can be reliably characterized by applying the permutation entropy and permutation statistical complexity as a function of the embedding delay to the LFPs. The scale is explicitly incorporated in this approach by changing the embedding delay *τ*. The location of the ‘relative synchronization time’, using the multi-scale entropy–complexity plane, allows us to infer useful information about the underlying dynamics of the complex LFP time series, and to characterize the system’s dynamics as AS/DS. The obtained numerical and experimental results confirm that this multi-scale symbolic information-theory approach provides a conceptually simple and computationally efficient tool for characterizing complex time series of brain circuits. The existence of AS mediated by a dynamical inhibition unveils several possibilities in the investigation of synchronized activity in the brain. Controlling delay-induced instabilities through AS synchronization may also be a learning mechanism used by the brain after processing the information that is available in the system. Further experimental verification of anticipating chaotic synchronization can be of ultimate help in future research of brain dynamics, and to understand further how information is processed by the brain.

Although several brain regions show significant specialization, higher functions such as cross-modal information, integration, abstract reasoning and conscious awareness are viewed as emerging from interactions across distributed functional networks. Indeed, most brain functions are thought to rely on the interrelationship between segregation and integration. The coexistence of these two principles is considered to be the origin of neural complexity [63]. Neural connectivity is a way by which neurons could generate diverse patterns of response and mutual statistical dependence. Synaptic connectivity allows neurons to act both independently and collectively. Thus, the brain function is fundamentally integrative; it requires that components and elementary processes work together, giving rise to complex patterns. Connectivity is essential for integrating the actions of individual neurons and therefore for enabling cognitive processes, such as perception, attention and memory. Connectivity translates unitary events at the cellular scale into large-scale patterns produced by neuronal ensembles. As non-causal mutual information fails to distinguish information that is actually exchanged from shared information due to common history and input signals [64], the current approach based on the multi-scale entropy–complexity can be very powerful to investigate processing between brain areas. This is important not only from a theoretical point of view—it might also help to determine which areas of the cortex could have a higher level of information, and to evaluate how causal interactions in neural dynamics would be modulated by behaviour. We believe that this will become an important tool for future research on the encoding capacity of biologically realistic neural networks.

## Authors' contributions

All authors contributed equally to the paper.

## Competing interests

We declare we have no competing interests.

## Funding

Research supported by PIP 0255/11 CONICET, Argentina (F.M.). Research project FIS2012-30634 INTENSE@COSYP MINECO (Spain) and Feder (C.R.M.). C.R.M. and F.S.M. thank the Brazilian project Cincia Sem Fronteiras PVE 88881.068077/2014-01 (CAPES, Brazil).

## Acknowledgements

O.A.R. and F.M. are members of the National Research Career of CONICET Argentina and acknowledge support by CONICET, Argentina.

## Footnotes

One contribution of 13 to a theme issue ‘Topics on non-equilibrium statistical mechanics and nonlinear physics (II)’.

- Accepted July 28, 2015.

- © 2015 The Author(s)