Estimating complexity of spike-wave discharges with largest Lyapunov exponent in computational models and experimental data

1.   Introduction

Pathological states of physiological systems are considered to be different from the normal ones. The first studies of various biological rhythms using tools of nonlinear dynamics [1] showed that processes in the brain are normally characterized by the presence of irregular components with a high degree of complexity. Such a dynamics gives many functional advantages, since chaotic systems are capable of operating over a wide range of conditions, and thus easy to adapt to changes. Clearly expressed periodicity appears in aging and many pathological states is accompanied by a decrease in the degree of chaos and complexity [2].

In case of absence epilepsy the pathological synchronization of thalamo-cortical system is usually mentioned [3], [4]. Also, the transition to the synchronous state is usually considered to be a transition to more regular regime. The spike-wave discharges, which are the main encephalographic manifestation of absence epilepsy, are often considered as periodic oscillations as well [5], [6]. However, this regularity is not strict from the point of view of nonlinear dynamics.

The most usual and direct way to estimate, whether the current regime of behavior is regular or not, is to calculate the first (largest) Lyapunov exponent [7], which characterizes the behavior of the two initially very close points in a phase space [8], [9]. This was already done for other diseases like schizophrenia (see [10]) and for different normal physiological conditions (see [11], including sleep [12]). The distance between nearby trajectories varies exponentially with time, just with the magnitude of largest Lyapunov exponent. Therefore, if the largest Lyapunov exponent (Λ) is negative, the two initially close points of the phase space converge over time, which indicates the presence of a periodic regime. If these points are going apart, then the largest Lyapunov exponent is positive and this means that the system is in a chaotic regime.

The current work has two main objectives. First, to determine how much the SWDs differ from the baseline activity in term of regularity of oscillations and, therefore, in terms of predictability. This would help constructing, using and analyzing results of forecasting empirical models, which a popular tool both for automatic detection of SWDs and similar phenomena (the review of modern approaches one can find in [13]), and for coupling analysis (consider, e. g., Granger causality [14] with its adaptations to neurophysiological data [15], [16], transfer entropy [17] and partial directed coherence [18]). Second, to compare experimental data with simulations of SWD models in terms of signal regularity in order to determine how much these models can reproduce the observed dynamics and and whether they are suitable as a source of signals for testing time series analysis approaches.

2.   Method, models and data

Calculating Lyapunov exponent from equations is not a big problem. Estimating Lyapunov exponents from real data is more sophisticated and risky task. There are a number of known methods. The original Wolf's approach [19] was presented to estimate two largest Lyapunov exponents, but demands a sufficient amount of data (at least 2¹³ data points for Rössler system, with the main timescale being about 60 points length). The Eckmann's method [20] is able to estimate the whole Lyapunov spectrum, but is even more demanding for data amount. It is also known to work better for maps, rather than for systems with continuous times. The time series of spike-wave discharges are relatively short due to their average length of about 5–6 s [21]) and nonstationary at the same time, since spectral and coupling characteristics are changing over time from the seizure onset to its termination [22]. Therefore, the direct use of Wolf's and Eckmann's methods occurs to lead to unsafe estimates. So, we used the Rosenstein's et al. algorithm [23], which can be considered as a simplified version of the method [19] (only one largest Lyapunov exponent can be estimated) because this method requires fewer data, with being more noise robust.

2.1.   Rosenstein's algorithm in brief

We used the method proposed in [23]. This method is easy to implement and fast because it uses a simple measure of exponential divergence that circumvents the need to approximate the tangent map. Furthermore, the method is accurate for small data sets because it takes advantage of all available data. The first step of this method involves reconstructing dynamics on the attractor from a scalar time series. For this purpose we used the method of time delays with lag chosen equal to the first minimum of the autocorrelation function or a quarter of the characteristic oscillation period [24] and dimension estimated with false nearest neighbor algorithm [25].

Trajectory, as reconstructed by the method of delays, can be represented as a matrix in which each row is a state vector in a phase space, that is,

where x_n is a state of the system at time moment n. Using the original time series {xn}n=1N measured with sampling step Δt, the method of delays provides x_n = (x_n, x_n−l, …, x_n−(D−1)l), where D is an embedding dimension and l is time lag. Thus, X^ is the M × D matrix, where M = N − (D − 1)l. After the reconstruction of the phase space it is necessary to find the nearest neighbor for each point of the trajectory. We assume that the nearest neighbor to a vector x_n is another vector x_j, for which the Euclidean distance between them d_n,j is minimal at the additional assumption that these two vectors are not close in time, i. e. | j − n | > l in order not to consider consequent vectors as neighbors.

Then, the largest Lyapunov exponent can be estimated as a function of the average difference between the nearest neighbours. Following the definition of Lyapunov exponent the following formula describes the evolution of distance between vectors x_n and x_j:

where θ is the discrete time corresponding to the real time of trajectory evolution θΔt.

Theoretically, the equation (2.3) determines a set of approximately parallel lines, the slope of which is proportional to Λ. To obtain the robust estimate for Λ, let us average (2.3) over all possible n and consider the left side of the averaged equation (2.3) as a function of the discrete time θ:

The slope of φ(θ), if it is proposed to be a linear function, is Λ. Therefore, the approach is to fit φ(θ) to experimental data numerically. Note, that one can omit the second term in the equation (2.4) without any consequences since this term is constant and does not affect the slope.

Dependencies φ(θ) are evaluated for different embedding dimensions D (see Figure 1(a)) and, if it is possible to find linear parallel sites for nearby dimension values (see Figure 1(b)), the slope of the approximating line is an estimate of the largest Lyapunov exponent.

2.2.   Experimental data

Data were obtained from a previously used and published data set [26] in which male WAG/Rij rats (6–9 months) were used as experimental subjects. The Ethical Committee on Animal Experimentation of Radboud University Nijmegen (RU-DEC) approved the experiment. Local field potentials (LFPs) were recorded with a self-constructed electrode system for multi-site LFP recording at specified and verified brain locations. Stainless steel electrodes insulated with poliamide (ø = 127 µm) were fixated in a Teflon block, which contained small holes located at the relative A/P & M/L coordinates of the multiple electrode target structures as determined by the rat brain atlas of Paxinos and Watson [27]. Twelve electrode wires including those from reference and ground were glued to the teflon-block and fixed at the top-site to a connector pin, which was entered into an electrode pedestal suitable for the connection to a multi-lead electrode cable, which was connected to a swivel allowing long term recording in freely moving and well-adapted rats. LFP signals were amplified with a physiological amplifier (TD 90087, Radboud University Nijmegen, Electronic Research Group), filtered by a band pass filter in the range [1; 100] Hz with 50 Hz reduction by means of a Notch filter, and digitalised with a constant sample rate of 2048 Hz by WINDAQ-recording-system (DATAQ-Instruments). EEG of each rat was recorded for a period of 4 hours during the dark phase. In this work we studied recordings from the perioral region of the somatosensory cortex, layer 6 and two thalamic nuclei (the ventroposterial medial nucleus and the caudal part of the reticular thalamic nucleus) from 10 rats. Time series of 10 SWD fragments of 4 s length, taken starting immediately from the onset of SWD, were analyzed for each rat.

2.3.   Mathematical models of SWDs

2.3.1.   Neural mass model

The neural mass models are the most popular type of SWD models. Here, we consider one of the most recent of them [28], which is a simplification of previously published [29] model. In neural mass models each population is modelled as a lumped oscillatory system described by several differential equations. Four ordinary differential equations (ODEs) were used in the considered model (2.5): one for thalamic reticular cells (RE), one for cortical interneurons (IN), one for cortical pyramidal cells (PY) and one for thalamo-cortical cells (TC).

In this study the model is placed in a bistable state and was implemented using a neural population version of the Amari neural field equations [30]. The parameters were as follows τ_PY = 26, τ_IN = 32.5, τ_TC = 2.6, τ_RE = 2.6, h_PY = −0.35, h_IN = −3.4, h_TC = −2, h_RE = −5, C₁ = 1.8, C₃ = 1.5, C₉ = 1, C₂ = 4, C₇ = 1.5, C₆ = 0.6, C₈ = 3, C₄ = 10.5, C₅ = 0.2, ε = 2·10⁵, α = 2.8, b = 0.5. The equations were solved numerically with Euler method using constant time step 2⁻¹¹.

2.3.2.   Neural network model

The neural network model we used is a mesoscale model [31] which is a development of the previously published model [32]. Each neuron in the model actually represents a group of nearby neurons, which are divided into four populations. The populations are the same as in the model (2.5): PY, TC, IN and RE. Also, each population was split into 2 parts denoted as focal and surrounding. There were 200 of PY neurons, 120 of TC and RE neurons, and 50 of IN neurons in the used version of the model. Model neurons of each type were presented by the FitzHugh–Nagumo equations [33], [34]. Connections between neurons were set using the scheme, plotted in the Figure 2, which mostly follow the work by Suffczynski et al. [29].

Communication matrices were generated as follows. First, matrices of size N_PY = 40, N_IN = 10, N_TC = 40, N_RE = 40 were found, which demonstrated the appearance of high-amplitude generation in response to an attempt to initiate and subsequently maintain a discharge (a short-term gradual increase in the internal connections between neurons of the PY population, and then between neurons PY and RE). Such matrices were taken as a pathological subnetwork responsible for the appearance of spike-wave discharges in the model. Then, a matrix responsible for the normal dynamics of the system (i. e., not able to respond by generating high-amplitude oscillations to an attempt to cause them) was selected for each of the pathological matrix. The size of these matrices was N_PY = 160, N_IN = 40, N_TC = 80, N_RE = 80.

The integrated signal of neurons in the PY and IN populations is an analogue of the LFP signal of the neocortex, TC is a collection of cells from the VPM, and RE is from the RTN. The equations of the model were solved numerically by the Euler method with a step of 0.5.

3.   Results

3.1.   Largest Lyapunov exponent for experimental data

Consideration of series simulated from classical nonlinear systems showed the efficiency of the described method to be mainly determined by the embedding dimension D and the lag l [23]. When one selects these parameters by the above-mentioned methods, the curve for evaluation of the largest Lyapunov exponent should look like Figure 1 (a). However, processing experimental LFP signals showed that it is not the case, see Figure 3 (a).

It was found that a linear region (though not so long as for classical nonlinear systems) of the curves can be achieved (see Figure 3 (b)) between θ = 60 and θ = 100 (i. e. for θΔt ∈ [30;50] ms) for D ≥ 5, if lag is chosen in the following way. Let l′ be a lag, which is used for reconstruction of the phase space, and l is used for nearest neighbour search. Then, it is possible to achieve dependence φ(θ) as illustrated in Figure 3 (b), if l is equal to first minimum of mutual information, and l′ is significantly smaller, e. g. l′ = l/7. This can be explained by the fact that the signal has a number of independent time scales, some of which are much smaller than the main timescale — detailed spectra of LFPs can be found in many papers, see e. g. [21] or [35], where averaged spectra of both ictal and preictal activities were plotted in the same scale for different channels.

In total, we processed 300 four-second segments of LFP recordings at the beginning of SWD onsets: 10 time series from each of the 10 selected rats from 3 channels for absence activity and 10 time series from each of the 10 selected rats from 3 channels for baseline activity. Estimates of the largest Lyapunov exponent were averaged first over all discharges for each rat separately, and then over all rats. The average values of estimates for the largest Lyapunov exponent and their standard errors of the mean for all considered channels are shown in Figure 4 (a). The achieved values do not vary a lot across different animals: SEMs are small in comparison with the absolute values.

The first look at Figure 4 shows that baseline and SWDs are not much different, while Λ for baseline are always a little bit larger. If we consider the distribution of Λ estimates to be close to normal, we have to decide that the difference in Λ is significant (errorbars do not overlap) only for cortex for p-value p < 0.05.

3.2.   Largest Lyapunov exponent for the neural mass model

The model (2.5) allows to simulate both a background signal and SWDs at the same set of parameters due to switching between two coexistent attractors, but there is no possibility to simulate different animals, since different parameters of the model correspond to different severity of epilepsy. This limitation is general for neural mass models. Therefore, we generated 20 four-second epochs corresponding to normal activity and 20 epochs corresponding to epileptiform one from different initial conditions. The average values of estimates for the largest Lyapunov exponent and their standard errors of the mean for all considered channels are shown in Figure 4 (b). The parameters of the Rosenstein's method were chosen in the same way as for the experimental data, but separately.

We can see from the Figure 4 (b) that the Lyapunov exponents calculated for the model (2.5) are much smaller than for the experiment. Also, Λ values for SWDs are significantly smaller than for baseline activity. It is interesting that for the cortex this difference is mostly at the same level p = 0.05 as for experimental data, but for the thalamic nuclei this difference is much more significant (at the level p = 0.001) in contrary to experiment, where there is no significant difference between baseline and SWDs for both considered thalamic nuclei.

3.3.   Largest Lyapunov exponent for the network model

In the network model the dynamics is determined not only by the parameters of individual elements, but also by the communication matrices. There could be many matrices suitable for modeling, as shown in [32]. So, we can say that different matrices are models of individual animals. Therefore, the approach used for the network model is the same as we used for the experimental data: first we averaged the estimates of the largest Lyapunov exponent over all SWDs and background activity episodes for each matrix separately, and then we averaged over all matrices, calculating SEMs. The parameters of the Rosenstein's method were chosen in the same way as for the experimental data, but separately.

We can see from the Figure 4 (c) that Λ estimates obtained for the network model are smaller than for the experimental data, but there is no so large gap between Λ for SWDs and for baseline episodes as in neural mass model. In all cases SWDs are more predictable than baseline.

4.   Discussion and conclusions

The estimation of the largest Lyapunov exponent from the time series of both experimental data and the two considered numerical models turned out to be positive and rather large. In all cases Λ was larger for SWDs than for baseline activity, but for the experiment and for the network model this difference was rather small. High Λ for SWDs means that the signal during absence seizures is not so regular as usually considered. This additional complexity can appear due to the existence of another, more fast time scales in the signal, not related to the main one, which behave irregular [36].

High values of Λ for the baseline activity means that long range prediction of SWDs by means of any deterministic approaches is not possible. The Λ ∼ 16 corresponds to the Lyapunov time of 0.0625 s. If we consider that the measurement errors are at the level of discretization which is usually 16 or 24 bits in modern devices, the initial small perturbation can grow exponentially to the whole dynamical range at less than 1 s. However, this is not likely, since such a growth would take place only for small perturbations, as it is determined by the definition of Lyapunov exponent, and then larger and more regular time scales would start to dominate, providing lesser growth of perturbation. However, the predictability of SWDs is still very limited. The similar results were recently achieved by analyzing Shannon entropy, the prediction horizon was found to be not larger than 4 s [37]. In [16], [22] the changes in couplings were found to occur not earlier than 3.3 s before SWD onset.

Our analysis has shown Lyapunov exponents to be very close for all three considered channels, including both thalamic nuclei and somatosensory cortex. This is not surprizing, since the whole cortico-thalamo-cortical circuit is responsible for SWD generation. The same result can be achieved from the network model, with absolute values being smaller. But the neural mass model provides very different estimates for RTN and VPM during baseline. It also shows the significant difference between the cortex and the VPM during SWDs. The reason of this discrepancy is hard to understand in detail, but in general we can hypothesize that this would be a general, “genetic” problem of all neural mass models. All such models consider each type of cells as one solid oscillator. This means that the direct links between neurons, existing both in experiment and in network models, are replaced by coupling via mean field. These mean field couplings, however, cannot provide the same flexibility and connectivity as many individual direct links. Therefore, each oscillator occurs to be more independent and self-sustained in its dynamics than an ensemble of real or model neurons in the experiment or network models. There is also the problem of comparison between neural mass models and experimental data, since there is no direct possibility to provide a “population of models” similar to the population of rats in the experiment. In network models this problem can be solved via using different connectivity matrices.