Abstract
In the epileptic brain, phase amplitude cross-frequency coupling (CFC) features have been used to objectively classify seizure-related states, and the inter-seizure state has been demonstrated as being random, in contrast to the seizure state being predictable; however, the excitatory and inhibitory networks underlying their dynamics remain unclear. Therefore, the objectives of this study are to classify the dynamics of seizure sub-states labeling seizure-like event (SLE) onset and termination intervals using CFC features and to obtain their underlying excitatory/inhibitory cellular correlates. SLEs were induced in mouse neocortical brain slices using a low-magnesium perfusate, and were recorded in Layer II/III using simultaneous local field potential (LFP) and whole-cell voltage clamp electrodes. Classification of onset and termination of SLE transitions was investigated using CFC features in conjunction with an unsupervised two-state hidden Markov model (HMM). γ-Distributions of their durations indicated that both are predictable. Furthermore, omitting 4 Hz from the HMM classifier switched both SLE sub-states from statistically deterministic to random without changing the dynamics of the SLE state. These results were generalized to 4-aminopyridine (4-AP)-induced SLEs and human seizure traces. Only during these sub-states, both excitatory and inhibitory currents coupled with the field. Where excitatory currents phase locked to a broad range of frequencies between 1 and 12 Hz, inhibitory currents dominantly phase locked at 4 Hz. We conclude that inhibition underlies the predictability of neocortical CFC-defined SLE transition sub-states.
Significance Statement
To date, the underlying excitatory and inhibitory bases of objectively identified seizure sub-states have not been determined. Here, we identify these sub-states using phase amplitude cross-frequency coupling (CFC) features. This is tested in both human and rodent seizure models. Then, using simultaneous field and whole-cell recordings in a mouse model, we investigate the potential for explaining seizure sub-state dynamics from interactions of the excitatory and inhibitory currents with local field network activity. We found that the frequency ranges at which these currents are coherent with field oscillations have a dramatic effect on the predictability of seizure state sub-states. Such information is critical for identifying the mechanisms underlying the dynamics of epileptic seizures.
Introduction
Intractable epilepsy patients make up 30% of clinical cases, necessitating further investigation into potential treatment options (Kwan and Brodie, 2000). In studies of the outcomes from surgical resection for these intractable epilepsy patients, seizure freedom after surgery was markedly lower for neocortical epilepsy patients as compared to those whose seizures originated from the temporal lobe (Spencer and Huh, 2008; de Tisi et al., 2011). Yet, neocortical circuitry is not examined as widely as temporal lobe structures, which could contribute to why there are fewer successful targeted therapies for these patients.
An epileptic seizure is a period of complex electrical brain signals characterized by excessive pathologic neuronal activity (Fisher et al., 2014). Previously, the waveforms of high- and low-frequency oscillations (HFOs and LFOs) were examined independently as biomarkers of the seizure-prone brain (Ren et al., 2011; Jacobs et al., 2012); however, the pathologic rhythms were not separable from the physiologic. Yet, this separation is crucial for identifying underlying local networks involved in the pathology. Phase amplitude cross-frequency coupling (CFC) between LFO and HFO is a marker of neuronal excitability and communication (Lakatos et al., 2005; Mormann et al., 2005; Cohen et al., 2009; Canolty and Knight, 2010; Tort et al., 2010). Yet in the epileptic brain, it has been found that the ranges of coupled frequencies are different (Nariai et al., 2011; Guirgis et al., 2013, 2015; Samiee et al., 2018). There is, in particular, a marked increase in δ (0.5–4 Hz)-HFO CFC, which has been used to successfully identify resection zones for surgery and classify the seizure state (Ibrahim et al., 2013; Guirgis et al., 2013, 2014).
The seizure state is composed of multiple sub-states that include both the onset and termination transitions, but the underlying network processes that govern these sub-states remain unclear. Based on γ-distribution fits of interseizure and seizure duration histograms, seizure onset has been projected as randomly occurring (Suffczynski et al., 2006), whereas seizure termination has been demonstrated to be predictable (Bauer et al., 2017). However, there may be features that govern seizure sub-states that have predictable underlying patterns for both the onset and termination of the ictal events (Jacobs et al., 2018). In particular, using CFC to classify the ictal state may be more accurate and provide greater insight into state dynamics than visual inspection or amplitude-based measurements. Therefore, one of the main aims of this study was to identify the underlying dynamics of seizure sub-states using CFC features in the epileptic cortex.
A second aim of this study was to identify the underlying excitatory and inhibitory contributions to the CFC-identified seizure sub-states. During the onset of neocortical seizure-like events (SLEs), there is a large burst of inhibition synchronized with excitation (Kawaguchi, 2001; Timofeev and Steriade, 2004; Trevelyan et al., 2006; Lévesque et al., 2016; Librizzi et al., 2017; Chang et al., 2018). Furthermore, at termination, there may be a weak recovery of inhibitory currents, with a concomitant slowing down of large excitatory field bursts (Köhling, 2014). Yet, it is unknown how these excitatory and inhibitory patterns relate to the dynamics of the seizure sub-state durations. Furthermore, although combinations of LFO and HFO have been previously used to classify multiple seizure state sub-states (Guirgis et al., 2014), only extracellular recording techniques were applied. Hence, their underlying excitatory and inhibitory currents were not investigated. Therefore, simultaneous extracellular and intracellular measurements provide a substantial advancement in identifying the underlying cellular mechanisms of seizure sub-states.
In brief, our first aim is to identify the relation between CFC-identified seizure sub-states and the predictability of these classified intervals, and our second aim is to investigate the excitatory and inhibitory currents associated with these CFC features. To achieve this, simultaneous measurements of extracellular local field potential (LFP) and whole-cell recordings were made in the superficial layers of the mouse somatosensory cortex. A novel CFC based feature set was applied to a two-state hidden Markov model (HMM), and was used to classify neocortical SLE states, with onset and termination sub-states indicated by the marginal errors in the HMM. This was tested in low Mg2+, and 4-aminopyridine (4-AP), and then the results were generalized to human intracranial EEG (iEEG) signals. Finally, the excitatory and inhibitory current patterns underlying the different HMM-identified SLE state transitions were determined.
Materials and Methods
Animals
C57/Bl6 (RRID:IMSR_JAX:000664) male and female mice aged P10 to P16 (Charles River Laboratories, RRID:SCR_003792) were humanely killed for these experiments in accordance with the Canadian Animal Care Guidelines and with the United States Public Health Service’s Policy on Humane Care and Use of Laboratory Animals. All surgical procedures were approved and done in accordance with the guidelines of the Animal Care Committee of the University Health Network.
Cortical slice preparation
To prepare the cortical tissue for resection, mice were anesthetized using 50-mg/kg pentobarbital. The pedal reflex was used to assess the depth of anesthesia. Once deeply anaesthetized, the mice were decapitated, and then their whole brain was swiftly removed. The brain was placed in an ice cold, oxygenated sucrose solution made up of 248 mM sucrose, 2 mM KCl, 3 mM MgSO4, 1 mM CaCl2, 26 mM NaHCO3, 1.25 mM NaH2PO4, and 10 mM D-glucose (Florez et al., 2015). Then, to limit the contribution of distal projections to induced epileptiform activity, the lateral 20% of both cerebral cortices was excised, and the remainder of the cortex was fixed to a vertical block using a cyanoacrylate adhesive. Coronal slices (500 µm) were sectioned using a Leica 1200 V vibratome. Then, each slice was hemisectioned and transferred to a solution of artificial CSF (ACSF; 95% O2, 5% CO2) containing 123 mM NaCl, 25 mM NaHCO3, 10 mM glucose, 3.5 mM KCl, 1 mM MgSO4, 1.2 mM NaH2PO4, and 1.5 mM CaCl2. The slices were incubated for 30 min at 32 ± 0.5°C followed by 1 h at room temperature before the start of the experiments.
Electrophysiology
The slices were transferred individually to a submerged recording chamber and were perfused with ACSF (95% O2, 5% CO2; 10 ± 1 ml/min; 34 ± 0.5°C). Glass electrodes (1.5 mm, World Precision Instruments) were used for both LFP and whole-cell recordings. The LFP recording electrodes (∼2-MOhm resistance) contained ACSF and were placed within 200 µm of the intracellular electrodes in cortical Layers II and III. An Olympus BX51 microscope (OLY-150IR camera–video monitor unit) was used with an infrared filter to visualize the neurons. The whole-cell recording electrodes (3–5 MOhm resistance) contained a solution of 135 mM K-gluconate, 1 mM MgCl2, 10 mM NaCl, 2 mM Na2ATP, 0.3 mM NaGTP-Tris, 10 mM NaHEPES, 0.5 mM EGTA, and 0.0001 mM CaCl2; pH 7.2–7.3. Signal acquisition and storage were performed using an amplifier (Multiclamp 700B), a digitizer (Digidata 1322A) and PClamp software (version 10.2; Molecular Devices/Molecular Devices Corporation). Whole-cell recordings were done in voltage clamp. Neurons were held at –70 mV for spontaneous EPSC (sEPSC) acquisition and 0 mV for spontaneous IPSC (sIPSC) acquisition. In 3/6 experiments for sIPSCs, CsCl replaced K-gluconate to prevent potassium current contamination of the intracellular signal. In these cases, putative pyramidal neurons were solely identified based on their appearance under the infrared filter. Furthermore, the addition of the CsCl did not affect the temporal correlation between the inhibitory currents and field activity; therefore, all inhibitory current data were combined for signal analysis.
Putative pyramidal neurons were identified based on their appearance under the infrared filter and their electrophysiological features. A series of hyperpolarizing and depolarizing current pulses (200-ms duration) were injected through the whole-cell recording electrode. We computed the following to confirm putative pyramidal neuron identity: an average input resistance of 150.75 ± 10.63 MΩ (n = 16 cells) was calculated from the slope of the current−voltage curve. A spike half-width of 1.59 ± 0.07 ms was obtained from the first spike at rheobase, and defined as the duration of the action potential at the membrane voltage 50% from threshold to peak. Finally, a membrane time constant of 19.61 ± 1.41 ms was obtained from fitting an exponential to the first 50 ms of the minimally injected current. These values reported are consistent with those reported in the literature (Avermann et al., 2012; van der Velden et al., 2012; Tyler et al., 2015). However, we did not correct for liquid junction potential since all results central to the hypothesis were measured against a change from baseline in voltage clamp. To assure only healthy neurons were used, only those whose resting membrane potential before SLE induction greater than –65 mV were studied; furthermore, any recording with a change in the input resistance >20% during the inter-SLE phase relative to the ACSF baseline was omitted from further analysis.
Tissue viability was assessed in two ways. First, we observed that all slices exhibited greater than one SLE. Second, SLE amplitude did not decrease over time. Therefore, only tissue with sustained intrinsic seizure control mechanisms and good slice viability was used for this study (Voss et al., 2013). For SLE induction, MgSO4 was omitted during the preparation of the ACSF (low Mg2+), or 4-AP (100 µM) was applied to ACSF containing normal levels of MgSO4 (Wang et al., 2016).
Patient description
Nineteen seizures recorded from iEEG grids placed on the neocortex were obtained from seven patients with intractable epilepsy. Seizures were obtained from patients at Toronto Western Hospital (Toronto, Canada) and Phramongkutklao Hospital (Bangkok, Thailand). Data from these patients were used in previous studies examining HFO coherence, CFC in seizure onset zone classification and seizure predictability (Cotic et al., 2015; Guirgis et al., 2015; Li et al., 2016; Jacobs et al., 2018). Patients and the ethics committees of the affiliated institutions provided informed consent and approval. All patients were seizure free (Engel class 1) after surgery. For each seizure, a bipolar recording from an 8 × 8 iEEG multi-electrode recording taken from electrodes nearby the resected zone was chosen for the analysis. This onset zone had been previously defined by two expert neurologists. These recordings were acquired at sampling frequencies of 500, 1000, and 2000 Hz. The data obtained from Canada was notch-filtered at 60 Hz and its harmonics up to Nyquist frequency, and the data obtained from Thailand were notch-filtered at 50 Hz and its harmonics up to Nyquist frequency. Before feature extraction, all signals were decimated to 500 Hz.
Data analysis
CFC-based feature set
For an objective way to classify the SLE state and SLE sub-states, a novel CFC feature set was applied to a two-state HMM. First, all LFP recordings were decimated to 1000-Hz sampling rate then FIR notch filtered at 60 Hz and its first five harmonics. The instantaneous phase of the LFOs (1–12 Hz) and amplitude of the HFOs (30–250 Hz), were computed by convolution of the filtered LFP signal with a complex Morlet wavelet as defined by:
(1)
(2)
These were computed at fb = 3-Hz bandwidth, with fc = 1-Hz center frequency. To avoid edge effects, a 5-s window on either side of the signal of interest was kept during the computation of the continuous wavelet transform (cwt), and then subsequently removed for calculation of the CFC indices (ICFC). Then the ICFC was estimated on 2-s sliding windows. For every 2-s window, the instantaneous phase signal was divided into 20° equal intervals. In each phase interval, the amplitude, Ahigh, was averaged then normalized by the sum of the average amplitudes (Tort et al., 2010):
(3)where N is the number of phase intervals and j is a discrete phase interval. The Kullback–Leibler (KL) distance was computed between the Anorm and a uniform distribution. Then, the KL distance was normalized to an interval between 0 and 1. This normalized KL distance was the phase amplitude CFC index (Icfc; also known as modulation index) for each data point on the CFC image. Then, the ICFCs used for training and testing the HMM were obtained from the ICFC matrix on first modifying the vertical axis to be logarithmic, then re-shaping the image to a 4 × 4 matrix (Fig. 1; features 1–16). The remaining features were the derivatives of the initial set of indices [Fig. 1B, features (1’)–(16’)]. The features used for training were obtained using 2-s windows for the ICFC computation and 0% overlap on the sliding window function. The features used for testing were obtained similarly, using 2-s windows and 95% overlap.
Building the feature set for SLE state classification using ICFC matrices. A, LFP recording from the neocortex showing a representative SLE induced in low-Mg2+ ACSF. Regions expanded for comparison are highlighted by gray bars. A1, Expanded pre-SLE region unfiltered (upper), and filtered (lower). The lower trace was band pass filtered at 3–6 and 30–40 Hz, indicated by the black and superimposed gray lines, respectively. This range was chosen to represent the frequency space shown to have highest ICFC values in that instance. The right image is the ICFC matrix of a 2-s window including the expanded region, demonstrating low ICFC over the full frequency spectrum, with peaks at 3–6 and 30–40 Hz. A2, Expanded region of the SLE state unfiltered (upper), and filtered (lower). The lower trace was band pass filtered at 1–3 and 30–40 Hz, indicated by the black and superimposed gray lines, respectively. This range was chosen to represent the frequency space shown to have highest ICFC values in that instance. The right image is the ICFC matrix of a 2-s window including the expanded region, demonstrating peak ICFC over a small range of the frequency spectrum. The numbers on the image indicate the features unwrapped to make the feature set. B, The feature set obtained from the signal in A used for classifying the SLE state. The numbers on the vertical axis correspond to the frequency ranges indicated by the numbers on the pseudocolor plot in A2 (1–16) and their time derivatives [(1’)–(16’)]. Each row was normalized to the maximum and minimum ICFC values across the full trace.
ICFC control analysis
To identify significant ranges of CFC, surrogate data were generated by converting the original time series into a complex series using the cwt (1- to 12-Hz intervals in 1-Hz increments). The phase was then random block-point shuffled with blocks of 1% of the sampling frequency (10 data points, randblock.m). Then the CFC indices were computed on 5-s windows at the classified ro onset of the SLEs and classified re end of the SLEs. The amplitude coefficients of the high-frequency cwt (30–250 Hz, in 5-Hz increments) were used, with the surrogates, to compute the CFC indices. This was repeated 200 times. The original CFC was compared to the surrogate data either at onset or termination. The CFC was z score normalized to the mean and standard deviation of the distribution of the surrogate values. Ranges considered significant were those whose distributions were above 3× the SEM of the surrogate dataset. This test for significance provided the boundaries of significant CFC for further control analysis.
HMM state classification
We used a first order Q = two-state HMM which, for our application, delineated SLE and inter-SLE states S1 and S2 at a time t, where 1 ≤ t ≤ T, and T indicates the length of the feature set (Guirgis et al., 2014). The HMM requires estimating the following distributions {s, A, B}: for i = 1,2 and j = 1,2 (1) s = Pr(q1 = j); the probability of initially being in state Q = qj. (2) At = Pr(qt = j | qt-1= i); the probability of transitioning from one state to another at time t. (3) Bt = Pr(Xt = x | qt = j) = {bj(O)}; the probability of observing x, a given value of the feature set, when in the state q at a given time. O is a 2-dimensional set of observed features for state j. The model values for {s, A, B} were initially estimated using the k-means clustering algorithm, and then their accuracy was enhanced using the expectation-maximization (EM) algorithm. To compute B, each {O} set of features was fit using a mixture of Gaussians with the probability density function (PDF) defined as
(4)where µm is the mean, CVm the covariance and m = 1 or 2 is the number of Gaussian basis functions. B = {bj(O)} is the sum of the weighted Gaussian probability densities over all the basis functions.
The marginal posterior distributions (MPDs) of the HMM were computed as follows:
(5)
(6)
(7)where α is the joint probability of observing all data from time 1:t; and β is the conditional probability of all data from t + 1 to T at state i. Then, the log likelihood evaluated the goodness-of-fit between the feature set of training observations O, and the model parameters {s, A, B}. The process was repeated to the point where either the difference between log likelihood values of successive iterations dropped below 10−5 or 100 iterations took place. Generally, the log likelihood converged within 11–25 iterations and in no case did it reach 100 iterations. The stem plots are the MPD(1) and MPD(2) in testing conditions, which was re-named as S1 and S2 to represent the inter-SLE state and SLE state, respectively.
A γ-distribution parameterization of time histograms was performed using the gamfit function in MATLAB. The γ-distribution is characterized by a shape and a rate parameter. The maximum likelihood estimates of the shape parameter, α, led to the following inferences of the dynamics of the SLE state intervals (Bauer et al., 2017):
(5)
The SLE state S2 was subdivided into putative SLE transitions within state S2 at onset (ro) and end (re) of the SLE state (rs). These time windows were obtained from 0 → 1 and 1 → 0 fluctuations of the MPD S2. Duration and intensity of the SLE state S2 was computed (Fig. 2B,C) using the following operations: Duration(S2) = L(ro) + L(rs) + L(re); Intensity(S2) = (VL)2/L; where VL is the LFP signal within a time interval, L, of the sub-sections of S2 (ro, rs, re). A Pearson correlation coefficient was applied to measure the linear correlation between the duration of the sub-intervals of the SLE state S2 and that of the duration and intensity of the SLE state S2 (0 → no correlation, 0.1–0.3 → weak correlation, 0.31–0.69 → moderate correlation, 0.7–1 → strong correlation).
ICFC features classify low-Mg2+-induced SLE sub-states. A, LFP recording from neocortical Layer II/III showing a representative SLE induced in low-Mg2+ ACSF, and the MPDs S1 and S2. Boxes indicate putative SLE transitions at onset (ro) and end (re) of the SLE state (rs). These time windows are obtained from 0 → 1 and 1 → 0 fluctuations of the MPDs S1 and S2. B, Scatter plots showing a moderate positive correlation between the duration of ro, rs, and re and the duration of the full SLE state. Each data point is one SLE. C, Scatter plots showing a moderate positive correlation between state rs and the ro transition, and the intensity of the SLE state and weak correlation between the transitions re and the intensity of the SLE state. For Pearson correlation coefficients, see Table 1, lines 1–6. D, Histograms of the durations of ro, rs, and re with their respective γ-distribution fits. See Table 2 for shape parameters of the γ-distribution fits and average durations of each interval. All distributions fit to a shape parameter >1.
Performance of the HMM classifier was determined using receiver operator characteristic (ROC) curves. To compare to a standard approach in defining seizure onset and termination, these were identified using the LFP power. The LFP power was computed by squaring the absolute value of the cwt and averaging over the full SLE state. The SLE onset was defined as the time where the LFP signal rose to 3× the baseline signal. The SLE termination was defined as a return to baseline after a final field bursting event, followed by >10-s refractory period. If the area-under the curve (AUC) is 1, that is a perfect classification, whereas if it is 0.5 or less, it is a worthless classification (for an example, see Fig. 3A). During the inter-SLE interval, false positive (Fp) was obtained if, for each point on the MPD of state S2, state S2 was classified as higher probability than state S1, and the opposite was considered true negative (Tn). During the SLE interval, true positive (Tp) was when S2 remained in higher probability than S1, and if the probability distribution switched, this was false negative (Fn). The threshold for defining the SLE interval was varied between 0 and 1 on S2 to obtain the points on the ROC curve. Sensitivity, specificity and a specificity factor, ξ, were computed as follows:
(6)
(7)
(8)
Reconstructed data containing low-frequency activity from low-Mg2+ SLEs poorly classifies the SLE state. A, Example of LFP trace of a low-Mg2+ SLE (top). The signal power is used as an alternate definition of the SLE state (bottom) and is applied to the computation of ROC curve. Gray box indicates defined SLE state. a1, MPD of S2 of the LFP presented in A. a2, MPD of linear interpolated reconstructed data obtained from the envelope of the SLE in A. a3, MPD of cubic interpolated reconstructed data obtained from the envelope of the SLE in A. B, One example of a single burst from the SLE in A at 118 s under normal conditions (the unfiltered dataset), using a linear interpolant (reconstructed data) and using a cubic interpolant (reconstructed data). C, Example of the power spectra of the three datasets. The low-frequency power remains similar across models, but the high-frequency power is minimized. D, ROC curves for the three conditions for all SLEs tested. For a summary of the ROC characteristics, see Table 4 (n = 46 reconstructed SLEs, 14 slices, 10 subjects). Inset is a representation of the computation of the ROC sensitivity and specificity. The square wave shows a typical MPD of the HMM. The lines are depictions of the stem plots of the MPD of state S2. Fp is false positive, Tp is true positive, Fn is false negative, Tn is true negative.
Optimal sensitivity and specificity values were also reported from the specificity factor ξ-exponential curve fit. These were reported as the nearest point on the exponential function relative to the coordinate (0,1).
Controls for ICFC classified using HMM
To account for the possibility that the irregular shape of the low-frequency signal alone is sufficient to capture the dynamics of the SLE state using CFC features, data with minimal HFO contribution was simulated by reconstructing the original datasets. We obtained the envelope function of the original LFP dataset (envelope.m), and then added Gaussian white noise (wgn.m) of the same variance as the original dataset computed in a 10-s window of quiescence during the pre-SLE state. To account for variability in the wave form shapes during SLEs, this was repeated using two variations of the envelope interpolation: linear and cubic.
To identify whether large spiking transients would affect the classified state transitions, data were first bandpass filtered 30–250 Hz using 10,000 order finite impulse response (FIR) filters that were run in the forward and reverse direction (filtfilt.m). Then, a spike melting algorithm (spikinator.m) was applied to smoothly remove all spikes above a specified threshold (Zalay et al., 2009) using the following parameters: 10-ms window, threshold of 1.2, amplitude of 0.01, smoothing of 1, degree of spike cleaning of 0.8, and window frame margin of 0.1. These parameters were chosen on inspection of the amplitude of the large spiking transients. This high-frequency amplitude signal was used in conjunction with the phase of the low-frequency 1- to 12-Hz oscillation and feature sets were built as previously defined. We only used the CFC features directly (features 1–16, no derivatives) for this test. Each new signal was tested on an HMM that was trained on its corresponding original signal. The % difference between the S2s was defined as the sum of the difference between the MPD of S2 of the original signal and the MPD of S2 of the spike reduced signal, divided by the total number of points tested.
Charge transfer energy of the cellular currents
To evaluate the charge transfer energy of the sEPSCs and sIPSCs, the trapezoidal rule was applied. Because each SLE sub-state was of different duration, data were evaluated in amortized time. The SLE sub-state duration was divided into 10 time intervals of equal duration, and the energy was averaged over each 10% of the 100% of sub-state duration. All data were statistically tested relative to the 10% before the onset transition sub-state, ro.
Peak phase coherence between the cellular currents and the LFPs
For identifying the correlation between the timing of the synaptic events and the phase of the LFP oscillation, first, the timing of spontaneous PSCs (sPSCs) was identified using an amplitude threshold and the findpeaks.m function in MATLAB, with a minimum peak distance of 50 ms and minimum peak height of 3× SD of the 10% pre-SLE area. Then, at each time of the sPSC peak, the phase at which the event occurred was computed as the instantaneous phase of the cwt, as previously defined. Phase vector plots include the phase distribution of the peak sPSCs, mean phase direction (bold black lines), mean phase vector length (number at top) and scale (relative probability value). The mean phase vector direction is the mean phase at which the peaks are occurring. The mean phase vector length is the degree of variability in the timing of the peaks around the mean; therefore, the higher the mean phase vector length, the stronger the phase locking. The mean phase vector length is the phase locking value (PLV) used in a recursive algorithm.
Using a recursive algorithm in MATLAB, we varied the threshold for peak identification and calculated the phase of the low-frequency signal (between 1 and 12 Hz) until we observed significant elevated phase locking during the SLE transition sub-state. For each phase vector plot, the Rayleigh test was used as a statistical test of circular uniformity (CircStats toolbox), with a p < 0.05 threshold for significance. The algorithm included a criterion to identify PLVs >0.3 (indicating moderate phase locking) and p < 0.05. For visualization purposes, the filtered version of the LFP at the frequency range of interest was plotted superimposed on a stem plot where the stems are located at each detected sPSC.
Experimental design and statistical analysis
Statistical analyses were performed using MATLAB and Microsoft Excel. An unpaired t test was used when there were two groups and a normal distribution, Wilcoxon rank-sum test when there were two groups and non-normal distribution. A Bonferroni post-hoc correction was applied in each case for tests of multiple comparisons. Pearson Correlation was used for analyses of correlation between two variables. Details of numbers of SLEs, slices and subjects evaluated are reported in appropriate subsections of figure legends. The α threshold was set to 0.05 (for further details, see Table 1).
Statistics table
Results
CFC classifies distinct SLE sub-states
First, we asked whether CFC features were able to classify SLE sub-states in the neocortex of mice. SLEs were induced using low-Mg2+ ACSF and the ICFC matrices were computed on the resulting signal (Fig. 1A). The ICFC values between 1- to 12-Hz low-frequency and 30- to 250-Hz high-frequency were reduced, normalized, then applied to a two-state HMM. To distinguish state S1, the inter-SLE state, and state S2, the SLE state, we computed the MPDs of the HMM. Spontaneous SLEs began 5.4 ± 0.8 min after the addition of low-Mg2+ ACSF, and recurred at 0.61 ± 0.13 SLEs per minute (for Ns, see Table 2).
Summary of duration distributions for SLEs and seizures
The width of the MPD is characterized by marginal errors in the HMM. We examined these errors at the onset and termination of the SLEs, and hypothesized that the SLE sub-states are the timeframe of these errors (Fig. 2A). To test this hypothesis, we correlated the duration of the errors with state S2 defined SLE duration and intensity (Fig. 2B,C). See Table 2 for the average duration of the onset interval ro, SLE state rs, and end interval re. The sub-state intervals had positive moderate correlations with the SLE state S2 duration, whereas the SLE state rs had a positive strong correlation (Table 1, lines 1–6). In contrast, intensity of the SLE state S2 was moderately correlated with ro and rs. Hence, SLE severity was correlated primarily with state rs; however, the duration of the SLE was correlated with all three intervals. These data supported the hypothesis that MPD values between 0 and 1 at the onset and termination of the events were transition sub-states into and out of the SLE state, rs. Hence, CFC between 1–12 and 30–250 Hz classified distinct SLE sub-states.
SLE sub-states have deterministic underlying dynamical processes
To identify the type of dynamics underlying the SLE sub-states, we computed the γ-distribution fits of the distribution of the durations (Fig. 2D; Table 2). The shape parameter for the SLE state rs was >1, suggesting that this classified state was statistically deterministic. Furthermore, both ro and re duration histograms were fit to shape parameters >1, demonstrating that the putative transition sub-states were also statistically deterministic. To identify how these SLE sub-intervals compared to the interval between the end of the termination sub-state and the beginning of the onset transition sub-state (the inter-SLE interval), we repeated the γ fit distribution on the histogram of the durations of the inter-SLE interval. We observed the inter-SLE interval had a γ fit duration of 0.76 [0.54 1.1] (n = 46 inter-SLE intervals, 17 slices, 14 subjects). This suggested that the inter-SLE interval was random, whereas the sub-intervals of the SLE state S2 were statistically deterministic.
To confirm the relationship between our classified SLE state and that of an amplitude-based identification of onset and termination, we applied ROC based analysis (Fig. 3D). The threshold for defining the SLE state in S2 was varied to obtain sensitivity and specificity measurements for the ROC curve. We observed that our classification was not 100% perfect, but close to 1, as indicated by the AUC of the ROC < 1, and by using other tests of sensitivity and specificity (Table 4). On inspection of the traces we observed that field bursting activity that traditionally defines the onset and termination of seizures is sometimes observed before, or after, the classified onset and termination transition sub-states.
One possibility is that the CFC features that were used to classify the SLE sub-states were primarily associated with the wave form shape, and not due to the existence of true low to high-frequency CFC (Kramer et al., 2008; Aru et al., 2015; Jensen et al., 2016). To address this, we tested the HMM using reconstructed data that was obtained by the lower envelope of the original dataset (Fig. 3). This removed the high-powered HFOs while preserving the power of the LFOs (Fig. 3C). Then, we re-computed the CFC features from this dataset, and re-tested the HMM. To account for variability in wave form shape, data were reconstructed using 2 forms of interpolation, linear and cubic. Both linear and cubic interpolated datasets performed poorly when testing on the original low-Mg2+ trained HMM (Fig. 3D), as demonstrated by low ROC characteristics as compared to the original dataset (Table 4). These data suggest that the low-frequency wave form shape is not sufficient to classify the CFC-based of the SLE state.
To further identify the impact of spike shape on the transition classification, data were bandpass filtered at 30–250 Hz, and large amplitude spikes were smoothly excised using a spike melting algorithm (Fig. 4A,a1,a2; Zalay et al., 2009). The amplitude of this signal and the phase of the low-frequency signal was then used to generate a feature set. We trained each feature set on the original feature set, then tested it on the new feature set with spikes omitted (Fig. 4Aa3). We asked (1) how much does the number of excised spikes alter the MPD of the HMM? (2) what is the relative difference between the defined transition state durations with and without the large spiking transients? and (3) do these differences have an effect on the γ-distribution shape parameter fits of the duration histograms? We observed that there was a significant positive moderate correlation between the relative number of excised spikes and the overall difference in the MPDs of S2 (Fig. 4B; Table 1, line 7). This effect led to, on average, a 2-s difference between the classified sub-state durations, and no significant differences between the sub-states (Fig. 4C). Furthermore, the shape parameters of the γ-distribution fits of the SLE duration sub-states were >1 (Fig. 4D), as previously observed. These data demonstrate that the high amplitude transients impact the HMM classifier within a narrow margin of error, without altering the statistical distributions of the SLE sub-states.
Large spikes do not affect the statistical distribution of SLE sub-states. A, a1, Example of a low-Mg2+ SLE. Inset is the expanded trace indicated by the box. It shows the occurrence of a large spike transient in the dataset. a2, Example that has been filtered at 30–250 Hz. Black trace indicates original filtered signal, blue trace indicates filtered signal with minimized large spike transient. Inset is the expanded trace at time indicated by the box in a1. a3, MPD of classified state S2. Black is the original dataset and blue is the dataset with large spikes minimized. Insets are classified onset and termination transitions; dt is the difference between the classified transition sub-states. Boxes used to identify the onset and termination sub-states. B, Correlation between the number of spikes removed from the datasets and the difference between the MPD of original S2 and the S2 with large spikes minimized (r2 = 0.1435, m = 1.37; b = 3.37; for Pearson correlation coefficient, see Table 1 line 7). C, Boxplot showing the difference between the classified sub-state transitions as defined by the MPD of the original S2 and S2 with large spikes minimized. The whiskers are the data points that are not outliers. The outliers are plotted using the ‘+’ symbol. D, Duration histograms of the onset, mid and termination intervals with the following γ-distribution shape parameter fits: 1.5 [1.0 2.4], mid 1.6 [1.0 2.5], termination 1.6 [1.0 2.6] (n = 30 SLEs, 13 slices, seven subjects).
4-Hz rhythm underlies statistical determinism of SLE transition sub-states
To identify what range of frequencies contained information about the deterministic properties of the SLE sub-intervals, we varied the frequency ranges of interest that were used to classify the sub-intervals of the SLE. Then we re-tested and re-trained the two-state HMM and re-computed the shape parameter fits of the γ-distributions obtained on their histograms (Table 3). As a representative example, we removed the 3- to 6-Hz frequency range from the original HMM feature set, and the model was re-trained and re-tested (Fig. 5). In this case, the onset and termination sub-states were Poisson distributed and random, respectively, as indicated by shape parameters ≤1, whereas the rs state was statistically deterministic. Furthermore, we observed that the onset sub-state remained, on average, similar in duration to the original trained model. Yet, the termination sub-state was markedly reduced in the majority of tests. This was in parallel to a prolonged rs state duration.
Shape parameters for other combinations of training for low-Mg2+ SLEs
The statistical distributions of the SLE transition sub-states are altered on omission of 3–6 Hz from the CFC feature set. A, Example of LFP (top) with MPDs of S2 tested using two feature sets. S2, The original feature set including 1–12 and 30–250 Hz. S2 (no 3–6 Hz), The feature set composed of 1–2.5 and 7–12 Hz. ro’, rs’, and re’ are the updated durations of the SLE sub-states as defined by the MPD of S2 on omission of 3–6 Hz. B, Histograms of the durations of ro’, rs’, and re’ for all SLEs with their respective γ-distribution fits. For a summary of the shape parameters of the γ-distribution fits, see Table 2.
Of 395 combinations of frequencies tested on our training set, 128 were successfully able to classify the SLE state, at varying degrees of accuracy. Training our dataset on a subset of ten cases, we preserved the statistical determinism of the middle rs state. In all but one case where the 4 Hz was included in the model, the SLE transition sub-states had shape parameters >1. When this frequency band was not included in the model, the SLE transition sub-states had shape parameters ≤∼1. This suggested that information from the 4-Hz range was necessary to classify statistically deterministic durations of the SLE termination transition sub-state.
Dynamics of seizure sub-states in humans can be distinguished using their underlying CFC features
To extend these findings to human seizures, we classified seizures that were recorded from patients before undergoing epilepsy surgery (Table 5). Because of interspecies variability, we re-trained the model to a human seizure trace, and then tested the remainder of the seizures obtained across these patients (Fig. 6A). We asked whether similar dynamical changes would result if we tested this model using the feature set with which we omitted the 3- to 6-Hz frequency range, as compared to that which included this range (Fig. 6B,C). When the SLE states were classified using the whole ICFC feature set, the SLE and termination sub-states had shape parameter fits >1, whereas the SLE onset sub-state had a shape parameter < 1 (Fig. 6B; Table 2). When the 3- to 6-Hz range was omitted from the feature set, there was a marked difference in the shape parameters of termination distribution histogram (Fig. 6C; Table 2). The SLE termination transition sub-state duration shifted to a shape parameter fit < 1, suggesting information contained within the 3- to 6-Hz frequency range was correlated with the statistical determinism of the termination state.
The statistical distribution of the human seizure termination sub-state is altered on omission of 3–6 Hz from the CFC feature set. A, a1, a2, Two different examples of iEEG traces of human seizures (a1, patient 4; a2, patient 2; top) with MPDs of S2 (SZ) tested using two different feature sets. Top, The original feature set including 1–12 and 30–250 Hz. Bottom, The feature set composed of 1–2 and 7–12 Hz. B, Histograms of the durations of ro, rs, and re with their respective γ-distribution fits. C, Histograms of the durations of ro’, rs’, and re’ with their respective γ-distribution fits. For a summary of the shape parameters of the γ-distribution fits, see Table 2 (n = 19 seizures, seven patients).
Testing in 4-AP-induced SLEs
To further generalize the results, we classified SLEs that were induced using bath application of 4-AP, 100 μM. Application of 4-AP in neocortical slice preparations resulted in spontaneous SLEs that began 5.2 ± 0.8 min after the addition of the 4-AP solution, and then recurred at a rate of 0.42 ± 0.04 per minute (for Ns, see Table 2). We asked whether the low-Mg2+ model was sufficient to train the network for classifying the 4-AP SLEs, the two-state HMM was trained using a feature set from the low-Mg2+ model, then tested on 4-AP (Fig. 7). The threshold for defining the SLE state in S2 was varied to obtain sensitivity and specificity measurements for the ROC curve. The 4-AP dataset was classified with high specificity and sensitivity (Fig. 7B; Table 4). Hence, the two-state HMM trained on low-Mg2+ SLEs classified the 4-AP-induced SLE state.
The statistical distribution of the 4-AP SLE termination sub-state is altered on omission of 3–6 Hz from the CFC features using low-Mg2+ SLEs for training the model. A, Representative trace of a SLE induced on addition of 4-AP, 100 µM. Middle trace is the state S2 of the MPD of the HMM from the LFP trace using 1–12 and 30–250 Hz. Bottom trace is the state S2 of the MPD of the HMM from the LFP trace using 1–2.5 and 7–12 Hz features for model training. ro’, rs’, and re’ are the updated durations of the SLE sub-states as defined by the MPD of S2 on omission of 3–6 Hz. B, Receiving operating characteristics curves applied to both in vitro datasets as a result of training the HMM on one low-Mg2+ SLE trace and testing on both low-Mg2+ and 4-AP SLEs. Inset shows expanded exponential fits of a segment of the ROC curve. Refer to Table 2 for a complete characterization of the ROC curve parameters (low-Mg2+ SLEs: n = 46 SLEs, 14 slices, 10 subjects; 4-AP SLEs: n = 28 SLEs, seven slices, seven subjects). C, Histograms of the durations of ro, rs, and re with their respective γ-distribution fits. D, Histograms of the durations of ro’, rs’, and re’ with their respective γ-distribution fits. See Table 2 for shape parameters of the γ-distribution fits and duration averages for each interval.
Summary of statistics for ROC curves
iEEG patient details
We asked whether similar dynamical changes would result if we tested this model using the feature set with which we omitted the 3- to 6-Hz frequency range, encompassing the 4-Hz rhythm not necessary. When the SLE states were classified using the whole ICFC feature set, the SLE intervals had a shape parameter fit >1 (Fig. 7C). When the 3- to 6-Hz range was omitted from the feature set, there was no difference in the duration the of SLE intervals as compared to the first feature set used (Table 2). However, there was a marked difference in the shape parameters of duration distributions histograms (Fig. 7D). The SLE termination transition sub-state duration had a shape parameter fit <1, suggesting information contained within this frequency range was correlated with the statistical distribution of the termination state. Yet, two distinctive peaks emerged in the rs and re state, as indicated by the arrows in Figure 7D. This suggested that there were two separate dynamical activities present in this 4-AP model.
We theorized that the presence of two underlying dynamical distributions suggested a time-dependent change in the 4-AP SLE characteristics that was not present in the low-Mg2+-induced SLEs (Fig. 8). We grouped the first four SLEs along a given time interval within a single slice as numbers 1 through 4 for both datasets. For low-Mg2+, we observed that the re’ sub-state was of short duration, and remained of short duration, over time (Fig. 8C, left). In the case of 4-AP, the re’ sub-state duration was short for the first SLE relative to the 4th SLE (Fig. 8C, right). These data suggested that the temporal dynamics of the 4-AP SLE model changed over the course of an experimental protocol, and, at least initially, a similar dynamical mechanism that contained information from the 4-Hz rhythm existed between both models.
Temporal relationship of the SLE termination sub-state. A, Representative trace of the low-Mg2+ SLE model with its respective MPD of the HMM. Boxes indicate SLE termination transition classification. B, Representative trace of the 4-AP SLE model with its respective MPD of the HMM. Boxes indicate the SLE termination transition classification. C, Boxplots showing the average duration of the termination transition of 4-AP and low-Mg2+ SLEs (*p < 0.05; 4-AP n = 24 SLEs, six slices, six subjects; low-Mg2+ n = 28 SLEs, seven slices, seven subjects; see Table 1, lines 8–15). The whiskers are the data points that are not outliers. The outliers are plotted using the ‘+’ symbol. For 4-AP, one subject was omitted from this calculation as only the first two SLE transition sub-states were classified accurately.
Identifying the 4-Hz HFO-coupled oscillation during the transition states in low-Mg2+ SLEs
Given that the dynamics of the SLE sub-state transitions were dependent on the 4-Hz rhythm, we next wanted to confirm that this rhythm is a valid phase-amplitude CFC signal (Fig. 9). First, we inspected the unfiltered signal and observed that high-frequency activity of high amplitude occurred transiently on field bursting activity during the transition periods, some of which was oscillating at 4 Hz. Then, we applied surrogate analysis to distinguish the range of high-frequency signals that were significantly correlated with the 4-Hz rhythm. During a 5-s window at the beginning of the onset, and during a 5-s window at the beginning of termination, this range was variable, peaking between 30 and 100 Hz. Using this significant range for each SLE, we computed the correlation between the ICFC and the power of the 4-Hz rhythm within the full transition states using 2-s windows to compute the CFC indices with 1-s overlap. If there was significant coupling, then this correlation would be either negatively correlated or not correlated with the 4-Hz signal power (Jensen et al., 2016). As a control, we also identified the correlation between the 4-Hz signal power and the average CFC indices from the 200- to 250-Hz range, which was not significantly coupled to the 4-Hz rhythm during the first 5 s of the onset transition and the first 5 s of the termination transition. For the 43 SLEs tested, 91% had a range of significant 4-Hz rhythm within the onset state, whereas at the termination state, 72% were significant. This is in comparison to the 200- to 250-Hz range with minimal significance during the onset and termination states (Table 1, lines 16–19). These data suggest that there is significant coupling between the 4-Hz oscillation at specific ranges of high-frequency rhythms, both at the onset and termination of the SLE state. During the termination sub-state, there is minimal coupling between the 4-Hz rhythm and the 200- to 250-Hz oscillation, as indicated by a p > 0.05. Because the average p value for the Pearson correlation coefficient during the onset transition was >0.05, likely the 4-Hz rhythm couples transiently to the 200- to 250-Hz range closer to the initiation of rs.
CFC coupled 4-Hz rhythm present during onset and termination phases of the SLE state. A–C, Three examples LFP recordings of low-Mg2+-induced SLEs. Boxes indicate 5-s windows at which the CFC features were tested for significance. Inset in B is expanded trace from LFP segment showing HFOs present on LFOs. a1–c1, CFC indices for the 5-s windows indicated by boxes. Regions of significance designated by black contour lines. a2, Expanded traces from A demonstrating HFOs at specific phases of LFOs. b2, c2, Examples of correlation fits of the ICFC and the cwt power of the 4-Hz oscillation during the ro sub-states (left) and re sub-states (right) for representative examples (p values indicated for Pearson correlation coefficients; for average correlation coefficients, see Table 1, lines 16–19 n = 43 SLE, 14 slices, 10 subjects).
Excitatory and inhibitory current properties during the SLE transition sub-states
To track the temporal dynamics of excitatory and inhibitory currents during the SLE transition sub-states, pyramidal neurons from Layer II/III were whole-cell patched and voltage clamped at –70 mV (reversal potential for IPSCs) and 0 mV (reversal potential for EPSCs; Figs. 10, 11). Then, we computed the charge transfer energies for the SLE transition sub-states as previously defined using the HMM. Because each SLE had a different transition time, data were binned into 10 windows of analysis, representing a percentage of the full transition state. The onset transition sub-state also included 10% before the onset, which was used for significance, and 10% after the transition, which defined the beginning of the SLE state rs.
Excitatory currents (sEPSCs) in the onset (ro) and end (re) sub-states of the SLE. A, LFP (top) and intracellular voltage clamp trace (bottom) of a low-Mg2+-induced SLE. The SLE transition sub-states indicated by the boxes. The cells are held at –70 mV for sEPSC acquisition. The data were binned into 10 equal segments that on a scale of 0–100 corresponding to a percentage of the SLE transition sub-states. B, Expanded traces of onset and termination sub-states from the representative example in A. C, Average charge transfer energy of the sEPSCs with horizontal lines at the top showing significance relative to pre-SLE baseline. The whiskers are the data points that are not outliers. The outliers are plotted using the ‘+’ symbol. Significance calculated relative to the 10% before the onset sub-state (*p < 0.05; 29 SLEs, 10 slices, nine subjects; Table 1, lines 20–41).
Inhibitory currents (sIPSCs) in the onset (ro) and end (re) sub-states of the SLE. A, LFP (top) and intracellular voltage clamp trace of a low-Mg2+-induced SLE. The SLE transition sub-states indicated by the boxes. The cell is held at 0 mV for sIPSCs acquisition. The data were binned into 10 equal segments that on a scale of 0–100 corresponding to a percentage of the SLE transition sub-state. B, Expanded traces of onset and termination sub-states from representative example in A. C, Average charge transfer energy of the sIPSCs (*p < 0.05; 19 SLEs, six slices, five subjects, see Table 1, lines 42–63). The whiskers are the data points that are not outliers. The outliers are plotted using the ‘+’ symbol. Significance calculated relative to the 10% before the onset transition sub-state.
During the onset sub-state there was a large rise in the charge transfer energy of the sEPSCs. During the termination sub-state these currents remained elevated, then recovered to baseline levels at 80–90% of the termination sub-state (Fig. 10). In contrast, for the sIPSCs, there was a temporary rise in sIPSCs, which peaked, then dropped, at varying temporal windows during the transition interval (Fig. 11). Furthermore, the sIPSCs were at pre-SLE baseline levels during the SLE termination sub-state except for the initial 10–30% of the transition. Hence, large sEPSCs mark both the onset and termination SLE transitions, whereas sIPSCs dynamically peaked at varying intervals during both sub-states.
Peak phase coherence between spontaneous excitatory and inhibitory currents and the low-frequency rhythm
Next, we identified a population of currents that were coherent with the field rhythms (Figs. 12,13). We used a threshold measurement to identify peaks in sEPSCs and sIPSCs, and then varied the threshold until a PLV > 0.3 was observed (between 1 and 12 Hz). The threshold for detecting phase coherent currents (PCCs) was different for each SLE. For the sEPSCs (Fig. 12), during the onset and termination sub-states, 22/29 of the SLEs showed intermediate to high PLVs peaking at varying frequencies between 1 and 12 Hz for both onset and termination. In both sub-states, ∼70% of the cases were within the intermediate range, and 30% in the high range of PLVs We also used a Rayleigh test to test the cases for significance. In this case, during the onset, only 31% of the SLEs passed, which was consistent with those that had high PLVs Whereas, during the termination, 62% passed this test. Furthermore, there was no clear phase advancement or phase delay from onset to termination of the SLEs at the 4-Hz rhythm of interest. We hypothesized that the non-significant currents with moderate to high PLVs may, nonetheless, be integral to the dynamics of the transition sub-states.
Phase coherent excitatory currents (PCCs) during the SLE transition sub-states. A, LFP (top) and intracellular voltage clamp trace of a low-Mg2+-induced SLE. The SLE transition sub-states are indicated by the boxes. The cell is held at –70 mV for sEPSC acquisition. B, Peak detection of the sEPSCs and the filtered LFP at 1–3 Hz during the onset (ro) and termination (re) transition sub-states. Gray circles indicate detected sEPSCs. These are superimposed at the phase of the LFP oscillation at the time of detection. Thick black lines show average of the timing of the sEPSCs relative to the phase of the oscillation. C, Circular histograms of phase at which peak of sEPSC was detected from example in B demonstrating PLVs during the onset and termination phases indicated by the length of the black direction arrow. D, Characterization of the PCCs as compared to baseline events (baseline). Baseline was chosen as 10 s just before the onset sub-state transition. The average amplitude of the PCC is 48 ± 2 pA (mean ± SEM), whereas the average amplitude of the baseline currents is 16 ± 2 pA (*p = 1.1E-94; 7402 baseline events, 1544 PCC events, 29 SLEs, 10 slices, nine subjects; see Table 1, lines 64–65). E, Boxplots showing the median frequencies of sEPSCs with the 25th and 75th percentile ranges (*p < 0.001; see Table 1, lines 66–67). The whiskers are the data points that are not outliers. F, Histogram of the frequency of maximal PLVs for the events showing moderate to high phase locking.
Phase coherent inhibitory currents during the SLE transition sub-states. A, LFP (top) and intracellular voltage clamp trace of a low-Mg2+-induced SLE. The SLE state transition sub-states indicated by the boxes. The cell is held at 0 mV for sIPSC acquisition. B, Peak detection of the sIPSCs and the filtered LFP at 1–3 Hz during the onset (ro) and termination (re) state transition sub-states. Gray circles indicate detected sIPSCs. These are superimposed at the phase of the LFP oscillation at the time of detection. Thick black lines show the pattern of the timing of the sEPSCs relative to the phase of the oscillation. C, Circular histograms of phase at which peak of sEPSC was detected from example in B demonstrating PLV during onset and termination phases indicated by the length of the black direction arrow. Boxplot shows the frequency of maximal PLVs during SLE transition sub-states with the 25th and 75th percentile ranges (n = 19 SLEs, six slices, five subjects). D, Average sIPSCs detected during the onset and termination sub-states. Baseline was chosen as 10 s just before the onset sub-state transition. The average amplitude is 128 ± 3 pA (mean ± SEM), the average amplitude of the baseline events is 15 ± 4 pA (p = 2.2E-98; 3155 baseline events, 606 PCC events, 19 SLEs, six slices, five subjects; Table 1, lines 68–67). E, Boxplots showing the median frequencies of sIPSCs with the 25th and 75th percentile ranges (*p < 0.001; Table 1, lines 70–71). The whiskers are the data points that are not outliers. The outliers are plotted using the ‘+’ symbol. F, Histogram of the maximal PLVs for the events showing moderate to high phase locking.
We identified that the amplitude of these events was higher than those at baseline (Fig. 12D). Furthermore, the frequency of the currents at baseline was 13 Hz. This was also the same for the average frequency of all the spontaneous events during the transition sub-states. Yet, the PCCs had a median frequency of 2.3 Hz (Fig. 12E), and were spread across multiple frequency bands. This was further indicated by the histogram of maximal PLVs for all of the events showing moderate to high phase locking, which peaks at 1–2 and 9–10 Hz (Fig. 12F).
These findings were in contrast with the sIPSCs (Fig. 13). For both transition sub-states, 15/19 SLEs showed intermediate to high PLVs. From this, 63% of these were within the intermediate range, whereas 36% were within the high range. In contrast, during the termination, 54% of SLEs showed intermediate PLVs (0.36–0.44). The remainder were of low PLVs (<0.3). These were also phase coherent to a range of frequencies between 1 and 12 Hz. Using a Rayleigh test of significance, 73% of sIPSCs passed in the SLE onset, and 40% during the termination. Furthermore, the subset of common high PLVs (54%) had a phase advancement of 70 ± 11.2° phase angle difference from the onset to the termination. These results were not correlated with slice or subject number.
As with the sEPSCs, the amplitude of these events was higher than those at baseline (Fig. 13D). Furthermore, the frequency of the currents at baseline was ∼14 Hz. As with the sEPSCs, this was the same for the average frequency of all the spontaneous events during the transition sub-states. The sIPSC PCCs had an average frequency of 3.7 Hz (Fig. 13E), consistent with a correlation between the sIPSCs and the 4-Hz rhythm. This was further indicated by the histogram of maximal PLVs for all of the events showing moderate to high phase locking, which peaks at 1–4 and 11–12 Hz (Fig. 13F).
Specificity of PCCs to SLE transition sub-states
To determine whether the PCCs were specific to the transition sub-states, we performed similar tests for identifying PCCs on the period 10 s before the SLE onset sub-state, and on the SLE state rs. During the inter-SLE period, only 12% sEPSC-SLE pairings showed significant and high PLVs, with a total of 5 detected events, and only 9% of the sIPSC-SLE pairings showed moderate PLVs with 19 total detected events. These results were similar during the ictal state, rs, where only 19% sEPSC-SLE pairings showed significantly moderate PLVs with a total of 35 detected events, whereas only 18% sIPSC-SLE pairings showed moderate PLVs with a total of 63 detected events. This suggested that the PCCs were specific to the onset and termination sub-state.
sIPSCs show phase locking to the 4- to 5-Hz oscillation
To quantify the coherence between the intracellular currents and the low-frequency field oscillations, we calculated the frequency dependent correlation between the PLVs and the excitatory and inhibitory currents over the low-frequency range (Fig. 14). At the threshold for defining the PCCs, we examined the average PLVs at varying frequency bands over both onset and termination. For sEPSCs, we observed phase locking maximal at 2 Hz, consistent with the average frequency of the PCCs. For sIPSCs, we observed an elevated 4- to 5-Hz PLV, consistent with the average frequency of the sIPSCs. At 5 Hz, the median PLV breached the moderate phase locking threshold, suggesting sIPSCs, on average, moderately phase locked to 5 Hz. We then compared the PLVs between sEPSCs and sIPSCs (Table 1, lines 88–99). In this case, we observed elevated phase locking for both 4 and 5 Hz for sIPSCs as compared to sEPSCs. This suggests stronger phase locking between sIPSCs and the 4- to 5-Hz frequency range as compared to sEPSCs.
PLVs for sEPSCs and sIPSCs. Boxplot showing the median PLVs for sEPSC PCCs and sIPSC PCCs with their 25th and 75th percentiles. Whiskers indicate extreme data points that are not considered outliers. Outliers are indicated by the ‘+’ symbol. (*p < 0.05; Table 1, lines 72–99; sEPSCs n = 29 SLEs, 10 slices, nine subjects; sIPSCs n = 19 SLEs, six slices, five subjects).
Discussion
To the best of our knowledge, we are the first to demonstrate how specific CFC features demark statistical determinism of the SLE sub-states. We ask the following questions: what states do CFC features classify? Can these states be predicted? What is the underlying excitatory and inhibitory current bases for the predictability of the classified states?
CFC classified SLE sub-states
We observe that CFC between the phase of 1–12 Hz and the amplitude of 30–250 Hz classify SLE state, and SLE onset and termination transition sub-states. A positive correlation between the SLE onset sub-state, and the SLE state, and SLE duration/intensity, suggests these intervals are distinct sub-states of the SLE. The termination state, which is comprised of bursting discharges, only positively correlates with the SLE duration. Field bursting activity that traditionally defines the onset and termination of seizures (Dreier and Heinemann, 1991; Boido et al., 2014; Ellender et al., 2014), is sometimes observed before the classified onset sub-state and after the classified termination sub-state. Field bursting activity is composed of low-frequency and high-amplitude events; however, to make the feature set, we used only information from the phase of the LFOs. An HMM is known to capture essential patterns within features, with the assumption that each pattern performs a distinct role in the network (Rabiner and Juang, 1986); hence, the output is highly dependent on our choice of features. Using the CFC features, with only low-frequency phase information to run the HMM, classifies sub-intervals of the SLEs that were not previously observed.
Can seizure sub-states be predicted?
The CFC chosen for the HMM classifies three distinct SLE intervals that are statistically deterministic, in vitro. The histogram of the duration of the three SLE sub-intervals had a shape parameter >1, suggesting that both the onset and termination of SLEs are statistically deterministic. This is in line with previous results demonstrating that a deterministic process underlies the seizure termination state; hence, it can be predicted (Bauer et al., 2017). This is also consistent across both low-Mg2+ and 4-AP SLE models, and in human iEEG seizures. Furthermore, like with low-Mg2+-induced hippocampal SLEs (Suffczynski et al., 2006), the SLE onset is statistically deterministic, as indicated by a shape parameter >1. Yet, the previous work suggested that this is an exception in the low-Mg2+ model, as the interseizure intervals from human seizures have a shape parameter <1, and therefore are randomly occurring and unpredictable. These previous results computed the duration of the full interseizure interval from both human and animal models. We also obtained a shape parameter <1 for the human iEEG onset sub-state, identifying this sub-state as random in humans. As our results focused on the seizure state and its respective sub-states, future studies may examine interictal human iEEG durations to determine if a predictable interictal period will emerge on removal of the unpredictable seizure onset sub-state.
The random occurrence of the SLE state is also observable in our model, in vitro, if we remove the transition sub-states from the computation of the inter-SLE interval. These results are also consistent with the analysis done in in vivo data (Colic et al., 2013), and further confirmed in a computational model (Grigorovsky and Bardakjian, 2018). Hence, the beginning of the onset transition sub-state may emerge randomly; however, the transition into the SLE state, rs emerges predictably, in vitro. Therefore, CFC features classify predictable sub-intervals of the SLEs that were not previously observed.
We consider alternatives when making the assertion of predictability. First, a γ-distribution shape parameter >1 can also be the addition of multiple random variables that are exponentially distributed with the same averages (Cox and Lewis, 1978). Yet, the consistency across models and species, suggests that underlying dynamical mechanisms are at work. Second, we assume that statistical determinism implies predictability. Third, this predictability could be an artifact of the HMM, as it has previously been demonstrated that the HMM can predict human seizure onset using the iEEG signal directly (Wong et al., 2007). Yet, the transition to seizures can be detected in advance of the clinical onset using CFC and a different classification approach, suggesting the seizure onset is predictable and this is not an artifact of the HMM (Jacobs et al., 2018). Changes in CFC features have been demonstrated as biomarkers of hyperexcitable brain states (Guirgis et al., 2015; Weiss et al., 2016; Grigorovsky and Bardakjian, 2018; Samiee et al., 2018). Therefore, we have demonstrated that seizure transition probability can be monitored by changes in the underlying network excitability (Kalitzin et al., 2005). As such, we have found a novel SLE interval, in vitro, that has predictable underlying dynamics independent of the inter-SLE interval.
The 4-Hz phase underlies SLE transition sub-state predictability
Statistical determinism of the low-Mg2+ SLE onset and termination sub-states is observable with information contained in the 4 Hz frequency band. We observe a 4 Hz signal which, when omitted, results in ∼1 (confidence interval (CI): [0.81 1.9]) for the onset, and <1 (CI: [0.5 1]) for the termination. Although a slightly larger interval that contains the 3 Hz leads to random underlying dynamics <1 Hz (CI: [0.4 0.8]) for seizure termination, the striking shift to a near-Poisson distribution on omission of the 4 Hz is of interest, as it provides information about a network that shifts from predictable to unpredictable.
Previous work has demonstrated that there are multiple sub-states within the SLE and during pre-SLE state (Guirgis et al., 2014). The sub-states of the SLE were previously classified using 4–8 Hz, which contained sufficient network information to classify three sub-states of the SLE state in the hippocampus, which is associated with physiologic θ-γ coupling in the neocortex (Canolty et al., 2006). The transition sub-state histogram distributions are deterministic only when information contained within the 4-Hz range is added to the HMM. Hence, we need only this range to capture the information of statistical determinism of the transition sub-states, overlapping with these previous results. In contrast, for the 4-AP SLE model, omitting this range only affects the mid-SLE state and the termination transition sub-state in 4-AP, rather than the onset state. Furthermore, omitting 3- to 6-Hz range from human SLEs only effects the termination sub-state. This suggests the underlying dynamical mechanisms of SLE onset in the 4-AP model, and in human seizures, are not dependent on 4 Hz. Hence, the 4-Hz range most likely underlies the statistical determinism of seizure termination across models.
This is also confirmed when we observed that the removal of a range which contains the 4-Hz oscillation results in a large reduction in the duration of the classified termination sub-state, with only minor differences in the onset sub-state duration, yet a high proportion of SLEs have PCCs that are phase locked to this rhythm during the onset. We ask why the 4-Hz rhythm only has a major effect on the duration of the termination. It has previously been observed that a 400- to 800-Hz signal emerges during the onset transition of low-Mg2+ SLEs in the hippocampus (Lasztóczi et al., 2004). These findings were further confirmed in a similar SLE model (Khosravani et al., 2005). Also, in the cortex, fast ripples >250 Hz have been implicated in the initiation of seizures (Bragin et al., 1999; Jacobs et al., 2008). Therefore, this 4-Hz oscillation may have effects on the onset sub-state, but coupled with a higher frequency range that was not added to the model features; hence, not sensitive in the classifier.
This higher frequency activity in the fast ripple range is likely a summation of pyramidal cell burst spiking in the hippocampus (Dzhala and Staley, 2004). The neocortex may have a similar underlying mechanism. This is observed in our model as a large rise in charge transfer energy of sEPSCs during the onset phase. But these are coherent with a range of frequencies 1–12 Hz within the onset sub-state, clustering at 1–2 and 9–10 Hz. This suggests that a large proportion of the onset duration transition is correlated with the activity of pyramidal neuron spiking, and only a small proportion shifts the underlying dynamics of the onset.
Excitatory and inhibitory current networks in SLE transition sub-states
The correlation between synaptic currents and the onset of low-Mg2+ SLE model has been identified previously (Traub et al., 1994). The faster onset time that is a possibility is unlikely related to the CFC based transitions observed because we observe a consistent slow rise of excitatory currents during the classified transition sub-state, and slow rise, then fall, of inhibitory currents. Furthermore, this pattern parallels multiple previous reports in vitro and in vivo (Kawaguchi, 2001; Timofeev and Steriade, 2004; Trevelyan et al., 2006). However, the traditionally defined onset that consists of the largest burst of inhibition synchronous with excitation previously observed (Trevelyan et al., 2006), did not consistently correlate with rs SLE onset. Furthermore, the PCCs we observe that are specific to these transition sub-states, and are of larger amplitude than baseline, are related to bursts of inhibition and excitation temporarily coherent with the field activity. This is also confirmed by the observation that these PCCs are not following the average spontaneous event frequency during the transition sub-states and pre-SLE baseline (∼13 Hz) but are following a much lower frequency (2–4 Hz).
Phase locking of the GABAergic network to the 4-Hz oscillation
Because of the higher phase locking of the sIPSCs around 4 Hz, and scattered phase locking of the sEPSCs across multiple frequencies 1–12 Hz, we observed significantly elevated PLVs of sIPSCs compared to sEPSCs in the 4-Hz frequency bands. Hence, we suggest that the sIPSCs make a stronger contribution to the 4-Hz rhythm than the sEPSCs. It is also in line with previous reports of phase locking of action potentials during the SLE onset transition in the low-Mg2+ model (Quilichini et al., 2012). Our results validate these previous findings in this older animal model in the neocortex. Differences in the overall phase locking strength could also be the timing and amount of excitatory and inhibitory connections received by our recorded pyramidal neurons (Huh et al., 2016). Also, it could be the number of detected events. Yet, we only observed PLVs significant and above 0.3 within distinct frequency bands at onset and termination, and we used a large number of events for each case. Therefore, our results are likely due to the underlying network activity rather than the number of events.
Given that we did not block sodium conductance it is likely that these PCCs are, in part, neuronal action potentials (Buzsáki et al., 2012). Our examination of currents rather than action potentials limited the peak detection algorithm to the refractoriness of the currents, which we set at a maximum of 20 Hz. This suggests that although there may be higher frequencies of PCCs (in the range of HFO amplitude detection) we cannot test this using our current approach. Also, given we used a threshold method for detection of peak currents, there may be PCCs during the SLE state that are not distinguishable using our approach (Ziburkus et al., 2006; Cotic et al., 2011). Nonetheless, another lab has observed that HFOs maximally coupling to the 4-Hz oscillation reflect the summation of IPSPs in pyramidal neurons (Jefferys et al., 2012). The summative effects may be why the PCCs we observe are of much higher amplitude than the baseline events. Another lab has observed phase locking of action potentials occurs during multiple phases of the SLE state (Lévesque et al., 2016). Given the specificity of the PCCs to the SLE transition sub-states, action potential spiking during the various phases of the SLE state does not reflect the postsynaptic neuronal response.
It has been previously observed that interneurons which fired preferentially at θ (4–10 Hz) increase their frequency during the transition to the ictal event (Karunakaran et al., 2016); whereas there is a concomitant slowing (or static activation) of pyramidal neuron action potentials (Grasse et al., 2013; Lévesque et al., 2016). Yet, we observed that the postsynaptic neuron is receiving a static population of both high amplitude excitatory and inhibitory currents during the transition sub-states. This suggests that a non-action potential dependent mechanism synchronizes the inputs from the population of excitatory and inhibitory neurons during the transitions of SLEs.
Unlike previous results that focused on either SLE onset or termination, we report the presence of PCCs during multiple phases of the SLE, therefore can compare the phase relationships between the state transitions. Specifically, when focusing on the 4-Hz frequency to calculate instantaneous phase at which the PCCs occur, we observe a phase advancement for the sIPSC PCCs in the termination sub-state relative to the onset sub-state, whereas this observation was not consistent with the excitatory currents. This further validates our hypothesis that sIPSCs play a dominant role in the 4-Hz network. Others have demonstrated that inhibitory cells tend to fire action potentials during the trough of the 5- to 15-Hz oscillation, suggesting the currents may exhibit a particular phase preference (Lévesque et al., 2016). Yet, we did not observe a consistent phase relationship between SLEs during the onset or the termination of the SLEs. These previous results were reliant on observations made from bandpass filtered signals, whereas our instantaneous phase from the cwt, may have resulted in different conclusions. This may also be due to intrinsic cellular mechanisms affecting the timing of the detected peaks, such as dynamic changes in ionic currents, or cell swelling. Yet, it may also suggest that the inhibitory current population active at termination is different from the onset (Ellender et al., 2014; Wen et al., 2015).
This work supports the concept that GABAergic cells are able to synchronize local neuronal networks (Chagnac-Amitai and Connors, 1989; Cobb et al., 1995; Quilichini et al., 2012; Chang et al., 2018). We observed that inhibition predominantly phase locks to the 4-Hz oscillation. We observe this in both the frequency of the PCCs at ∼4 Hz, and the PLV of these currents significantly higher at the 4 Hz for inhibition relative to excitation. The presence of a 4-Hz oscillation has been previously observed as part of the θ rhythm in the neocortex (Siapas et al., 2005; Sirota et al., 2008; Fujisawa and Buzsáki, 2011). But, the 4 Hz on its own has been demonstrated as a task-specific rhythm phase locked with GABAergic networks (Fujisawa and Buzsáki, 2011). This previous study observed a larger fraction of interneurons than pyramidal neurons phase locked to the 4-Hz oscillation in the pre-frontal cortex. We observe a similar phase locking bias in the somatosensory cortex under epileptic conditions. In this θ range, interneurons have been observed to synchronize network oscillations via gap junctions (Galarreta and Hestrin, 1999; Gibson et al., 1999; Beierlein et al., 2000; Tamás et al., 2000; Blatow et al., 2003). Hence, the synergistic action of chemical synapses between interneurons and pyramidal neurons and electrical synapses between interneurons may not only synchronize the neuronal populations in the context of maintaining a physiologic 4-Hz rhythm; it may further play a role in shaping the dynamics of SLE sub-intervals.
Future directions
Our results uncover a population of large amplitude currents, both excitatory and inhibitory, that show peak phase-coherence with LFOs emerging during the transition sub-states of SLEs. Our observations lead to the following conclusions: There is higher phase locking of inhibition to the 4-Hz oscillation. The omission of the 4-Hz oscillation leads to Poisson distributed SLE transition sub-states, and the Poisson distribution follows that of the naturally occurring action potential firing frequency (Brunel, 2000). This suggests that inhibition maintains the predictability of the SLE transition sub-states. Future studies may consider targeting the coherent inhibitory events for direct modulation of the seizure termination phase. This work merges the underlying cellular mechanisms to fast and slow rhythms in the neural-glial network of the brain. In doing so, it can help unravel mechanism behind the activation of the local network during seizures and will guide future studies in the mechanisms underlying endogenous seizure transition dynamics.
Acknowledgments
Acknowledgements: We thank Dr. Martin Del Campo (Toronto Western Hospital) and Dr. Yotin Chinvarun (Bangkok hospital) for providing human iEEG seizure data. We also thank Dr. Liang Zhang for advice on the approach, Dr. Mirna Guirgis for providing the HMM MATLAB code, Dr. Lihua Wang for providing 4-AP data to test the model, and Dene Ringuette, Vasily Grigorovsky, and Daniel Jacobs for support in the formation of this manuscript.
Footnotes
The authors declare no competing financial interests.
This work was supported by the Canadian Institute of Health Research and Natural Sciences and Engineering Research Council of Canada.
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
References
In this issue
