## Abstract

Neurons tend to fire a spike when they are near a bifurcation from the resting state to spiking activity. It is a delicate balance between noise, dynamic currents and initial condition that determines the phase diagram of neural activity. Many possible ionic mechanisms can be accounted for as the source of spike generation. Moreover, the biophysics and the dynamics behind it can usually be described through a phase diagram that involves membrane voltage versus the activation variable of the ionic channel. In this paper, we present a novel methodology to characterize the dynamics of this system, which takes into account the fine temporal ‘structures’ of the complex neuronal signals. This allows us to accurately distinguish the most fundamental properties of neurophysiological neurons that were previously described by Izhikevich considering the phase-space trajectory, using a time causal space: statistical complexity versus Fisher information versus Shannon entropy.

## 1. Introduction

From a theoretical point of view, the brain constitutes a complex and dynamic system whose state variables represent information both external (e.g. perception of certain stimuli) and internal (e.g. reverberation of memories). Neurons are the main processing and transmission units in the brain and thus are critical to cognition. Therefore, in order to detect subtle changes in brain activity, we have to investigate the intrinsic dynamics of the neurons. This requires the study of the spontaneous activity potential of cell membranes, known as action potentials or ‘spikes’, which allows us to describe in a fairly comprehensive way neural behaviour.

The modelling of neurons and neural circuits on the basis of cellular, synaptic biophysics and spiking models includes the dynamic changes of the neuron’s membrane potential [1]. The simplest models representing a minimal biophysical interpretation for an excitable neuron are conductance-based models. The first model of spiking neurons was proposed by Alan Lloyd Hodgkin and Andrew Huxley in 1952 [2]. It described the ionic mechanisms underlying the initiation and propagation of the action potentials in the squid giant axon. Hodgkin and Huxley (HH) explained the ionic mechanisms underlying the initiation and propagation of action potentials, through a set of nonlinear ordinary differential equations that approximate the electrical characteristics of the excitable cell [2]. In its simplest version, the HH model represents a neuron by a single isopotential electrical compartment, neglects ion movements between subcellular compartments, and represents only ion movements between the inside and outside of the cell. Applying nonlinear dynamic theory, Izhikevich proposed a classification of neurons depending on bifurcation and resting state. While there are a huge number of possible ionic mechanisms of excitability and spike generation, there are just four bifurcation mechanisms that can result in a transition from resting state to spiking. These bifurcations divide neurons into four categories: integrators or resonators, monostable or bistable [3–5].

Izhikevich proposed a simpler model in which the number of variables is considerably smaller in comparison to the HH model. Then the behaviour of a neuron can be described by the ‘simple spiking model’ [4], which can reproduce 20 of the most fundamental neurocomputational features. Although not all the parameters of the model are feasible to be measured directly in the laboratory, as they have no direct biological interpretation, the model reproduces the essential characteristics observed in neurons. In contrast to the HH model, the Izhikevich simple model is much more efficient from a computational point of view, which proves useful in large-scale computations [6–9]. Furthermore, in the framework of this model, several analytical methods have been developed in order to study these relevant features [6]. However, quantification of the typical oscillatory activity patterns using an information theoretical approach accounting for the causality of the signals is still missing.

In particular, many important aspects of brain activity, such as rhythms, are of functional importance to understand how information is processed in the mammalian brain. Neural oscillatory activity patterns are rhythmic neural activities in the brain and can be generated either by mechanisms within individual neurons or by interactions between neurons. At the level of individual neurons, patterns of oscillatory activity can appear either as oscillations in the membrane potential or as rhythmic patterns of action potentials, which then produce the oscillatory activation of post-synaptic neurons. The relationship between these devices and behaviour, therefore, could provide a novel understanding about the functions of brain oscillations [10–14]. We use the Bandt–Pompe permutation methodology for the evaluation of the probability distribution function (PDF) associated with a time series [15]. Based on the quantification of the ordinal ‘structures’ present in the neuron’s membrane potential and their local influence on the associated probability distribution, we incorporate the time series’ own temporal causality through an algorithm that is easy to implement and compute. More specifically, in this paper we propose a versatile method to quantify the 20 most fundamental neurocomputational features of oscillatory patterns of biological neurons, considering subtle measures accounting for the causal information: Shannon permutation entropy [16,17], Fisher permutation information [18,19] and Martín–Plastino–Rosso (MPR) permutation statistical complexity [16,17]. Our approach allows us to estimate the ‘clustering properties’ of these different neuronal structures, to quantify the causality of the signal, and to infer the emergent dynamical properties of the system through two- and three-dimensional representations.

## 2. Simple model of spiking neurons

Bifurcation methodologies [5] allow us to accurately reproduce the biophysical properties of HH neuronal models just by taking a two-dimensional system of ordinary differential equations and four different parameters. This model is named the *simple model of spiking neurons* [4]:
2.1and
2.2with the auxiliary after-spike resetting
2.3Here *v* is the membrane potential of the neuron, and *u* is a membrane recovery variable, which accounts for the activation of K^{+} ionic currents and inactivation of Na^{+} ionic currents, and gives negative feedback to *v*. Thus, we are just considering a system of ordinary differential equations of two variables *u* and *v*, and all the known types of neurons can be reproduced by taking different values of the four parameters *a*, *b*, *c* and *d*. After the spike reaches its apex at +30 mV (not to be confused with the firing threshold), the membrane voltage and the recovery variable are reset according to equation (2.3). The variable *I* accounts for the inputs to the neurons [3].

In figure 1, we use the simple model of spiking neurons to reproduce the 20 most fundamental neurocomputational properties of biological neurons by changing the different parameters *a*, *b*, *c*, *d* and *I* [4]. Figure 1 shows different neuronal behaviours after taking the current inputs: (A) tonic spiking, (B) phasic spiking, (C) tonic bursting, (D) phasic bursting, (E) mixed model (bursting then spiking), (F) spike frequency adaptation, (G) class 1 excitability, (H) class 2 excitability, (I) spike latency, (J) subthreshold oscillations, (K) frequency preference and resonance, (L) integration and coincidence detection, (M) rebound spike, (N) rebound burst, (O) threshold variability, (P) bistability of resting and spiking states, (Q) depolarizing after-potentials, (R) accommodation, (S) inhibition-induced spiking and (T) inhibition-induced bursting. The corresponding parameter values for the different neuronal behaviours denoted by (A)–(T) in figure 1 are given in the electronic supplementary material, table S1.

In this paper, we use this biologically plausible and computationally efficient model to generate a simulation of oscillatory activity patterns of neurons. In the following section, we summarize the basic concepts of causal information quantifiers.

## 3. Information theory quantifiers

### (a) Shannon entropy, fisher information measure and Martín–Plastino–Rosso statistical complexity

Sequences of measurements (or observations) constitute the basic elements for the study of natural phenomena. In particular, from these sequences, commonly called time series, one should judiciously extract information on the dynamical systems under study. We can define an information theory quantifier as a measure that is able to characterize some property of the PDF associated with these time series of a given raw signal (i.e. a neuron’s membrane potential).

Entropy, regarded as a measure of uncertainty, is the most paradigmatic example of these quantifiers. Given a continuous PDF *f*(*x*) with and , its associated *Shannon entropy* *S* [20] is
3.1Given a time series , a set of *M* measures of the observable and the associated PDF, given by *P*≡{*p*_{j};*j*=1,…,*N*} with and *N* the number of possible states of the system under study, the Shannon logarithmic information measure [20] is defined by
3.2This 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. Such is not the case with the *Fisher information measure* [21,22]
3.3which constitutes a measure of the gradient content of the distribution *f* (continuous PDF), thus being quite sensitive even to tiny localized perturbations.

The Fisher information measure can be variously interpreted as a measure of the ability to estimate a parameter, as the amount of information that can be extracted from a set of measurements, and also as a measure of the state of disorder of a system or phenomenon [22,23], its most important property being the so-called Cramer–Rao bound. It is important to remark that the gradient operator significantly influences the contribution of minute local *f* variations to the Fisher information value, so that the quantifier is called a ‘local’ one. Note that Shannon entropy decreases with a skewed distribution, while Fisher information increases in such a case. Local sensitivity is useful in scenarios whose description necessitates an appeal to a notion of ‘order’ [24–26]. The concomitant problem of loss of information due to discretization has been thoroughly studied (see, for instance, [27–29] and references therein) and, in particular, it entails the loss of Fisher’s shift invariance, which is of no importance for our present purposes.

For Fisher information measure computation (discrete PDF), we follow the proposal of Dehesa and co-workers [30] based on the amplitude of probability *f*(*x*)=*ψ*(*x*)^{2}:
3.4Its discrete normalized version (0≤*F*≤1) is now
3.5Here the normalization constant *F*_{0} reads
3.6

If our system lies in a very ordered state, we can take it to be described by a PDF given by *P*_{0}={*p*_{k}≅1;*p*_{i}≅0 ∀ *i*≠*k*;*i*=1,…,*N*} (with *N* the number of states of the system). In consequence we have a Shannon entropy *S*[*P*_{0}]≅0 and a normalized Fisher information measure . On the other hand, when the system under study is represented by a very disordered state, one can take this particular state to be described by a PDF given by the uniform distribution *P*_{e}={*p*_{i}=1/*N*, ∀ *i*=1,…,*N*}. We obtain while *F*[*P*_{e}]≅0. One can state that the general behaviour of the Fisher information measure is opposite to that of the Shannon entropy [31].

It is well known, however, that ordinal structures present in a process are not quantified by randomness measures and, consequently, measures of statistical or structural complexity are necessary for a better understanding (or characterization) of the system dynamics represented by their time series [32]. The opposite extremes of perfect order (a periodic sequence, a regular crystal, for example) and maximal randomness (i.e. a fair coin toss, an ideal gas) are very simple to describe because they do not have any structure. The complexity should be zero in these cases, that is *C*[*P*_{0}]=*C*[*P*_{e}]=0. At a given distance from these extremes, a wide range of possible degrees of physical structure exists. The complexity measure allows one to quantify this array of behaviour in between ‘perfect order’ (e.g. regular crystal) and ‘complete disorder’ (e.g. ideal gas). Complexity can be characterized by a certain degree of organization, structure, memory, regularity, symmetry and patterns [33]. The complexity measure does much more than satisfy the boundary conditions of vanishing in the high- and low-entropy limits. In particular, maximum complexity occurs in the region between the system’s perfectly ordered state and the perfectly disordered one. Complexity measures allow us to detect essential details of the dynamics, to discriminate different degrees of periodicity associated with the period doubling bifurcation route to chaos, and more importantly to characterize the correlational structure of the orderings present in the time series.

We consider the *MPR statistical complexity* [34] 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.* [35], this statistical complexity measure (SCM) is defined through the product
3.7of the normalized Shannon entropy
3.8with (0≤*H*≤1), and the disequilibrium defined in terms of the Jensen–Shannon divergence
3.9with
3.10The 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 one 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 [36]. Note that the above-introduced SCM depends on two different probability distributions, one associated with the system under analysis, *P*, and the other with uniform distribution, *P*_{e}. Furthermore, it was shown that, for a given value of *H*, the range of possible *C*_{JS} values varies between a minimum and a maximum , restricting the possible values of the SCM in a given complexity–entropy plane [37]. 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 SCM. In order to calculate the three information theory-derived quantifiers mentioned previously, a probability distribution should be estimated from the time series of the system.

The study and characterization of time series by recourse to information theory tools assume that the underlying PDF is given *a priori*. By contrast, part of the concomitant analysis involves extracting the PDF from the data and there is no univocal procedure with which everyone agrees. Almost 10 years ago, Bandt & Pompe (BP) introduced a successful methodology for the evaluation of the PDF associated with scalar time-series data using a symbolization technique [15]. For a didactic description of the approach, as well as its main biomedical and econophysics applications, see [38].

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 [26,39,40] and can also yield information about temporal correlation [16,17]. It is clear that this type of analysis of time series entails losing some details of the original series’ amplitude information. Nevertheless, just by referring to the series’ intrinsic structure, a meaningful difficulty reduction has indeed been achieved by Bandt & Pompe 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 [15]. Furthermore, the ordinal patterns associated with the PDF are invariant with respect to nonlinear monotonic transformations. Accordingly, nonlinear drifts or scaling artificially introduced by a measurement device will not modify the estimation of quantifiers, a nice property if one deals with experimental data (e.g. [41]). These advantages make the BP methodology more convenient than conventional methods based on range partitioning (i.e. PDF based on histograms).

Additional advantages of the method reside in (i) its simplicity, as 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 [42]. 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 stationary assumption (that is, for *k*≤*D*, the probability for *x*_{t}<*x*_{t+k} should not depend on *t* [15]).

To use the Bandt & Pompe [15] 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:
3.11In the last expression, the symbol ♯ stands for ‘number’. Thus, an ordinal pattern probability distribution *P*={*p*(*π*_{i}),*i*=1,…,*D*!} is obtained from the time series.

Consequently, it is possible to quantify the diversity of the ordering symbols (patterns of length *D*) derived from a scalar time series, by evaluating the so-called permutation Shannon entropy, Fisher permutation information measure and permutation MPR statistical complexity. Of course, 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 [39].

Regarding the selection of the parameters, Bandt & Pompe suggested working with 4≤*D*≤6 and specifically considered an embedding delay *τ*=1 in their cornerstone paper [15]. Nevertheless, it is clear that other values of *τ* could provide additional information. It has been shown recently that this parameter is strongly related, if it is relevant, to the intrinsic time scales of the system under analysis [43–45].

The local sensitivity of Fisher information measure for discrete PDFs is reflected in the fact that the specific ‘*i* ordering’ of the discrete values *p*_{i} must be seriously taken into account in evaluating the sum in equation (3.5) [18,19]. The pertinent numerator can be regarded as a kind of ‘distance’ between two contiguous probabilities. Thus, a different ordering of the pertinent summands would lead to a different Fisher information value. In fact, if we have a discrete PDF given by *P*={*p*_{i},*i*=1,…,*N*} we will have *N*! possibilities for the *i* ordering.

The question is, which is the arrangement that one could regard as the ‘proper’ ordering? The answer is straightforward in some cases, a histogram-based PDF constituting a conspicuous example. For such a procedure, one first divides the interval [*a*,*b*] (with *a* and *b* the minimum and maximum amplitude values in the time series, respectively) into a finite number of non-overlapping sub-intervals (bins). Thus, the division procedure of the interval [*a*,*b*] provides the natural order sequence for the evaluation of the PDF gradient involved in the Fisher information measure. In our current paper, we chose for the BP PDF the lexicographic ordering given by the algorithm of Lehmer (http://www.keithschwarz.com/interesting/code/factoradic-permutation/FactoradicPermutation), amongst other possibilities, because it provides a better distinction of different dynamics in the Fisher versus Shannon plane [25,26].

## 4. An information theoretic characterization of the fundamental neuro- computational features of biological neurons

Most neurons are quiescent but can fire spikes when properly stimulated, and they typically respond by producing complex spike sequences. This shows the intrinsic dynamics of the neuron and to some extent the temporal characteristics of the stimulus. Understanding the features of neuronal responses encoding the variations in the stimuli is an important challenge in neuroscience. Oscillations in the brain may be generated by non-invasive brain stimulation, either by intrinsic mechanisms or by interactions between them [10–14], and may have a crucial role in feature binding, information transmission and the generation of rhythmic motor output even in isolation from motor and sensory feedback [46]. They are indeed manifested as a variety of rhythms that differ in their frequency, origin and reactivity to changes in sensory input and external stimuli [10–14]. A periodic stimulation can be used for interventions into the timing of biological rhythms to reveal their causal implication in behaviour [14]. The significance of rhythmic patterns for information processing has been pointed out by Kayser *et al.* [47], who presented naturalistic auditory stimuli to monkeys and used mutual information to quantify the amount of information about the identity of each stimulus encoded in either phase or amplitude per oscillatory frequency. Interestingly, the oscillatory phase of local field potential (LFP) oscillations contained significantly more information than the amplitude, with highest information at low frequencies (48 Hz). Moreover, oscillations can modulate information processing, as rhythmic inhibition plays an important role in oscillations throughout the brain, eliciting rebound bursts and resulting in re-excitation of neurons, regulating the neuronal firing in the brain [48,49].

Indeed, in the brain, a variety of rhythms have been described that differ in their frequency, origin and reactivity to changes in sensory input and task demands [10–14]. Brain oscillations are characterized by rhythmic changes in LFPs. Feedback loops across neurons contribute to the synchronization of cortical activities and modulation of oscillatory phase relationships among neuronal populations, providing deeper insights into how information is processed in the brain [14]. At the level of individual neurons, rhythmic patterns are produced by oscillations in membrane potential or as rhythmic patterns of action potentials, which cause oscillatory activation of post-synaptic neurons. As oscillatory rhythmic patterns can carry relevant information about external stimuli, detecting dynamic changes in neural systems is one of the most important tasks in theoretical neuroscience. Here, our proposal is to quantify a variety of oscillatory activity patterns at single neuron level that are generated using the ‘simple spiking model’ of Izhikevich [4]. We use a versatile method to quantify the 20 most fundamental neurocomputational features of biological neurons, by means of an information theoretical approach. More specifically, we consider measures accounting for the causal structure of the signal of a neuron’s membrane potential: the Shannon permutation entropy, Fisher permutation information and MPR permutation statistical complexity. As we mentioned in §2, these features of biological neurons are shown in figure 1. They were obtained with the simple model of Izhikevich [4]; the simulation time used to produce figure 1 is 100 ms, with a resolution of 0.25 ms. In order to perform analyses within the BP formalism, we need to have a much larger number of points in the simulation of membrane potential responses (*M*≫*D*!). We took 1800 trials (repetitions) to obtain 180 000 simulation points for each case. We used the Bandt & Pompe [15] methodology for evaluating the PDF, *P*, associated with the time series, considering an embedded dimension *D*=6 and time lag *τ*=1. This embedding dimension (pattern length) is enough to efficiently capture the information causality of the ordinal structure of the time series [15].

Figure 2*a* shows the informational causal plane of entropy versus complexity, *H*×*C*. Note that the MPR statistical complexity grows linearly as the Shannon entropy becomes higher when considering the different kinds of dynamics. The continuous lines represent the curves of maximum and minimum statistical complexity, and , respectively, as functions of the normalized Shannon entropy [37]. The degree of order decreases as entropy increases, and thus a system with a lower degree of entropy is characterized by a higher degree of order. Figure 2*b* shows Fisher permutation information, equation (3.5), versus the MPR permutation statistical complexity equation (3.9), i.e. the causal information plane *F*×*C*, for the same neurocomputational signals used above. Note that figure 2*b* presents a better distinction between the 20 neuronal rhythms. That is, permutation Fisher information behaves nonlinearly as a function of the MPR permutation statistical complexity. A similar behaviour can be observed in figure 2*c*, which shows Fisher permutation information, equation (3.5), versus permutation Shannon entropy, equation (3.8), i.e. the causal plane *F*×*H*. This is a reasonable result, as we observed that MPR statistical complexity grows linearly as the Shannon entropy becomes higher.

Applying nonlinear dynamic theory, Izhikevich proposed a classification of neurons depending on bifurcation and resting state. Izhikevich presented a simple model in which the number of variables is considerably smaller in comparison to the HH model. The behaviour of a neuron is then described by this ‘simple spiking model’. In order to do so, we take different values suggested for the four parameters *a*, *b*, *c* and *d*, and different shapes proposed for the current *I* that accounts for the inputs to the neurons. By doing so, we can reproduce 20 of the most fundamental neurocomputational features of a neuron. If we were to take the same current in all 20 cases we would not reproduce the rhythmic activities presented, which is the objective of this work, in order to characterize them using an information theoretical approach.

Thus, we follow Izhikevich’s recipe to generate the 20 different activities, which includes also considering different input currents (stimuli). Taking *τ*=1 is the main reason why some of the rhythmic activities have similar values in the information plane *H*×*C*. However, notice that we have three information planes in figure 2: *C*×*H*, *F*×*C* and *F*×*H*. Figure 3 shows the emergent dynamical properties of the system causality information through a three-dimensional representation: MPR statistical complexity versus Fisher information versus Shannon entropy. Remarkably, our approach allows us to classify the most fundamental neural dynamics by means of an information theoretical approach. That is, we quantify the causality of the signal, and infer the emergent properties of each oscillatory activity pattern within a three-dimensional informational representation by assigning a ‘cluster’ to each of these features. Hence, the results shown in figure 3 are the key to accurately quantify and distinguish the neurocomputational properties, previously described by Izhikevich considering the phase-space trajectory [4], by using instead a versatile method based on the causal space representation. The current approach could be a useful tool for online classification of neuronal oscillations when analysing the neuron’s membrane potential, and monitoring brain activity.

## 5. Discussions and conclusion

In 1929, Hans Berger observed rhythmic variations in the human electroencephalogram (EEG), and, more than 50 years later, intrinsic oscillatory behaviour has been found in mammalian neurons [46]. Periodic stimulations can induce the timing of neuronal rhythms and their causal implication in behaviour, and are characterized by rhythmic and periodic changes in LFPs. In this paper, we consider 20 of the most fundamental neurocomputational features of biological neurons, taking 1800 trials in order to obtain 180 000 simulation points. Thus, we generate repetitive neural activity or oscillatory activity patterns using each of the most prominent characteristics of neurons, and we quantify their relevant features using a versatile method by means of an information theoretical approach. More specifically, we apply a robust approach to time-series analysis on the basis of counting ordinal patterns, in the membrane potential of a neuron, by introducing the concept of permutation quantifiers to quantify the complexity of the system behind the series. This could prove useful for online classification of neuronal activity patterns when analysing oscillations in the LFP signals that are also known as brain waves, and might help us to understand the intrinsic differences between healthy and unhealthy tissues.

A characteristic of oscillatory patterns is their periodicity, that is, the organization of temporal dynamics into cycles. In this paper, we have considered different neuronal oscillation patterns with a phase-locking rhythm. This is also a common scenario when taking into account a small group of neurons, where a stable phase difference of membrane voltage oscillations can be observed [10,14,47,50]. When considering many neurons, the coincidence of electrical activities would lead to coordinated changes that would be detected as variations of the LFP and that might be reflected also on the EEG. If individual neurons fire action potentials periodically, and these events are synchronized among many cells, one would then observe periodic LFP oscillations. Hence, global changes of electrical activity associated with EEG rhythms would usually take place if synchronization occurs among distinct brain areas. All these models rely on the hypothesis that the oscillatory cycle establishes recurrent temporal patterns. This allows the brain to generate a coding of temporal relations between groups of neural elements and between neural elements and the environment. Importantly, oscillations in the brain are non-stationary and are indeed subject to dynamic changes due to phase resetting. Oscillations become potentially powerful neurocomputational tools when their phase, and also amplitude and frequency, change due to dynamic modifications in their generating system or in the input [10,14,47,50]. Thus, a much more realistic neurocomputational model should take phase resetting into account. Hence, it would be interesting to apply the current approach to models in which it is possible to generate mechanisms of oscillatory phase resetting due to external stimuli, affecting neuronal behaviour.

Bandt & Pompe [15] suggested working with 4≤*D*≤6 and specifically considered an embedding delay *τ*=1 in their cornerstone paper. 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 can be considered by changing the embedding delays of the symbolic reconstruction. The underlying chaotic or stochastic nature of a system may depend on the resolution of the data record. Thus, it is more appropriate to define the concept of deterministic or stochastic behaviour on a certain range of scales. That is to say, a scale-dependent scheme should be considered when dealing with complex multiscale neuronal data. The main idea is therefore to generalize the estimation of the symbolic quantifiers, permutation entropy and permutation statistical complexity, accounting for different embedding delays. We use the term multiscale entropy–complexity causality plane to refer to the parametric curve described by the permutation quantifiers estimated from a time series with the embedding delay *τ* as a parameter [45], and considering a fixed embedding dimension *D*. The importance of selecting an appropriate embedding delay *τ* in the estimation of the permutation quantifiers (*H* and *C*) resides in estimating the intrinsic time scales of the system. The complexity is maximized when the embedding delay *τ* of the symbolic reconstruction matches the intrinsic time scale *τ*_{0} of the system, and therefore important additional information related to the correlational structure between the components of the physical system is provided by evaluating the SCM. Figure 2*a* shows that the statistical complexity grows almost linearly as the normalized Shannon entropy value increases, when fixing *τ*=1. Preliminary results show that, by taking values of *τ* that maximize the complexity measure, the statistical complexity grows nonlinearly as the Shannon entropy becomes higher. Moreover, the 20 different neuronal rhythmic activities are more much distinct when considering the informational plane *H*×*C*. This problem will be considered in detail in a future publication.

Cracking the neural code involves finding synchronization patterns and meaning in the noisy activity of cell ensembles. Despite the fact that neuronal ensembles may show oscillatory cycles allowing one to establish a coding of temporal relationships between groups of neural elements and between neural elements and the environment, this system is far away from being stationary, as it is subject to dynamic changes of external stimuli and thus phase resetting. A proper quantification of the dynamics of the neuronal activity is important for understanding the essential mechanisms of brain encoding, and to gain knowledge of how information is processed through the brain. As non-causal mutual information fails to distinguish information that is actually exchanged from shared information due to common history and input signals [49], the current approach based on the permutation statistical complexity versus Shannon entropy versus Fisher information three-dimensional representation can be a powerful tool to investigate information processing across brain areas. This versatile causal information tool can be used to provide an appropriate quantification of the neuronal oscillation patterns that could have an important role in shaping theories of perception, cognition and neural computation. This is not only important from a theoretical point of view, but also it might 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 this paper.

## Competing interests

We declare we have no competing interests.

## Funding

This research was supported by PIP (Proyecto de Investigación Plurianual) 0255/11 CONICET, Argentina (F.M.).

## Acknowledgements

F.M. and O.A.R. are members of the National Research Career of CONICET Argentina. F.M. and O.A.R. 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 21, 2015.

- © 2015 The Author(s)