Directional functional coupling during limbic seizures in rats revealed by nonlinear Granger causality

1. Introduction

Spontaneous seizures involving the hippocampus and parahippocampal structures are the hallmark of temporal lobe epilepsy. Epileptic activity may propagate widely along synaptic pathways and can become highly synchronized between limbic and other regions [1, 2]. In the most widely used experimental model of temporal lobe epilepsy, spontaneous limbic seizures develop after an initial pharmacologically or electrically induced status epilepticus [3]. However, in non-epileptic rats, limbic seizures may also occur during chronic exposure to antagonists of cannabinoid receptor type 1 (CB1) [4].

The endocannabinoid system participates in the regulation of brain sensitivity to epileptic activation [5, 6]. In the pilocarpine rat model for epileptogenesis CB agonists reduce early post-status epilepticus seizure manifestations and subsequent mortality [7]. Moreover, in the same model, agonists abolish the occurrence of late spontaneous epileptic seizures [5]. CB1 antagonists on the other hand diminish the brain resistance to seizures [5]. It was found that they facilitate the spreading of seizures from the midbrain to limbic structures in rats prone to midbrain-driven audiogenic seizures [8]. What is more, healthy rats treated long term with CB antagonist become prone to spontaneous seizures [4]. These seizures are topic of the present paper.

Endocannabinoids are synthesized on demand after activation of postsynaptic neurons and they retrogradely suppress neurotransmitter release through presynaptic metabotropic CB1 receptors. In this way, endocannabinoids contribute to maintaining the physiological level of neuronal excitability and synaptic function [9]. Epileptic excitation quickly activates this defense mechanism against overexcitation [10]. Endocanabinoid CB1 receptors [11, 12] play the critical role in the endocannabinoid-dependent protection against seizures, and local deletion of these receptors in the hippocampus strongly exacerbates seizures [13].

Widespread extralimbic regions, including both cortical and subcortical structures, are increasingly recognized to be involved in secondary generalization of limbic seizures in epileptic rats [14, 15] and complex partial seizures in patients with temporal lobe epilepsy (TLE) [15, 16]. The present study was aimed to analyze directional functional coupling between the neocortex, hippocampus, thalamus and midbrain during ictal period of SLV326-induced limbic seizures. The hippocampus is particularly prone to epileptic excitation. The midbrain region contains triggering areas for certain seizure types, primarily for generalized tonic-clonic seizures [17-19]. The thalamus is an important brain site for seizure propagation and its connectivity with the hippocampus and cortex is altered in patients with TLE [15, 20, 21].

Establishing results of the study on time series analysis, it is important to use methods which are capable to resolve a number of issues: 1) nonlinearity, 2) coupling directionality, 3) cause of similarity between signals, including direct links, common source, driving mediated through some intermediate structure, or even simple random coincidence of oscillations [22]. In this study the time varying adapted nonlinear Granger causality method [23] was chosen. This technique takes into account features of experimental data: time scales, nonlinearity, effective dimension. Therefore, it has a good temporal resolution with fine sensitivity and specificity, as it was shown [24, 25], being mostly insensitive to mediated coupling [24].

2. Material and Methods

2.1. Subjects and local field potential recording

The seizures analyzed in the present study occur during chronic treatment with the CB1 antagonist SLV326, as it was reported previously in [26]. Since the method to elicit the seizures was reported extensively there, here it will be repeated only in short.

The study was performed in accordance with the guidelines of the European Community for the use of experimental animals and was approved by the ethical committee for animal studies (RUDEC-2007-161).

The electroencephalographic (EEG) experiments were performed on 24 Crl:WI Wistar rats (Charles River Laboratories, Sulzfeld, Germany), aged 7 months. From an age of 8–9 weeks, SLV326 was administered daily by oral gavage, 2-3 mg/kg, dissolved in a semi-solid solution. Solutions of SLV326 were provided by Solvay Pharmaceuticals (Weesp, Netherlands).

The 24 animals were surgically provided with two tripolar EEG electrodes (Plastics One MS-332/2-A) under complete isoflurane anesthesia. Electrodes were placed in the following brain structures (distances provided in mm from bregma: anterior, lateral, and depth): frontal cortex (FC): +2, -2, -1; hippocampus (HP): -4.2, -3.6, -4.1; thalamus (TH): -2.6, -2.7, -7.3; and midbrain (MB): -8.8, -1.7, -5.2. Ground and reference electrodes were placed over the cerebellum bilaterally, with reference electrode located on the side ipsilateral to the recording electrodes. Animals were allowed to recover for at least two weeks.

LFP and video recordings were made during 24 hours. Local field potential (LFP) signals were amplified, filtered between 1 and 100 Hz and digitalized at the sample rate equal to 512 Hz using the Windaq system (DATAQ Instruments, Akron, Ohio, USA).

In the 24 hour recording period, 6 out of 24 rats demonstrated generalized limbic seizures. In the present study, 30 seizures recorded in three rats were analyzed, because for these rats all electrodes demonstrated the appropriate signal level during the whole recording.

2.2. Application of time-variant adapted nonlinear Granger causality to LFP data

To analyze the coupling between four measured signals from different brain areas, the adapted nonlinear time varying Granger causality method was used 23. The method is based on the construction of two predicting models: a univariate model, which predicts values from the series  based on only itself with the average error (1), and a bivariate one, which uses also data from the series  (its prediction error (2) is usually less). The prediction improvement (3) is a main coupling measure and lies in between 0 and 1.

where  is a measured value,  is a corresponding value predicted with a univariate model, and  is a value predicted with a bivariate model.

Absolute values of  are usually not informative, as it was shown by Smirnov and Mokhov [27]. However, increase or decrease of  makes sense under the condition that the operator of the evolution of the considered systems did not change structurally (but its parameters could change [24]). Changes in coupling strength can be detected this way, if Granger causality is used in a moving time window [28]. Here, a time window of length 1 s (512 data points) was used with window overlap of 0.9 s. The method parameters, including time scales were adjusted in accordance with recommendations and criteria developed in [25, 29, 30].

2.3. Mutual information function in moving window

The mutual information function (MI) was calculated for the same channel pairs and in the same conditions as for the Granger causality, including time window length and shift, and averaging, by means of the method proposed in [31]. Note, that MI is an undirected nonlinear measure of simultaneous similarity in time series. This means that it cannot show the direction of driving or reveal the reason of changes, as almost all simple measures such as cross-correlation (linear or nonlinear) or coherence. The MI was applied to scalar series. While there are fine modern works [32], showing the possibility to improve the measure performance significantly by using the vector series with non-uniform embedding, this demands too many data to be applied in a moving window. Due to the same reason, the transfer entropy [33, 34] was not used.

2.4. Statistical analysis of coupling estimates

For every seizure the individual dependencies PI(t) and MI(t) were calculated, including 10 s before and 10 s after the seizure. Accordingly, each considered EEG fragment containing a limbic seizure was split into three parts:

1. Initial: 10 s prior to and 10 s after seizure onset: pre-ictal period and 10 s of the first high frequency stage;
2. Middle: 10 s prior and 10 s after the transition from the first high frequency stage to the second low frequency stage;
3. Final: 10 s before and 10 s after the seizure termination: low frequency stage and postictal period.

The resulting dependencies PI(t) and MI(t) were averaged across all seizures, matching onset, start of high frequency stage and seizure termination. Then, for each averaged dependency PI(t) and MI(t), the background levels PI_bg and MI_bg, respectively, were established as an average value over the 3 s time interval (baseline period, from 10 to 7 s before the seizure onset). Using the mean baseline levels, the normalized dependencies were calculated as PI₀(t) = PI(t) – PI_bg and MI₀(t) = MI(t) – MI_bg. The value of PI₀=0 (and MI₀=0) corresponds to the baseline level; positive values of PI₀ and MI₀ correspond to a larger coupling than in baseline and negative — to a lower one.

The values of PI(t) and MI(t) obtained for different seizures for the same time point were used as a sampling, and a single sampling t-test was performed to identify the statistical difference of the mean value from PI_bg or MI_bg respectively. Since this test was performed for every time point, a Bonferroni-like correction for multiple testing was implemented: the achieved p-values were multiplied by a number of independent time intervals, from which PI or MI were calculated. If after correction a resulting p-value occurred to be less than 0.05, the results were considered as significantly different from baseline and plotted in color (blue if PI<PI_bg or MI<MI_bg and red if PI>PI_bg or MI>MI_bg) in Figures 2 and 3. Otherwise, the results were plotted in gray as insignificant.

3. Results

3.1. Time-frequency analysis

All seizures started suddenly, with a sharp spike of small or medium amplitude detected in all recorded channels.  The total duration of the seizures was from 35.0 s to 115.3 s with the mean duration of 53.1 s. Two main stages can be distinguished (Figure 1):

1. a high frequency stage, starting from a frequency of about 20 Hz and decreasing to 15 Hz in the first 5 s;
2. a low frequency stage, at which the main seizure frequency dramatically falls to 2 Hz.

The high frequency stage lasted for 23.3–55.7 s, with a mean duration of 36.0 s. The amplitude of oscillations was rising during this stage. The whole stage was very nonstationary. Well pronounced higher frequencies can be seen in different time epochs.

This low frequency stage lasted for 5.7-59.6 s, with the mean longitude equal to 17.1 s. The high frequency dynamics of the previous stage was also present in the first part of the second stage, albeit less pronounced. The low frequency oscillations were losing amplitude gradually and stopped suddenly. This power fall maked a very sharp marker of the seizure termination.

Figure 1. Local field potentials recorded during a limbic seizure and their spectrograms: (a, b) — from frontal cortex, (c, d) — from thalamus, (e, f) — from hippocampus, (g, h) — from midbrain. Three black vertical lines mark: 1) start of seizure at time point 10 s, 2) margin between two stages  at time point 38 s , and 3) end of seizure at 50 s.

3.2. Granger causality analysis

Results of the coupling analysis, both using Granger causality and mutual information function are shown on Figures 2 and 3.

Figure 2. Dependencies of mean normalized prediction improvement PI₀ and mutual information function MI₀ on time for pair FC-MB, FC-Th, FC-Hp, calculated in 1 s time window. Red and blue dots indicate values significantly (p<0.05 with Bonferroni correction) different (red – larger, blue – smaller) from baseline level, gray dots — not different. Black dashed vertical lines indicate seizure onset, transition to the second stage and termination, gray line before them — length of time window. Baseline, first (high frequency) stage, second (low frequency) stage and postictal period are subscribed on the top of each subfigure.

Figure 3. Dependencies of mean normalized prediction improvement PI₀ and mutual information function MI₀ on time for pair Th-MB, Hp-MB, Hp-Th, calculated in 1 s time window. Red and blue dots indicate values significantly (p<0.05 with Bonferroni correction) different (red – larger, blue – smaller) from baseline level, gray dots — not different. Black dashed vertical lines indicate seizure onset, transition to the second stage and termination, gray line before them — length of time window. Baseline, first (high frequency) stage, second (low frequency) stage and postictal period are subscribed on the top of each subfigure.

3.2.1. Seizure initiation

About 2 seconds before the start of a seizure, a drop in coupling is seen in all considered channel pairs. This decoupling continues during at least first 10 seconds of the high frequency stage (see Figure 2, blue dots) except hippocampus – midbrain pair, in which coupling restores to normal level at the seizure onset. A new bidirectional drop in coupling in this pair starts about 6 seconds after seizure onset.

 3.2.2. Transition to low frequency stage

The driving from the FC to the HP remains lower than baseline during the whole seizure (Figure 2c upper panel). For all other channel pairs the coupling restores in the course of the high frequency stage with some global transient increases: unidirectional from the MB and HP to the FC (Figure 2a and Figure 2c middle panels) and from the TH to the MB (Figure 3a upper panel), bidirectional in the HP–MB loop (Figure 3b upper and middle panel).

Towards the transition from the high frequency stage to the low frequency one, all couplings are fairly equal to baseline except from for the FC–TH loop, which decreases about 5 seconds before the transition (Figure 2b, upper and middle panel), and except the above mentioned FC–HP loop, which shows decreased coupling during the whole seizure (Figure 2c, upper panel).

During the first few seconds of the low frequency stage a lower than baseline coupling is seen from FC and HP to MB (Figures 2a and 2e, upper panels) and from TH to FC (Figure 2b, middle panel). The couplings from FC and TH to HP are also lower than baseline, but these decouplings last the whole low frequency stage (Figure 2c, upper panel, and Figure 3c, middle panel).

3.2.3. Seizure termination

Towards the end of the seizures, some global transient increased couplings are observed – the same increases which were observed towards the end of the high frequency stage: unidirectional from the MB and HP to the FC (Figure 2a and Figure 2c middle panels), and from the TH to the MB (Figure 3a upper panel), but the bidirectional increase in the HP–MB loop present in the high frequency stage, is now only significant for the direction from the HP to the MB, (Figure 3b upper panel). The MB to HP increase is absent, while this increase was very pronounced in the end of high frequency stage (Figure 3b middle panel). Moreover, the drivings from FC and TH to HP are still low (Figure 2c, upper panel, and Figure 3c, middle panel).

When the moving window covers the seizure termination, there is a sharp increase in coupling in all pairs of channels in which the thalamus is involved as a driven structure (Figures 2b, 3a and 3c). However, such an increase was mentioned previously as a method artifact [24]. In the postictal stage all couplings are generally rather low compared to baseline.

3.3. Mutual information analysis

Results of the mutual information analysis are shown in the lower panel of each subfigure. The curves of normalized mutual function MI₀(t) vastly oscillate and most values are insignificant. Distributed along the whole seizure period, for some time points, the mutual information is significantly higher than baseline for the FC–MB pair (Figure 2a), while the MI is smaller than the baseline for the FC–HP and TH–HP pairs (Figures 2c and 3c). In the 5 s before the termination of the seizure, a number of values of MI are significantly higher than baseline for the TH–MB pair (Figure 3a).

4. Discussion

The complete scheme of complex changes in the network was plotted in Figure 4. The main outcomes are formulated as follows:

1. During limbic seizures two stages were distinguished: the high frequency (15–20 Hz) stage, followed by the low frequency (2 Hz) stage;
2. About 2 seconds before the seizure start, there is a drop in coupling in all considered channel pairs. This decoupling continues during at least first 10 s of the high frequency stage;
3. During the whole seizure the driving from the frontal cortex to the hippocampus remains decreased. For all other channel pairs, the coupling restores during the high frequency stage with some global transient increases;
4. In the postictal period all couplings are rather low compared to baseline.

Figure 4. Schemes of significant changes in coupling at different stages of seizure. Blue/red arrows show significant decay/increase of prediction improvement PI₀, significant changes in the mutual information function are shown near the arrows by signs MI↑ (red) and MI↓ (blue).

4.1. Discussion on Methods

Results of coupling analysis from complex time series are always method dependent [22]. Simple methods, such as linear correlation and coherence function can be easily calculated using standard toolboxes such as MATLAB or SciPy. However, it is generally understood that these methods are not entirely adequate due to the following reasons: 1) only linear similarities can be analyzed, 2) coupling direction cannot be revealed, 3) the cause of similarity remains unknown (unidirectional driving, bidirectional one, common source, mediated driving through some intermediate structure, or simply a random coincidence of oscillations). Some nonlinear measures, like mutual information function [31], nonlinear correlation coefficient [35], and phase synchronization index [36] overcome the first problem, but not the others.

The mutual information function was already used for studying mechanisms of epilepsy in animals [37] and humans [38] before. The time shift between series at which the measure demonstrated the maximum was considered as a way to determine the coupling direction. However, recent studies clearly showed that undirected by nature nonlinear measures of similarity, including phase synchronization indeces [39] and mutual information [40] are unsafe for determining the coupling directionality, since very often the results occur to be random and misleading even for the relatively simple simulated series. In the current study, the mutual information function calculated using the most advantageous approach [31] was applied but with very limited success: a lot of results are either insignificant or cannot be linked with the results of the spectral analysis or Granger causality.

The most advantageous methods, allowing to address all three mentioned problems, are Granger causality [41], partial directed coherence [42], phase dynamics [43], and transfer entropy [33]. All these methods are more or less parametric. Generally, parametrization sufficiently decreases data requirements. In principle, it is possible to extract information about directed coupling from relatively short series: 30 oscillations for phase dynamics [43] and 4–8 oscillation for Granger causality [23]. However, the risk to obtain false positive (insufficient specificity) and false negative (insufficient sensitivity) results due to wrong parametrization also rises [25, 44, 45]. Till now, there is no possibility to obtain an absolute warranty against false results; one can only reduce risks by means of adaptation of used methods to specifics of studied data.

In the present study the pairwise analysis was preformed. It is not safe, since direct and indirect couplings cannot be completely separated. But, the classical conditional Granger causality [41] as well as more advanced approaches [46], which aim to eliminate the redundancy, cannot be applied directly, since they demand too many data. For example, in [46] 4,000 data points (10 s of recording, 400 Hz sampling frequency) were used with linear Granger causality, while to study limbic seizures one has to operate with nonlinear causality (so, much more coefficients have to be estimated), otherwise the method specificity becomes very low, as it was shown in refs 24, 41, and 40. Also, relatively short epochs, not larger than 1 s, have to be considered due to signal nonstationarity. However, the intracranially recorded local field potentials are much less redundant than surface EEGs, especially in the case if they were obtained from different and distinct brain structures as in the present study.

4.2. Discussion on coupling changes during seizures

The performed analysis was based on averaging over 30 seizure recorded in three SLV-treated rats (the whole data set from 3 animals). In [47] the large variability of results of statistical models constructed from EEGs was found. Therefore, in additiona to already reported results the results of Granger causality were also averaged independently for individual animals, and animal based results were found to follow the whole set based averaged results, but with less significance. Therefore, the conclusions were made based on the whole data set.

Reduced driving from the frontal neocortex to subcortical structures and increased driving from subcortical sites to the neocortex, described during generalized limbic seizures in the present study, are in line with the results of the imaging studies in both epileptic rats 48 and patients with TLE [15, 16]. These studies have shown that limbic seizures are accompanied by activation of subcortical structures and neocortical deactivation. It has been hypothesized that abnormally increased activity of the thalamus and upper brainstem prevents normal activation of the cortex and leads to abnormally reduced function of the fronto-parietal cortex, which underlies abnormal motor behavior and impaired consciousness during complex partial seizures [14, 15]. Our results show that, despite involvement of the frontal cortex in the seizure expression, it plays only a passive role, being driven by subcortical structures.

The main active players seem to be the hippocampus and midbrain. Both structures show intrinsic epileptogenicity and contain triggering areas for certain seizure types: the hippocampus for limbic seizures, the midbrain for reflex audiogenic seizures [18, 19]. These two structures send direct and indirect glutamatergic projections to the cortex, thalamus and to each other. The hippocampus and midbrain may represent seizure drivers and their reciprocal excitatory interaction during the limbic seizures may represent a positive feedback loop driving the seizures.

While the central role of the hippocampus in TLE is well known, a role of the midbrain in limbic seizures is surprising. Nevertheless, participation of the midbrain in mechanisms of secondary generalized limbic seizures [49] and primary generalized tonic-clonic seizures [17] has been shown. It has been previously shown that chronic treatment with another CB1 receptor antagonist (rimonabant) significantly facilitates seizure spreading from the midbrain to cortex during audiogenic kindling [8]. It is interesting to investigate, whether the active participation of the midbrain in limbic seizures is a common phenomenon, or it is observed only during limbic seizures induced by cannabinoid antagonists.

5. Conclusion

The results of our study indicate that the cortico-thalamic loop, underlying the maintenance of absence seizures [50, 51], is decoupled during the whole generalized limbic seizure, that corresponds to reduced thalamo-cortical functional connectivity in patients with TLE [20]. On the other hand, during limbic seizures the frontal cortex is driven by the hippocampus that is in contrast to cortico-hippocampal decoupling during cortico-thalamic absence seizures [52].

Conflict of interest: none declared.

Funding

This research was funded by Stipendium of President of Russian Federation for support of young scientists СП-3605.2018.4 and Russian Foundation for Basic Research, Grants No. 17-02-00307 and 18-015-00418.

Ethical approval

All applicable international, national, and institutional guidelines for the care and use of animals were followed.