Thermodynamic Model of Criticality in the Cortex Based On EEG/ECOG Data

Criticality in the cortex emerges from the seemingly random interaction of microscopic components and produces higher cognitive functions at mesoscopic and macroscopic scales. Random graphs and percolation theory provide natural means to de- scribe critical regions in the behavior of the cortex and they are proposed here as novel mathematical tools helping us deciphering the language of the brain.


Introduction
Von Neumann emphasized the need for new mathematical disciplines in order to understand and interpret the language of the brain [53]. He has indicated the potential directions of such new mathematics through his work on self-reproducing cellular automata and fi-nite mathematics. Due to his early death in 1957, he was unable to participate in the development of relevant new theories, including morphogenesis pioneered by Turing [60] and summarized by Katchalsky [39]. Prigogine [56] developed this approach by modeling the emergence of structure in open chemical systems operating far from thermodynamic equilibrium. He called the patterns dissipative structures, because they emerged by the interactions of particles feeding on energy with local reduction in entropy. Haken [34] developed the field of synergetics, which is the study in lasers of pattern formation by particles, whose interactions create order parameters by which the particles govern themselves in circular causality.
Principles of self-organization and metastability have been introduced to model cognition and brain dynamics [36], [40], [35]. Recently the concept of self-organized criticality (SOC) has captured the attention of neuroscientists [4]. There is ample of empirical evidence of cortex conforming to the self-stabilized, scale-free dynamics of the sand pile during the existence of quasi-stable states [1,3,55]. However, the model cannot produce the emergence of orderly patterns within a domain of criticality [44,32]. Bonachela [10] describe SOC as pseudo-critical and suggest to complement self-organization with more elaborate, adaptive approaches.
Extending on previous studies, we propose to treat cortices as dissipative thermodynamic systems that by homeostasis hold themselves near a critical level of activity that is far from equilibrium but steady state, a pseudo-equilibrium. We utilize Freemans neuroscience insights manifested in the hierarchical brain model: the K (Katchalsky) sets [39,18]. Originally, K-sets have been described mathematically using ordinary differential equations (ODE) with distributed parameters and with stochastic differential equations (SDE) using additive and multiplicative noise [11,13]. This approach has produced significant results, but certain shortcoming have been pointed out as well [19,20,22]. Calculating stable solutions for large matrices of nonlinear ODE and SDE that closely correspond to chaotic ECoG activity are prohibitively difficult and time-consuming to model on both digital and analog platforms [57]. In addition, there are unsolved theoretical issues in constructing solid mathematics with which to bridge the difference across spatial and temporal scales between microscopic properties of single neurons and macroscopic properties of vast populations of neurons [21,23].
In the past decade, neuropercolation approach has proved to be an efficient tool to address the above shortcomings by implementing K-sets using concepts of discrete mathematics and random graph theory [41]. Neuropercolation is a thermodynamics-based random cellular neural network model, which is closely related to cellular automata (CA), the field pioneered by Von Neumann who anticipated the significance of CA in the context of brain-like computing [54]. The present study is based on applying neuropercolation to systematically implement Freeman's principles of neurodynamics. Brain dynamics is viewed as a sequence of intermittent phase transitions in an open system with synchronization-desynchronization effects demonstrating symmetry breaking demarcated by spatio-temporal singularities [24,48,26,29].
This work starts with the description of the basic building blocks of neurodynamics. Next, we develop a hierarchy on neuropercolation models with increasing complexity in structure and dynamics. Finally, we employ neuropercolation to describe critical behavior of brains and to interpret experimentally observed ECoG/EEG dynamics manifesting learning and higher cognitive functions. We conclude that criticality is a key aspect of the operation of brains and it is a basic attribute of intelligence in animals and in man-made devices.

Freeman K Models: Structure and Functions
We propose a hierarchical approach to spatio-temporal neurodynamics, based on K sets. Low-level K sets were introduced by Freeman in the 70s, named in the honor of Aharon Kachalsky, an early pioneer of neural dynamics [18]. K sets are multi-scale models, describing increasing complexity of structure and dynamical behaviors. K sets are mesoscopic models, and represent an intermediate-level between microscopic neurons and macroscopic brain structures. The basic K0 set describes the dynamics of a cortical microcolums with about 10 4 neurons. K-sets are topological specifications of the hierarchy of connectivity in neuron populations. When first introduced, K-sets have been modeled using a system of nonlinear ordinary differential equations (ODE) [18,22]. K-dynamics predict the oscillatory waveforms that are generated by neural populations. K-sets describe the spatial patterns of phase and amplitude of the oscillations. They model observable fields of neural activity comprising EEG, LFP, and MEG. K sets form a hierarchy for cell assemblies with the following components [25]: K0 sets represent non-interactive collections of neurons with globally common inputs and outputs: excitatory in K0e sets and inhibitory in K0i sets. The K0 set is the module for K-sets.
KI sets are made of a pair of interacting K0 sets, both either excitatory or inhibitory in positive feedback. The interaction of K0e sets gives excitatory bias; that of K0i sets sharpens input signals.
KII sets are made of a KIe set interacting with a KIi set in negative feedback giving oscillations in the gamma and high beta range . Examples include the olfactory bulb and the prepyriform cortex.
KIII sets made up of multiple interacting KII sets. Examples include the olfactory system and the hippocampal system. These systems can learn representations and do match-mismatch processing exponentially fast by exploiting chaos.
KIV sets made up of interacting KIII sets are used to model multi-sensory fusion and navigation by the limbic system. KV sets are proposed to model the scale-free dynamics of neocortex operating on and above KIV sets in mammalian cognition.
K sets are complex dynamical systems modeling the classification in various cortical areas, having typically hundreds or thousands of degrees of freedom. KIII sets have been applied to solve various classification and pattern recognition problems [22,12,42]. In early applications, KIII sets exhibited extreme sensitivity to model parameters which prevented their broad use in practice [12]. In the past decade systematic analysis has identified regions of robust performance and stability properties of K-sets have been derived [63,38]. Today, K sets are used in a wide range of applications, including detection of chemicals [33], classification [12,28], time series prediction [5] and robot navigation [37,43].

Basic Building Blocks of Neurodynamics
The hierarchical K model-based approach is summarized in the 10 Building Blocks of Neurodynamics [18,21,23]: 1. Non-zero point attractor generated by a state transition of an excitatory population starting from a point attractor with zero activity. This is the function of the KI set.
2. Emergence of damped oscillation through negative feedback between excitatory and inhibitory neural populations. This is the feature that controls the beta-gamma carrier frequency range and it is achieved by KII having low feedback gain.
3. State transition from a point attractor to a limit cycle attractor that regulates steady state oscillation of a mixed E-I KII cortical population. It is achieved by KII with sufficiently high feedback gain.
4. The genesis of broad-spectral, aperiodic/chaotic oscillations as background activity by combined negative and positive feedback among several KII populations; achieved by coupling KII oscillators with incommensurate frequencies.
5. The distributed wave of chaotic dendritic activity that carries a spatial pattern of amplitude modulation AM in KIII.
6. The increase in nonlinear feedback gain that is driven by input to a mixed population, which results in the destabilization of the background activity and leads to emergence of an AM pattern in KIII as the first step in perception.
7. The embodiment of meaning in AM patterns of neural activity shaped by synaptic interactions that have been modified through learning in KIII layers.
8. Attenuation of microscopic sensory-driven noise and enhancement of macroscopic AM patterns carrying meaning by divergent-convergent cortical projections in KIV. 9. Gestalt formation and preafference in KIV through the convergence of external and internal sensory signals leading to the activation of the attractor landscapes leading to intentional action.
10. Global integration of frames at the theta rates through neocortical phase transitions representing high-level cognitive activity in the KV model.
Principles 1 through 7 describe the steps towards basic sensory processing, including pattern recognition, classification, and prediction, which is the function of KIII models. Principles 8 and 9 reflect the generation of basic intentionality using KIV sets. Principle 10 expresses the route to high-level intentionality and ultimately consciousness. The greatest challenge in modeling cortical dynamics is posed by the requirement to meet two seemingly irreconcilable requirements. One is to model the specificity of neural action even to the level that a single neuron can be shown to have the possibility of capturing brain output. The other is to model the generality by which neural activity is synchronized and coordinated throughout the brain during intentional behavior. Various sensory cortices exhibit great similarity in their temporal dynamics, including the presence of spontaneous background activity, power-law distribution of spatial and temporal power spectral densities, repeated formation of AM spatial patterns with carrier frequencies in the beta and gamma ranges, and frame recurrence rates in the theta range.
Models based on ODE and SDE have been used successfully to describe the mesoscopic dynamics of cortical populations for autonomous robotics [43,47]. However, such approaches suffered from inherent shortcomings. They falter in attempts to model the entire temporal and spatial range including the transitions between levels, which appear to take place very near to criticality. Neuropercolation and random graph theory offers a new approach to describe critical behavior and related phase transitions in brains. It is shown that criticality in the neuropil is characterized by a critical region instead of a singular critical point, and the trajectory of the brain as a dynamical systems crosses the critical region from less organized (gaseous) phase to more organized (liquid) phase during input induced destabilization and vise versa [31,30,27]. Neuropercolation is able to simulate these important results from mesoscopic ECOG/EEG recording across large spacial and temporal scales, as will be introduced in this essay.

Motivation of Neuropercolation Approach to Neurodynamics
We utilize the powerful tools of random graph theory (RGT) developed over the past 50 years [16,8,9] to establish rigorous models of brain networks. Our model incorporates cellular automata and percolation theory in random graphs that are structured in accordance with cortical architectures. We use the hierarchy of interactive populations in networks as developed in Freeman K models [18], but replace differential equations with probability distributions from the observed random graph that evolve in time. The corresponding mathematical object is called neuropercolation [41]. Neuropercolation theory provides a suitable mathematical approach to describe phase transitions and critical phenomena in large-scale networks. It has potential advantage compared to ODEs and PDEs when modeling spatio-temporal transitions in the cortex, as differential equations assume some degree of smoothness in the described phenomena, which may not be very suitable to model the observed sudden changes in neurodynamics. Neural percolation model is a natural mathematical domain for modeling collective properties of brain networks, especially near critical states, when the behavior of the system changes abruptly with the variation of some control parameters.
Neuropercolation considers populations of cortical neurons which sustain their metastable state by mutual excitation, and that its stability is guaranteed by the neural refractory periods. Nuclear physicists have used the concept of criticality to denote the threshold for ignition of a sustained nuclear chain reaction, i.e., fission. The critical state of nuclear chain reaction is achieved by a delicate balance between the material composition of the reactor and its geometrical properties. The criticality condition is expressed as the identity of geometrical curvature (buckling) and material curvature. Chain reactions in nuclear processes are designed to satisfy strong linear operational regime conditions, in order to assure stability of the underlying chain reaction. That usage fails to include the self-regulatory processes in systems with nonlinear homeostatic feedback that characterize cerebral cortices.
A key question is how cortex transits between gas-like randomness and liquid-like order near the critical state. We have developed thermodynamic models of cortex [30,27], which postulates two phases of neural activity: vapor-like and liquid-like. In the vapor-like phase the neurons are uncoupled and maximally individuated, which is the optimal condition for processing microscopic sensory information at low density. In the liquid-like phase the neurons are strongly coupled and thereby locked into briefly stable macroscopic activity patterns at high density, such that every neuron transmits to and receives from all other neurons by virtue of small-world effects [15,59]. Local 1/f fluctuations have the form of PM patterns that resemble droplets in vapor. Large-scale, spatially coherent AM patterns emerge from and dissolve into this random background activity but only on receiving a CS. They do so by spontaneous symmetry breaking [31] in a phase transition that resembles condensation of a rain drop, in that it requires a large distribution of components, a source of transition energy, a singularity in the dynamics, and a connectivity that can sustain interaction over relatively immense correlation distances. We conclude that the background activity at the pseudo-equilibrium state conforms to self-organized criticality, that the fractal distribution of phase patterns corresponds to that of avalanches, that the formation of a perception from sensory input is by a phase transition from a gas-like, disorganized, low-density phase to a liquid-like high-density phase [30].

Random Cellular Automata on a Lattice
In large networks, such as cortex, organized dynamics emerges from the noisy and chaotic behavior of a large number of microscopic components. Such systems can be modeled as graphs, in which neurons become vertices. The activity of every vertex evolves in time depending on its own state, the states of its neighbors, and possibly some random influence. This leads us to the general formulation of random cellular automata. In a basic two-state cellular automaton, the state of any lattice point x ∈ Z d is either active or inactive. The lattice is initialized with some (deterministic or random) configuration. The states of the lattice points are updated (usually synchronously) based on some (deterministic or probabilistic) rule that depends on the activations of their neighborhood. For related general concepts, see cellular automata such as Conway's Game of Life, cellular neural network, as well as thermodynamic systems like the Ising model [6,14,52,62]. Neuropercolation models develop neurobiologically motivated generalizations of cellular automata models and incorporate the following major conditions: Interaction with noise: The dynamics of the interacting neural populations is inherently non-deterministic due to dendritic noise and other random effects in the nervous tissue and external noise acting on the population. Randomness plays a crucial role in neuropercolation models and there is an intimate relationship between noise and the system dynamics, due to the excitable nature of the neuropil.
Long axon effects (rewiring): Neural populations stem ontogenetically in embryos from aggregates of neurons that grow axons and dendrites and form synaptic connections of steadily increasing density. At some threshold the density allows neurons to transmit more pulses than they receive, so that an aggregate undergoes a state transition from a zero point attractor to a non-zero point attractor, thereby becoming a population. In neural populations, most of the connections are short, but there are a relatively few long-range connections mediated by long axons. The effect of long-range axons are similar to small-world phenomena [61] and it is part of the neuropercolation model.
Inhibition: An important property of neural tissue that it contains two basic types of interactions: excitatory and inhibitory ones. Increased activity of excitatory populations positively influence (excite) their neighbors, while highly active inhibitory neurons negatively influence(inhibit) the neurons they interact with. Inhibition is inherent in cortical tissues and it controls stability and metastability observed in brain behaviors [40,38]. Inhibitory effects are are part of neuropercolation models.

Update Rules
A vertex vi ∈ V of a graph G(V, E) is in one of two states, s(vi), and influenced via the edges by d(vi) neighbors. An edge from vi to vj, vivj ∈ E, either excites and inhibits. Excitatory edges project the states of neighbors and inhibitory edges project the opposite states of neighbors, 0 if the neighbor's state is 1, and 1 if it is 0. Vertex's state, influenced by edges, is determined by the majority rule; when the most neighbors are active, the higher a chance for the vertex to be active, and when the most neighbors are inactive, the higher a chance for the vertex to be inactive. At time t = 0 s(vi) is randomly set to 0 or 1. Then, for t = 1, 2, ..., T − 1, a majority rule is applied simultaneously over all vertices. A vertex vi is influenced by a state of its neighbor vj ∈ N (vi), whenever a random variable R(vi, t) is less than the influencing excitatory edge vjvi strength, ωj,i, else the vertex vi is influenced by an opposite state of neighbor vj. If the edge vjvi is inhibitory, the vertex vj sends 0 when s(vj) = 1, and 1 when s(vj) = 0. Then, a vertex vi gets a state of the most common influence, if there is such, otherwise a vertex state is randomly set to 0 or 1, Fig. 1.1(a).
x f ollows majority more likely Formula for the majority rule:

Two-Dimensional Lattice with Rewiring
Lattice-like graphs are built by randomly re-connecting some edges in the graphs set on a regular 2-dimensional grid in and folded into tori. In the applied basic random rewiring strategy, n directed edges from n × 2 vertices are plugged out from a graph randomly. At that point, the graph has n vertices lacking an incoming edge and n vertices lacking an outgoing edge. To preserve the vertex degrees of the original graph, the set of plugged-out edges is returned back to the vertices with the missing edges in the random order. An edge is pointed to the vertex missing the incoming edge and is projected from the vertex missing the outgoing edge, Fig. 1.1(b). Such a two-dimensional graph models a KI sets with homogeneous population of excitatory or inhibitory neural population. It will be useful to reformulate the basic subgraph on a regular grid, as shown in Fig.  1.2. We list all vertices along a circle and connect the corresponding vertices to obtain the circular representation on the top right corner of Fig. 1.2. This circle is unfolded into a linear band by cutting it up at one location, as shown in the middle of the figure. Such representation is very convenient for the design of the rewiring algorithm. Note that this circular representation is not completely identical to the original lattice, due to the different edges when folding the lattice into torus. We have conducted experiments to compare the two representations and obtained that the difference is minor and do not affect the basic conclusions. An important advantage of this band representation that it can be generalized to higher dimensions by simply adding suitable new edges between the vertices. For example, in 2D we have 4 direct neighbors and thus 4 connecting edge, in 3D we have 6 connecting edges, and so on. generalization to higher dimensions is not used in this work as the 2D model is a reasonable approximation of the layered cortex.

Double-Layered Lattice
We construct a double-layered graph by connecting two lattices, see Fig. 1.3(a). The top layer G0 is excitatory, while the bottom layer G1 is inhibitory. The excitatory subgraph G0 projects to the inhibitory subgraph G1 through edges that excite, while G1 inhibitory layer projects toward G0 with edges that inhibit. An excitatory edge from G0 influence the vertex of G1 with 1 when active and with 0 when inactive. Conversely, an inhibitory edge from G1 influences the vertex of G0 with 0 when active and with 1 when inactive.

Coupling Two Double-Layered Lattices
By coupling two double-layered graphs, we have four subgraphs interconnected to form a graph G = G0 ∪ G1 ∪ G2 ∪ G3. Subgraph G0 and G1 are coupled into one double layer and G2 and G3 into the other. There are 3 excitatory edges between G0 and G2, G0 and G3, G1 and G2, G1 and G3, respectively, which are responsible to couple the two double-layered lattices, Fig. 1.3 Illustration of coupled layers with randomly selected connection weights between them. 1.3(a):Example of a double-layer coupled system that includes an excitatory layer G0 and an inhibitory layer G1. Within the layers the edges are excitatory; edges from G0 to G1 are excitatory and edges from G1 to G0 are inhibitory. 1.3(b): Example of two coupled double-layers, all together 4 coupled layers: G0 and G2 are excitatory layers, while G1 and G3 denote inhibitory layers. Each subgraph has three edges rewired. Edges from G1 to G0 and edges from G3 to G2 are inhibitory; the other edges are excitatory.

Statistical Characterization of Critical Dynamics of Cellular Automata
The most fundamental parameter describing the state of a finite graph G at time t is its overall activation level S(G, t) defined as follows: It is useful to introduce the normalized activation per vertex, a(G, t), which is zero if all sites are inactive and 1 if all sites are active. In general, a(G, t) is between 0 and 1: The average normalized activation over time interval T is given as: The latter quantity can be interpreted as the average mangetization in terms of statistical physics. Depending on their structure, connectivity, update rules, and initial conditions, probabilistic automata can display very complex behavior as they evolve in time. Their dynamics may converge to a fixed point, to limit cycles, or to chaotic attractors. In some specific models, it is rigorously proven that the system is bistable and exhibits critical behavior. Namely, it spends a long time in either low-or high-density configurations, with mostly 0 or 1 activation values at the nodes, before a very rapid transition to the other state [2,45]. In some cases, the existence of phase transitions can be shown mathematically. For example, the magnitude of the probabilistic component ω of the random majority update rule in Fig. 1.1(a) is shown to have a critical value, which separates the region with stable fixed point from the bistable regime [41].
In the lack of exact mathematical proofs, Monte Carlo computer simulations and methods of statistical physics provide us with detailed information on the system dynamics. The critical state of the system is studied using the statistical moments. Based on the finite size scaling theory, there are statistical measures which are invariant to system size at the state of criticality [7]. Evaluating these measures, the critical state can be found. It has been shown that in addition to the system noise, which acts as a temperature measure in the model, there are other suitable control parameters which can drive the system to critical regimes, including the extent of rewiring and the inhibition strength. [58].
Denote this quantity ω which will be a control parameter; it can be the noise, rewiring, connection strength or other suitable quantity. In this work, ω describes noise effects and related to the system noise level = 1 − ω. Distributions of a(G, ω) and a * (G, ω) depend on ω. At a critical value of ωi, the probability distribution functions of these quantities suddenly change. ωi values are traditionally determined by evaluating the 4th order moments (kurtosis, peakedness measure) u4, [7,17,49,50,51], but they can also be measured using the 3rd order moments (skewness measure) or u * 3 : According to finite-size scaling theory, for two structurally equivalent finite graphs G and G of different sizes, u4(G) = u4(G ) and u * 3 (G) = u * 3 (G ) for ω = ωi. However, for ω = ωi, u4(G) = u4(G ) and u * 3 (G) = u * 3 (G ). In this study, we use both the 3rd and 4th order momentums to determine the critical points of the system.

Dynamical Behavior of 2-D Lattices with Rewiring
Two-dimensional lattices folded into tori have been thoroughly analyzed in [45]; the main results are summarized here. In the absence of rewiring, the lattice is in one of the two possible states: ferromagnetic and paramagnetic, depending on the noise strength ω. The following states are observed: For ω < ω0, states 0 and 1 are equiprobable across the graph and a distribution is uni-modal. This state corresponds to paramagnetic regime.
For ω > ω0, one state dominates; vertex states are mostly 0 or mostly 1, and graph's a distribution is bi-modal. This is the ferromagnetic condition.
In the close neighborhood of the critical point ω0, the lattice dynamics is unstable. Immediately above ω0, the vertex mostly follows its majority influences.  For ω > ω0, the vertices cease to act individually and become a part of a group. Their activity level is determined by the population. This regime is defined by the condition that the vertices project more common activation than their receive, in the probabilistic sense. This transition is the first building block of neurodynamics and it corresponds to KI level dynamics. Fig. 1.4 illustrates critical behavior of a 2-D torus with ω0 ≈ 0.866; no rewiring. Calculations show that for the same torus having 5.00% edges rewired, the critical probability changes to ω0 ≈ 0.830. Rewiring leads to a behavior of paramount importance for future discussions: the critical condition is manifested as a critical region, instead of a singular critical point. We will show various manifestations of such principle in the brain. In the case of rewiring discussed in this section, the emergence of critical region is related to the fact that same amount of rewiring can be achieved in multiple ways. For example, it is possible to rewire all edges with equal probability, or select some vertices which are more likely rewired etc. In real life such as brains, all these conditions happen simultaneously. This behavior is illustrated in Fig. 1.5 and it has been speculated as a key condition guiding the evolution of the neuropil in the infant's brain in the months following birth [46,45,2] giving rise to a robust, broad-band background activity as a necessary condition of the formation of complex spatio-temporal patterns during cognition.

Narrow Band Oscillations in Coupled Excitatory-Inhibitory Lattices
Inhibitory connections in coupled subgraph G = G0 ∪ G1 can generate oscillations and multi-modal activation states [58]. In an oscillator, there are three possible regimes demarcated by two critical points ω0 and ω1: ω0 is the transition point where narrow-band oscillations start.
ω1 is the transition point where narrow-band oscillations terminate.
Under appropriately selected graph parameters, a(G, ω) distribution is uni-modal when ω < ω0, bi-modal when ω0 < ω < ω1, and quadro-modal when ω > ω1, Fig. 1.6. Each subgraph has 2304 (5%) edges rewired within the layer and 1728 (3.75%) edges rewired towards the other subgraph. Just below ω0, vertices may oscillate if they are excited. a(G, ω) values of temporarily excited graph oscillate, but return to the basal level in the absence of excitation. This form of oscillation is the second building block of neuron inspired dynamics. When ω is further increased, i.e., ω0 < ω < ω1, a(G, ω) exhibits sustained oscillations even in the absence of external perturbation. When the doublelayer is temporarily excited, it returns to the basal oscillatory behavior. Self-sustained oscillation is the third building block of neuron inspired dynamics and it corresponds to KII hierarchy.

Dynamics of Two Double Arrays
Two coupled layers with inhibitory and excitatory nodes exhibit narrow-band activity oscillations. These coupled lattices with appropriate topology can model any frequency range with oscillations occurring between ω0 and ω1. For example, an oscillator G(5%, 2.5%) = G0 ∪ G1 has 2.5% edges rewired across the two layers and 5% edges are rewired within each layer. G(5%, 2.5%) has ω0 ≈ 0.840 and ω1 ≈ 0.860, Fig. 1.7(b). Similarly, oscillator G(5%, 3.75%) has ω0 ≈ 0.845 and ω1 ≈ 0.865. Each of the coupled lattices in this example is of size 100×100. When coupling two such oscillators together, we have a total of 4 layers (tori) coupled into a unified oscillatory system. Clearly, the two oscillators coupled together may or may not agree on a common frequency. We will show that under specific conditions, various operational regimes may occur, including narrow band oscillations, broad-band (chaotic) oscillations, and intermittent oscillations, and various combinations of these. One of multiple ways to couple two oscillators is to rewire excitatory edges between subgraphs G0-G2, G0-G3, G1-G2, and G1-G3, respectively. Two coupled oscillators have additional behavioral regimes. When oscillators G(5%, 2.5%) and G(5%, 3.75%) are coupled with 0.625% edges into a new graph G(5%, 2.5%,. 3.75%, 0.625%), ω0 is shifted to lower level, Fig. 1.8. Every parameter describing a graph influences the critical behavior and under appropriate conditions. If two connected oscillators with different frequencies cannot agree on a common mode, but together they can generate large-scale synchronization or chaotic background activity.
In the graph made of two oscillators, there are four critical points ω0 < ω1 < ω2 < ω3, separating 5 possible major operational regimes.
Case of ω < ω0: The aggregate activation distribution is uni-modal and it represents the paramagnetic regime in lattice dynamics analogy.
ω0 < ω < ω1: The two coupled oscillators exhibit large-scale synchronization and the activation distribution is bi-modal. ω1 < ω < ω2: The two coupled oscillators with different frequencies cannot agree on a common mode, so together they generate aperiodic background activity.
ω2 < ω < ω3: only one oscillator oscillates in a narrow band and the activation distribution is okta-modal. ω > ω3: none of the components exhibit narrow-band oscillation and the activation distribution is a hexa-modal. This corresponds to ferromagnetic state.
The above regimes correspond to the fourth building block of neurodynamics as described by KIII sets.

Intermittent Synchronization of Oscillations in 3 Coupled Double Arrays
Finally, we study a system made of three coupled oscillators, each with its own inherent frequency. Such a system produces broad-band oscillations with intermittent synchronization-desynchronization behavior. In order to analyze quantitatively the synchronization features of this complex system with over 60,000 nodes in six layers, we merge the channels into a single two dimensional array at the readout. These channels will be used to measure the synchrony among the graph components. First we find the dominant frequency of the entire graph using Fourier transform after ensemble averaging across the ≈10,000 channels for the duration of the numeric experiment. We set a time window W for the duration of two periods determined by the dominant frequency. We calculate the correlation between each channel and the ensemble average over this time window W . At each time step, each channel is correlated with the ensemble over the window W . Finally, we determine the dominant frequency of each channel and for the ensemble average at time t and the phase shifts are calculated. The difference in phases of dominant frequencies measure the phase lag between the channel and the ensemble.
We plot the phase lags across time and space to visualize the synchrony between graph components, see Fig. 1.9. We can observe intermittent synchronization spatio-temporal patterns as the function of the ω, which represents the noise component of the evolution rule. For low ω values the phase pattern is unstructured, see Fig. 1.9(a). At a threshold ω value intermittent synchronization occurs Fig. 1.9(b). Further increasing ω, we reach a regime with highly synchronized channels. Intermittent synchronization across large cortical areas has been observed in ECoG and EEG experiments as well. The correspondence between our theoretical/computational results and experimental data is encouraging. Further detailed studies are needed to identify the significance of the findings.

Hebbian Learning Effects
According to principles of neurodynamics, learning reinforces Hebbian cell assemblies which respond to known stimuli by destabilizing the broad-band chaotic dynamics an leading to a attractor basin through a narrow-band oscillation using gamma band career   frequency. To test this model, we implemented learning in the neuropercolation model We use Hebbian learning, i.e., the weight from node i to node j is incrementally increased if these two nodes have the same state (1 − 1 or 0 − 0). The weight fromi tojj incrementally decreases if the two nodes have different activations (0 − 1 or 1 − 0). Learning has been maintained during the 20 step periods when input was introduced.
In the basal mode without learning, the 3 double layers can generate broad-band chaotic oscillations, see Fig1.10(a) where the lower row is zoomed-in version of top row. The spikes in Fig.1.10(a)   Coupling two KII oscillators into KIII, lowers ω0, the point above which large-scale synchronous oscillations are formed.
by flipping a few % of the input layer nodes to state 1 (active) for the duration of 20 iteration steps. During the learning stage, inputs are introduced 40 times (20 steps each), at regular intervals of 500 iteration steps. Without learning, the activity returns to the low-level chaotic state soon after the input ceases. Fig.1.10(b) shows that a narrowband oscillation becomes prominent during learning, when a specific input is presented. After learning, the oscillatory behavior of the lattice dynamics is more prominent, even without input, but the learnt input elicits much more significant oscillations. This is the manifestation of the 6th and 7th principles of Freemans neurodynamics, and it can be used to implement classification and control tasks using the neuropercolation model.
Further steps, toward the 8th and 9th principles involve the analysis of the amplitude modulation (AM) pattern, which is converted to a feature vector of dimension n. Here n is the number of nodes in our lattice which, after some course-graining, correspond to the EEG/ECoG electrodes. They provide our observation window into the operation of the cortex or cortex simulation. The feature vector is used to describe or label the cortical state at any moment. The AM pattern is formed by a low-dimensional, quasi-limit cycle attractor after synaptic changes with learning. The synaptic weights are described in a matrix, and the combination of different modalities to form Gestalts is done by concatenation of feature vectors. Research into this direction is in progress.

Conclusions
We employed neuropercolation model to demonstrate basic principles of neurodynamics. Neuropercolation is a family of stochastic models based on the mathematical theory of probabilistic cellular automata on lattices and random graphs and motivated by structural and dynamical properties of neural populations. The existence of phase transitions has been demonstrated in probabilistic cellular automata and percolation models. Neuropercolation extends the concept of phase transitions to large interactive populations of nerve cells near criticality. Criticality in the neuropil is characterized by a critical region instead of a singular critical point. The trajectory of the brain as a dynamical systems crosses the critical region from less organized (gaseous) phase to more organized (liquid) phase during input induced destabilization and vise versa. Neuropercolation simulates these important results from mesoscopic ECOG/EEG recording across large spacial and temporal scales.
We showed that neuropercolation is able to naturally model input-induced destabilization of the cortex and consequent phase transitions due to sensory input and learning effects. Broad-band activity with scale-free power spectra is observed for unlearnt conditions. On the other hand, Hebbian learning results in the onset of narrow-band oscillations indicating the selection of specific learnt inputs. These observations demonstrate the operation of the first 7 building blocks of neurodynamics. Our results indicate that the Learning effect in coupled oscillators. (a) Activity levels without learning; first row contains extended time series plots, 2nd row has zoomed in sections. Input spike is shown at every 500 steps. The activity returns to base level after the input ceases. 1.10(b): Activity levels with learning;The oscillations are prominent during learning and maintained even after the input step is removed (decayed).
introduced Hebbian learning effect can be used to identify and classify inputs. Considering that our model has 6 coupled lattices each with up to 10,000 nodes, massive computational resources have been involved in the simulations. Parallel chip implementation may be a natural way to expand the research in the future, especially when extending the work for multi-sensory Gestalt formation in the cortex.