**Extensive Four-Dimensional Chaos in a Mesoscopic Model of the Electroencephalogram**

#### 1.1 Background

Chaos that requires at minimum four degrees of freedom to be adequately described is known as “four-dimensional chaos” (FDC). Lorenz [2] was one of the first to describe this phenomenon in 1984, and his results were given a firm numerical underpinning by Sigeti [3]. Other examples of such behavior exist in the literature but are generally rarely reported [4], [5].

Despite the complexity of human cognition and behavior, the focus on nonlinearity and chaos to date in neuroscience has been on low-dimensional phenomena. Freeman [6] provides a rare exception to this, describing a hierarchy of attractors, from point attractor dynamics under deep anaesthesia, to chaos in perception. Nonetheless, the importance of attractor dynamics in cortical activity is well documented, with significant implications for the brain at rest and while performing tasks [7]–[10].

Previous work shows evidence for FDC in Liley’s mesoscopic model of the electroencephalogram, which is related to an inverse period doubling cascade [1]. That cascade also accounts for intermittent behavior, which is reminiscent of burst-suppression-like phenomena occurring in anaesthesia [11] and epileptic encephalopathies [12]. In this report, we extend those findings and reflect on the importance of high-dimensional chaos in mathematical neuroscience models.

#### 1.2 The Liley Model

Liley’s theory of neural dynamics is a spatiotemporal theory of the electroencephalogram (EEG). Its mesoscopic character implies that it does not focus on fine neuronal detail, instead concerning itself with the dynamics of populations of neurons. The macrocolumnar formulation of the model consists of ten first-order coupled nonlinear ordinary differential equations parameterized by a significant number of physiological constants, which describe the neural population properties in detail. For further discussion of this formulation and a complete derivation, see [1], [13]–[15].

#### 1.3 A Shortcut to Extensive Four-Dimensional Chaos Search

We assume that the *n* Lyapunov exponents of an *n*-dimensional dynamical system described by ordinary differential equations (ODE) are ordered and given by λ 1 ≥ λ 2 ≥ ⋯ ≥ λ n − 1 ≥ λ n

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i1.gif” alt=”View MathML“/>

. The largest integer *D* for which λ 1 + λ 2 + ⋯ + λ D ≥ 0

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i2.gif” alt=”View MathML“/>

is called the *topological dimension* of the attractor [16]. The next highest dimension ( D + 1

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i3.gif” alt=”View MathML“/>

) is the minimum integer dimension is which the attractor can exist.

Previously [1], the parameter set investigated provided the following first four *λ* values (base *e*, per second): λ 1 = 9.6

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i4.gif” alt=”View MathML“/>

, λ 2 = 0

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i5.gif” alt=”View MathML“/>

, λ 3 = − 6.4

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i6.gif” alt=”View MathML“/>

, and λ 4 = − 11.5

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i7.gif” alt=”View MathML“/>

. This means that λ 1 + λ 2 + λ 3 = 3.2

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i8.gif” alt=”View MathML“/>

, and λ 1 + λ 2 + λ 3 + λ 4 = − 8.3

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i9.gif” alt=”View MathML“/>

. So, in that case D = 3

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i10.gif” alt=”View MathML“/>

and D + 1 = 4

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i11.gif” alt=”View MathML“/>

, and the minimum integer dimension in which the attractor can exist is four: this is why the claim of FDC holds. It is important to realize the difference between the so-called hyper chaos, which has been studied extensively in maps and ODEs, and FDC. In FDC the topological dimension is at least three as in the case of hyper chaos, but FDC is characterized by only a single positive Lyapunov exponent, whereas hyper chaos is associated with two or more positive exponents (i.e. λ 1 > 0

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i12.gif” alt=”View MathML“/>

and λ 2 > 0

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i13.gif” alt=”View MathML“/>

).

Building on the insights from previous work, we here look for cases where λ 1 > | λ 3 |

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i14.gif” alt=”View MathML“/>

and λ 2 = 0

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i15.gif” alt=”View MathML“/>

. Since the exponents are ordered, we have λ 1 + λ 2 + λ 3 > 0

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i16.gif” alt=”View MathML“/>

and λ 1 + λ 2 + λ 3 + λ 4 < 0

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i17.gif” alt=”View MathML“/>

, with λ 3 < 0

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i18.gif” alt=”View MathML“/>

and λ 4 < 0

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i19.gif” alt=”View MathML“/>

. As such, *D* will always be at least three, and the minimum integer dimension that the attractor can exist in is at least four. Also note that the Kaplan–Yorke (or Lyapunov) dimension D KY

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i20.gif” alt=”View MathML“/>

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i21.gif” alt=”View MathML“/>

(1)

hence, with the above *λ*’s, D KY > 3

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i22.gif” alt=”View MathML“/>

always. In other words, for the purpose of our investigation it here suffices to examine the largest three Lyapunov exponents, dramatically decreasing the computational burden associated with the search. It is in fact possible to state that the chaos is at least four-dimensional, if λ 1 > | λ 3 |

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i23.gif” alt=”View MathML“/>

when λ 1 > 0

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i24.gif” alt=”View MathML“/>

and λ 2 = 0

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i25.gif” alt=”View MathML“/>

, without evaluating the full spectrum of ten exponents. This reduces the original computation from 120 coupled nonlinear ODEs to only 43. The calculation for the present report has been performed using the well-known Christiansen–Rugh algorithm for the partial Lyapunov spectrum [17], under the same boundary conditions and simulation lengths as discussed in previous studies [14].

#### 1.4 Extension of High-Order Chaotic Dynamics in the Liley Model

It is convenient to study the extension of the FDC respect to two important model parameters, p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i26.gif” alt=”View MathML“/>

and p e i

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i27.gif” alt=”View MathML“/>

, which are the excitatory and inhibitory input pulse densities to the modeled excitatory neural population. These are expected to vary most considerably and widely physiologically, capturing the effect of incoming thalamo-cortical input. The other physiological parameters of the model are the same as in Ref. [1].

Investigation of 433,316 different p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i28.gif” alt=”View MathML“/>

and p e i

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i29.gif” alt=”View MathML“/>

combinations, selected uniformly at random from a biologically relevant section of the p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i30.gif” alt=”View MathML“/>

– p e i

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i31.gif” alt=”View MathML“/>

plane (i.e. 0 < p e e ≤ 30

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i32.gif” alt=”View MathML“/>

, 0 < p e i ≤ 10

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i33.gif” alt=”View MathML“/>

), has been carried out. Selecting sets that show the largest Lyapunov exponent (LLE) being positive, 158,013 of the total points showed a LLE ≥ 1

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i34.gif” alt=”View MathML“/>

per second, base *e*. We use this threshold to avoid ambiguity with exponents that have a slow convergence to zero. Out of these clearly chaotic sets, 34,533 or 21.8 % of the chaotic instantiation of the model had λ 1 > | λ 3 |

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i35.gif” alt=”View MathML“/>

, exhibiting FDC, i.e. about 8 % of the overall points that have been simulated.

Figure 1 illustrates the behavior of the LLEs of the system as a function of p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i36.gif” alt=”View MathML“/>

and p e i

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i37.gif” alt=”View MathML“/>

. The color scale is such that blue represents an LLE ≤ − 20

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i38.gif” alt=”View MathML“/>

, red means LLE ≥ 20

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i39.gif” alt=”View MathML“/>

and green corresponds to a value of approximately zero. Again, this is a consequence of the slower rates of convergence for exponents associated with periodic orbit dynamics, which do not always correspond to exactly zero at the end of the simulation run. The extensive nature of chaos in this plane is evident, as are the limit cycles ( LLE = 0

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i40.gif” alt=”View MathML“/>

, green) and point attractors (negative LLE, blue). Note also the periodic windows interspersed among the chaotic fringes.

**Fig. 1.** Largest Lyapunov exponent values (LLE) plotted as a function of p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i41.gif” alt=”View MathML“/>

and p e i

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i42.gif” alt=”View MathML“/>

. *Blue* is for LLE ≤ − 20

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i43.gif” alt=”View MathML“/>

, *green* is for LLE ≈ 0

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i44.gif” alt=”View MathML“/>

, and *red* is for LLE ≥ 20

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i45.gif” alt=”View MathML“/>

, all values per second, base *e*. Thalamic inputs are expressed in units of ms ^{−1}

Superimposing the high-order chaotic points (in black), Fig. 2 gives an illustration of the extension of the FDC region. The overall set of FDC points is made of three different areas: one very extended and of irregular shape at large values of p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i46.gif” alt=”View MathML“/>

, one limited and more regularly shaped at intermediate values of p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i47.gif” alt=”View MathML“/>

and a very small, elongated collection of points at low p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i48.gif” alt=”View MathML“/>

, p e i

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i49.gif” alt=”View MathML“/>

, far from the chaotic (red) region and bordering with the sea of periodic activity (in green). This region is particularly interesting, since it corresponds to small values of the thalamic input, which are associated with ordinary thalamic activity. In our previous work, instead, we reported FDC for a parameter set with a high value of p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i50.gif” alt=”View MathML“/>

, possibly corresponding to a pathological or other abnormal state.

**Fig. 2.** LLE plotted as a function of p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i51.gif” alt=”View MathML“/>

and p e i

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i52.gif” alt=”View MathML“/>

, with four-dimensional chaotic points in *black*. The parameter set discussed in our previous work [1] is indicated by *a white circle* and corresponds to p e e = 24.2453 ms − 1

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i53.gif” alt=”View MathML“/>

, p e i = 2.299 ms − 1

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i54.gif” alt=”View MathML“/>

Using the bifurcation package AUTO [18], a continuation in the two model parameters for the saddle-node on invariant cycle (SNIC) points of the periodic orbits allows for some speculation on the genesis of such intriguing, high-order phenomena. Figures 3 and 4 show a partial and full analysis, in conjunction with the plot of LLEs. Two 1:1 resonance points [19], highlighted in red in Fig. 3 and light blue in Fig. 4, overlap with sections of the FDC area. Lines of SNICs also correspond to the finger-shaped gaps (in green) for LLE ≈ 0

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i55.gif” alt=”View MathML“/>

. Further strong resonances often occur at the end of such gaps. The complex bifurcation scenario that may arise from such strong resonances could be partially responsible for the FDC, in particular for the high-order chaos occurring at low p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i56.gif” alt=”View MathML“/>

values. Generation of FDC has also been associated with an inverse period doubling cascade in our previous work, and it could be possible that similar mechanisms control FDC for high p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i57.gif” alt=”View MathML“/>

values. A thorough analysis of the processes involved is beyond the scope of this report.

**Fig. 3.** Continuation of codimension-two points associated with SNICs in the full p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i58.gif” alt=”View MathML“/>

– p e i

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i59.gif” alt=”View MathML“/>

plane, i.e. not limited to biologically relevant values. *The inset* shows the accumulation of SNICs lines in the area of high-order chaos at low values of p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i60.gif” alt=”View MathML“/>

, p e i

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i61.gif” alt=”View MathML“/>

. Strong resonance occurring inside the four-dimensional chaotic set are shown as *red squares*. Other 1:1 resonances are shown as *white squares*

**Fig. 4.** Overlapping of SNICs lines (*light blue*) with LLE plotted as a function of p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i62.gif” alt=”View MathML“/>

and p e i

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i63.gif” alt=”View MathML“/>

. *Points in black* are where the chaos is four-dimensional. The strong resonances inside the chaotic areas (in the vicinity of ( 3.5 , 1.5 )

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i64.gif” alt=”View MathML“/>

and ( 19.5 , 0.2 )

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i65.gif” alt=”View MathML“/>

) are now indicated by *light blue squares*

#### 1.5 Conclusions

The results presented here show the extent of four-dimensional chaos within a biologically relevant parameter slice of Liley’s model. The fact that FDC is so extensive suggests extreme caution when performing visual inspections or even nonlinear time-series analyses to match experimental EEG traces with brain states or functions [20], [21]. In fact, the chaotic attractor associated with FDC dynamics for a specific parameter set at high p e e

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i66.gif” alt=”View MathML“/>

shows an amorphous appearance [1], which may make it hard to distinguish from noise. Combinations of noise and cortical activity already appear to be very difficult to untangle for dynamics simpler than the one discussed here [22].

A novel, important finding is that high-dimensional chaos is not limited to pathological or abnormal brain states but is present also for values of the thalamic input well inside the ordinary range of thalamic activity, i.e. 0 < p e e < 10 ms − 1

<img class=”mathimg” src=”http://www.mathematical-neuroscience.com/content/inline/s13408-015-0028-3-i67.gif” alt=”View MathML“/>

. Hence, activity associated with high-dimensional strange attractors could occur more frequently than so far assumed. This aspect has also relevance as regards multistable behavior, given that Liley’s model can support multistable dynamics induced by different attractors [15]. Multistability has in fact been shown to capture aspects of brain activity in a variety of important neurological settings, including perceptual decision making and critical behavior of the brain at rest [23]–[25]. We hope that our findings may inspire further research work into the role of high-order chaotic dynamics in brain activity.

**Source: **Extensive Four-Dimensional Chaos in a Mesoscopic Model of the Electroencephalogram

**Via:** Google Alert for Neuroscience