Abstract

Juvenile myoclonic epilepsy (JME) is a form of idiopathic generalized epilepsy. It is yet unclear to what extent JME leads to abnormal network activation patterns. Here, we characterized statistical regularities in magnetoencephalograph (MEG) resting-state networks and their differences between JME patients and controls by combining a pairwise maximum entropy model (pMEM) and novel energy landscape analyses for MEG. First, we fitted the pMEM to the MEG oscillatory power in the front-oparietal network (FPN) and other resting-state networks, which provided a good estimation of the occurrence probability of network states. Then, we used energy values derived from the pMEM to depict an energy landscape, with a higher energy state corresponding to a lower occurrence probability. JME patients showed fewer local energy minima than controls and had elevated energy values for the FPN within the theta, beta, and gamma bands. Furthermore, simulations of the fitted pMEM showed that the proportion of time the FPN was occupied within the basins of energy minima was shortened in JME patients. These network alterations were highlighted by significant classification of individual participants employing energy values as multivariate features. Our findings suggested that JME patients had altered multistability in selective functional networks and frequency bands in the fronto-parietal cortices.

Author Summary

We proposed an energy landscape method to quantify the occurrence probability of network states in magnetoencephalograph (MEG) oscillatory power during rest, which was derived from a pairwise maximum entropy model (pMEM). We compared the energy landscapes measures of three resting-state networks between patients with juvenile myoclonic epilepsy (JME) and healthy controls. The pMEM provided a good fit to the binarized MEG oscillatory power in both patients and controls. Patients with JME exhibited fewer local minima of the energy and elevated energy values than controls, predominately in the fronto-parietal network across multiple frequency bands. Furthermore, multivariate features constructed from energy landscapes allowed significant single-patient classification. Our results further highlighted the pMEM as a descriptive, generative, and predictive model for characterizing atypical functional network properties in brain disorders.

INTRODUCTION

Juvenile myoclonic epilepsy (JME) is the most common syndrome of the wider group of idiopathic generalized epilepsies (Wolf & Beniczky, 2014). Patients with JME often exhibit three main types of seizures: myoclonic, absence, and generalized tonic-clonic seizures (Wolf et al., 2015a). Typical JME characteristics are normal or close to normal clinical MRI of the brain and interictal EEG with irregular spike-waves or polyspike-waves with frontal predominance (Camfield, Striano, & Camfield, 2013a). JME patients are susceptible to seizure precipitation after sleep deprivation, alcohol usage, excise, or demanding cognitive processing (Delgado-Escueta & Enrile-Bacsal, 1984; Yacubian & Wolf, 2014). JME is a lifelong condition, and treatment with antiepileptic drugs is usually necessary.

Although the pathogenetic mechanisms of JME are still not fully understood (Berkovic, Howell, Hay, & Hopper, 1998), JME has been recognized as a network disorder affecting brain activity and connectivity that leads to cognitive impairments (Chowdhury et al., 2014; Wolf et al., 2015a) and personality traits similar to patients with frontal lobe lesions (Engel, 1997). BOLD functional MRI (fMRI) and diffusion-weighted imaging showed hyperconnectivity in the frontal lobe in JME (Caeyenberghs et al., 2015; Vollmar et al., 2012). Electrophysiological data suggests that JME has an impact on multiple functional networks, including the fronto-parietal network (FPN) (Wolf et al., 2015a), the default mode network (DMN) (McGill et al., 2012), and the sensorimotor network (SMN) (Clemens et al., 2013), which may be driven by dysfunctional thalamocortical circuitry (Betting et al., 2006; Gotman et al., 2005; Hamandi et al., 2006; J. H. Kim et al., 2007).

Several sensitive markers from resting EEG and magnetoencephalography (MEG) recordings have been identified for classifying patients with epilepsy and predicting seizure onsets, including information entropy (Kannathal, Choo, Acharya, & Sadasivan, 2005; Song, Crowcroft, & Zhang, 2012; Song & Zhang, 2013), Lyapunov exponent (Babloyantz & Destexhe, 1986; Iasemidis, Chris Sackellares, Zaveri, & Williams, 1990), and phase plane portraits (Iasemidis et al., 1990). These methods describe statistical regularities of electrophysiological signals from a dynamical system perspective, in line with the theoretical account of epileptic seizures as bifurcations from stable states (da Silva et al., 2003). In JME, however, it is yet unclear whether atypical statistical properties of network activation are present during rest, and if so, whether the changes are frequency specific.

This study addressed these problems by applying a pairwise maximum entropy model (pMEM) approach (Yeh et al., 2010) to source-localized, frequency-specific MEG resting-sate oscillatory activity (Figure 1). The pMEM is a statistical model of the occurrence probability of network states, with its parameters being constrained by the network’s regional activity and pairwise regional coactivation from empirical data. According to the principle of maximum entropy, the pMEM is the most parsimonious second-order model of a system with minimum assumptions (Jaynes, 1957), and it permits multistability in a system with metastability states (Cirillo & Lebowitz, 1998; Deco, Senden, & Jirsa, 2012). The pMEM has been successfully applied to the collective behavior of spiking neural networks (Bialek, 2017; Schneidman, Berry, Segev, & Bialek, 2006; Tang et al., 2008; Tkacik, Schneidman, Berry, Michael, & Bialek, 2006) and BOLD fMRI responses (Ashourvan, Gu, Mattar, Vettel, & Bassett, 2017; Ezaki, Sakaki, Watanabe, & Masuda, 2018; Watanabe et al., 2013; Watanabe et al., 2014a). Here, we extended this theoretical framework to MEG oscillatory activity in three functional networks: FPN, DMN and SMN. Furthermore, based on the pMEM fitted to individual participants, we depicted an energy landscape for each of the networks at theta (4–7 Hz), alpha (8–13 Hz), beta (15–25 Hz), and gamma (30–60 Hz) bands. The energy landscape is a graphical representation of all network states and their energy values (Ezaki, Watanabe, Ohzeki, & Masuda, 2017). We then compared several quantitative measures obtained from the energy landscapes between JME patients and controls.

Figure 1. 

Illustration of the energy landscape analysis on a network of four regions of interst (ROIs). (A), Selection of ROIs from the source-space signals. (B), Signal filtering in frequency bands of interest. (C), Envelope extraction using the absolute value of the analytical representation of the signal. (D), Binarization of the data. (E), Fitting the pMEM to match the empirical data distribution of binarized network states. (F), Determining the relationships between network states using the Dijkstra algorithm on energy values. (G), Interpretation of local minima of the energy on the anatomical level. (H), Simulation of the occurrence of network states belonging to different basins.

Figure 1. 

Illustration of the energy landscape analysis on a network of four regions of interst (ROIs). (A), Selection of ROIs from the source-space signals. (B), Signal filtering in frequency bands of interest. (C), Envelope extraction using the absolute value of the analytical representation of the signal. (D), Binarization of the data. (E), Fitting the pMEM to match the empirical data distribution of binarized network states. (F), Determining the relationships between network states using the Dijkstra algorithm on energy values. (G), Interpretation of local minima of the energy on the anatomical level. (H), Simulation of the occurrence of network states belonging to different basins.

Our results demonstrated that the pMEM provided a good fit to the statistical properties of functional networks in both JME and control groups. JME patients showed reduced numbers of local energy minima and elevated energy values in the theta-, beta-, and gamma-band FPN activity, but not in the SMN. We further demonstrated that the pMEM could be used as a generative model for simulating dysfunctional network dynamics in JME, and as a predictive model for single-patient classification. These findings suggest anatomically- and frequency-specific network abnormalities in JME.

METHODS

Participants

Fifty-two subjects participated in the experiment. Demographic and clinical features of the participants are summarized in Table 1. Twenty-six patients with JME were recruited from a specialist clinic for epilepsy at University Hospital of Wales in Cardiff. Consensus clinical diagnostic criteria for JME were used by an experienced neurologist (Trenité et al., 2013). Inclusion criteria were (1) seizure onset in late childhood or adolescence with myoclonic jerks, with or without absence seizures; (2) generalized tonic-clonic seizures; (3) normal childhood development as assessed on clinical history; and (4) generalized spike wave on EEG and normal structural MRI. Twenty-six healthy control participants with no history of significant neurological or psychiatric disorders were recruited from the regional volunteer panel. All testing was performed with participants taking their usual medication. The study was approved by the South East Wales NHS ethics committee, Cardiff and Vale Research and Development committees, and Cardiff University School of Psychology Research Ethics Committee. Written informed consent was obtained from all participants.

Table 1. 
Demographics of patients with JME and healthy control participants.
 PatientsControls
Number of participants 26 (8 males) 26 (7 males) 
Age median 27 27 
Age range 19–45 18–48 
Seizure type (number of patients) MJ (26)   
Absences (15) – 
GTCS (26)   
Antiepileptic drugs (Number of patients taking the drug) LEV (13), VPA (12),   
LTG (5), TPM (4), – 
ZNM (4) – 
 PatientsControls
Number of participants 26 (8 males) 26 (7 males) 
Age median 27 27 
Age range 19–45 18–48 
Seizure type (number of patients) MJ (26)   
Absences (15) – 
GTCS (26)   
Antiepileptic drugs (Number of patients taking the drug) LEV (13), VPA (12),   
LTG (5), TPM (4), – 
ZNM (4) – 

Note. MJ myoclonic jerks, GTCS generalised tonic clonic seizures, LEV leveiracetam, VPA sodium valproate, LTG lamotrigine, TPM topiramate, ZNM zonisamide.

MEG and MRI data acquisition

All participants underwent separate MEG and MRI sessions. Whole-head MEG recordings were made using a 275-channel CTF radial gradiometer system (CTF Systems, Canada) at a sampling rate of 600 Hz. An additional 29 reference channels were recorded for noise cancellation purposes, and the primary sensors were analyzed as synthetic third-order gradiometers (Vrba & Robinson, 2001). Up to three sensors were turned off during recording because of excessive sensor noise. Subjects were instructed to sit comfortably in the MEG chair while their head was supported with a chin rest and with eyes open focused on a red dot on a gray background. For MEG/MRI co-registration, fiduciary markers that are identifiable on the subject’s anatomical MRI were placed at fixed distances from three anatomical landmarks (nasion, left and right preauricular) prior to the MEG recording, and their locations were further verified using high-resolution digital photographs. The locations of the fiduciary markers were monitored before and after MEG recording. To ensure that the movement artifacts did not dominate the recording, the average Euclidean distance between fiducials was computed for every participant. There was no significant difference between head movements of the JME and control group (t(25) = −1.27, p = 0.22) with the mean head shift of 0.55 cm. Each recording session lasted approximately 5 minutes.

Whole-brain T1-weighted MRI data were acquired using a General Electric HDx 3T MRI scanner and an 8-channel receiver head coil (GE Healthcare, Waukesha, WI) at the Cardiff University Brain Research Imaging Centre with an axial 3D fast spoiled gradient recalled sequence (echo time 3 ms; repetition time 8 ms; inversion time 450 ms; flip angle 20°; acquisition matrix 256 × 192 × 172; voxel size 1 × 1 × 1 mm).

Data pre-processing

Continuous MEG data was first segmented into 2 s epochs. Before segmentation, MEG data was filtered with a 1 Hz high-pass and a 150 Hz low-pass filter to avoid DC step changes between epochs. Every epoch was visually inspected. Those containing major motion, muscle or eye-blink artifact, or interictal spike wave discharges were excluded from subsequent analysis. The artifact-free epochs were then reconcatenated. This artifact rejection procedure resulted in cleaned MEG data with variable lengths between 204 s and 300 s across participants, and the data lengths were comparable between JME patients and controls (t(50) = 1.38, p = 0.17). The 200 s of cleaned MEG data was used in subsequent analysis. For participants with longer than 200 s of cleaned MEG data, a continuous segment of 200 s during the middle of the recording session was used.

Source localization of oscillatory activity in resting-state networks

We analyzed the MEG oscillatory activity using an established source localization method for resting-state networks (Brookes et al., 2011; Hall, Woolrich, Thomaz, Morris, & Brookes, 2013; Muthukumaraswamy et al., 2013). For each participant, the structural MRI scan was coregistered to MEG sensor space using the locations of the fiducial coils and the CTF software (MRIViewer and MRIConverter). The structural MRI scan was segmented, and a volume conduction model was computed using the semirealistic model (Nolte, 2003). The preprocessed MEG data was band-pass filtered with a fourth-order zero phase lag Butterworth filter into four frequency bands: theta 4–8 Hz, alpha 8–12 Hz, beta 13–30 Hz, and low-gamma 35–60 Hz (Niedermeyer, 2005). For each frequency band, we downsampled the data to 250 Hz and computed the inverse source reconstruction using an LCMV beamformer on a 6 mm template with a local spheres forward model in Fieldtrip (version 20161101, http://www.fieldtriptoolbox.org). The atlas-based source reconstruction was used to derive virtual sensors for every voxel in each of the 90 regions of the Automated Anatomical Label (AAL) atlas (Hipp et al., 2012). Each virtual sensor’s time course was then reconstructed.

We focused our analysis on three resting-state networks (Figure 2): the FPN, the DMN, and the SMN in which electrophysiological changes had been reported in patients with epilepsy (Clemens et al., 2013; McGill et al., 2012; Wolf et al., 2015a). Each resting-state network comprised bilateral regions of interest (ROIs) from the AAL atlas identified in previous studies (Rosazza & Minati, 2011; Tewarie et al., 2013a). The FPN included 10 ROIs: middle frontal gyrus (MFG), pars triangularis (PTr), inferior parietal gyrus (IPG), superior parietal gyrus (SPG), and angular gyrus (AG). The DMN included 10 ROIs: orbitofrontal cortex (OFC), anterior cingulate cortex (ACC), posterior cingulate cortex (PCC), precuneus (pCUN), and AG. The SMN included 6 ROIs: precentral gyrus (preCG), postcentral gyrus (postCG), and supplementary motor area (SMA). For each ROI, its representative time course was obtained from the voxel in that ROI with the highest temporal standard deviation. The mean MEG activities of the ROIs of each network were not significantly different between JME patients and controls (FPN: F(1, 50) = 0.75, p = 0.39; DMN: F(1, 50) = 0.21, p = 0.65; SMN: F(1, 50) = 0.15, p = 0.70).

Figure 2. 

The regions of interest (ROIs) of three resting state networks: the fronto-parietal network (FPN), the default mode network (DMN), and the sensori-motor network (SMN). The ROIs were obtained from the 90 AAL atlas (Hipp et al., 2012).

Figure 2. 

The regions of interest (ROIs) of three resting state networks: the fronto-parietal network (FPN), the default mode network (DMN), and the sensori-motor network (SMN). The ROIs were obtained from the 90 AAL atlas (Hipp et al., 2012).

To calculate the oscillatory activity, we applied Hilbert transformation to each ROI’s time course, and computed the absolute value of the analytical representation of the signal to generate an amplitude envelope of the oscillatory signals in each frequency band.

Pairwise maximum entropy model of MEG oscillatory activity

During rest, different brain regions exhibit pairwise co-occurrence of oscillatory activity (Horwitz, 2003) and rapid changes of brain network states (C. J. Stam & Straaten, 2012). To obtain an estimate of network-state transitions and their probabilities, we fitted a pMEM to individual participant’s MEG data, separately for each resting-state network and each frequency band.

According to the principle of maximum entropy, among all probabilistic models describing empirical data, one should choose the one with the largest uncertainty (i.e., entropy), because it makes the minimum assumptions of additional information that would otherwise lower the uncertainty (Yeh et al., 2010). The pMEM estimates the state probability of a network, with its regional activity and regional co-occurrence to be constrained by empirical data. It is equivalent to the Ising model in statistical mechanics (Bialek, 2017). In neuroscience, the pMEM was first used in a seminal study for fitting the distribution of neuronal spiking activity across cells (Schneidman et al., 2006; Tkacik et al., 2006). More recently, the same method was used in fMRI study in which it was shown that the model is capable of estimating the underlying structural connectivity with a higher accuracy than other functional connectivity methods (Watanabe et al., 2013). Later studies using the pMEM for fMRI have identified key characteristics of brain state transition (Kang, Pae, & Park, 2019), perceptual metastability (Watanabe Masuda, Megumi, Kanai, & Rees, 2014b), and the effects of aging (Ezaki et al., 2018). A further advantage of using the pMEM is that various statistical physics theory of the model is available, potentially contributing to the understanding of multivariate data when they are fitted with the pMEM (Bialek, 2017). The current study used this approach to unveil differences between the JME patients and controls in large-scale brain networks (Figure 1). Below we outlined the theoretical background and the fitting procedure. A more detailed description of the pMEM modeling and subsequent energy landscape analysis is available elsewhere (Ezaki et al., 2017).

Consider a resting-state network consisting of N ROIs. For each ROI’s real-valued signal, we thresholded the ROI’s Hilbert envelope according to the median of the amplitude. Data points above the threshold were denoted as high oscillatory power (+1), and data points below the threshold were denoted as low oscillatory power (−1). The oscillatory activity in ROI i (i = 1, …, N) at time t was transformed to a binary time series ri(t), with ri(t) = +1 for high oscillatory activity and ri(t) = −1 for low oscillatory activity. The activity pattern of a N-dimensional binary vector s(t) = [r1(t), r2(t), …, rN(t)] represents the state of the network at time t.

The N-ROI network has a total of 2N possible states sk(k = 1, …, 2N). From the binarized oscillatory activity, we calculated the probability of occurrence of each network state, denoted by Pemp(sk). We further calculated the empirical average activation rate for each ROI 〈riemp and the pairwise co-occurrence between any two ROIs 〈rirjemp:
riemp=1Tt=1Tri(t),
(1)
rirjemp=1Tt=1Tri(t)rj(t),
(2)
where T denotes the number of time points in the data. The fitting procedure aimed to identify a pMEM model that preserves the constraints in Equations 1 and 2 and reproduces the empirical state probability Pemp(sk) with the maximum entropy. It is known that the pMEM follows the Boltzman distribution (Yeh et al., 2010), given by
PpMEM(sk|h,J)=exp(E(sk))k=12Nexp(E(sk)),
(3)
where E(sk) represents the energy of the network state sk, defined by
E(sk)=i=1Nhiri(sk)12i=1Nj=1NjiJijri(sk)rj(sk),
(4)
and ri(sk) refers to the ith element of the network state sk. h and J are the model parameters to be estimated from the data: h = [h1, h2, …, hN] represents the bias in the intensity of the oscillatory activity in each ROI; J = [J11, J12, …, JNN] represents the coupling strength between two ROIs. The average of the activation rate 〈rimod and pairwise co-occurrence 〈rirjmod expected by the pMEM are given by:
rimod=k=12Nri(sk)PpMEM(sk|h,J),
(5)
rirjmod=k=12Nri(sk)rj(sk)PpMEM(sk|h,J).
(6)

We used a gradient ascent algorithm to iteratively update h and J, until 〈rimod and 〈rirjmod match 〈riemp and 〈rirjemp from the observed data, with a stop criterion of 5 × 106 steps. In each iteration, the updates of the parameters were given by hinew = hiold + ϵ(〈riemp − 〈rimod) and Jijnew = Jijold + ϵ(〈rirjemp − 〈rirjmod). The learning rate ϵ was set to 10−8.

As in previous studies (Ezaki et al., 2017, 2018; Watanabe et al., 2014a), we used an accuracy index:
d=(D1D2)/D1
(7)
to quantify the goodness of fit of the pMEM (Figure 3), where
D2=k=12NPemp(sk)log2(Pemp(sk)/PpMEM(sk))
(8)
is the Kullback-Leibler divergence between the probability distribution of the pMEM and the empirical distribution of the network state. D1 represents the Kullback-Leibler divergence between the independent maximum entropy model (MEM) and data. By definition, the independent MEM is restricted to have no pairwise interaction (i.e., J = 0). Therefore, d represents the surplus of the fit of the pMEM over the fit of the independent model. The index d = 1 when the pMEM reproduces the empirical distribution of activity patterns and interactions without errors, and d = 0 when the pairwise interactions do not contribute to the description of the empirical distribution.
Figure 3. 

The pMEM fitting. (A), The occurrence probability of each network state of the FPN from the fitted pMEM (Pmod) was plotted against that from the empirical data (Pemp). Each data point was averaged across juvenile myoclounic epilepsy (JME) patients (red) and controls (blue). (B), The averaged accuracy index d in the JME and control groups for each network and frequency band. Error bars denote the standard errors across participants.

Figure 3. 

The pMEM fitting. (A), The occurrence probability of each network state of the FPN from the fitted pMEM (Pmod) was plotted against that from the empirical data (Pemp). Each data point was averaged across juvenile myoclounic epilepsy (JME) patients (red) and controls (blue). (B), The averaged accuracy index d in the JME and control groups for each network and frequency band. Error bars denote the standard errors across participants.

Energy landscape of resting-state network dynamics

The pMEM parameters h and J determine the energy E(sk) of each network state sk(k = 1, …, 2N), given by Equation 4. It is worth noting that the current study used pMEM as a statistical model to be constructed from the MEG data, not as its literal notion from statistical physics. We did not claim that E(sk) represents th metabolic or physical energy of a biological system. Instead, the concept of the energy of a resting-state network stems from the information theory. Here, E(sk) indicates the model prediction of the inverse appearance probability of the state sk under the empirical constraints of regional activity (parameter h) and regional interactions (parameter J). For instance, if E(si) < E(sj), the pMEM predicts that the network activity pattern is more likely to be at the state si than sj.

For each resting-state network and each frequency band, we depicted an energy landscape as a graph of the energy function across the 2N possible network states sk, characterizing state probabilities and state transitions from the perspective of attractor dynamics (Watanabe et al., 2014a). Because the computational cost increases dramatically with the size of a network, we estimated an energy landscape separately for each resting-state network.

The energy landscape of a network was defined by two factors: the energy E(sk) of each network state, and an adjacency matrix defining the connectivity between network states. Two states were defined to be adjacent, or directly connected, if and only if just one ROI of the network had different binarized oscillatory activity (high vs. low). In other words, two states are adjacent when they have a Hamming distance of 1 between their binary activity vectors. For example, for a network with four ROIs, states [−1, −1, −1, +1] and [−1, −1, +1, +1] are adjacent, and states [−1, −1, −1, +1] and [−1, −1, +1, −1] are not.

Quantitative measures of energy landscape

We used three measures to understand the differences in the energy landscape between JME patients and healthy controls: (1) the number of local energy minima of within-network dynamics, (2) the relative energy of the local minima, and (3) the generative basin duration at significant minima.

Number of local energy minima

A local energy minimum was defined as the network state with a lower energy value than all its adjacent states. Because lower energy corresponds to a higher probability of occurrence, network states of local minima can be likened as attractors in attractor dynamics. For each participant, we exhaustively searched through the 2N network states to identify all the local minima of the participant’s energy landscape. We then compared the number of local energy minima between JME and control groups (Figure 4).

Figure 4. 

The average number of local minima in the JME and control groups. Error bars denote the standard errors across participants.

Figure 4. 

The average number of local minima in the JME and control groups. Error bars denote the standard errors across participants.

Relative energy of the local minima

The number of local minima is determined by the energy difference between network states and their neighbors (i.e., a minimum has a lower energy level than all its neighbors). On the other hand, the energy value of a specific state is determined by its occurrence probability (Equation 3). Therefore, theoretically, the two measures had no direct dependency. The energy values of local minima on aggregated energy landscapes indicate the ease of transition from one stable state to another (Ezaki et al., 2017; Kang et al., 2019).

We calculated the mean energy Ē(sk) of each network state sk averaged across all participants. Then, we used the mean energy to depict an aggregated landscape, which allowed us to identify common energy minima shared between JME patient and control groups. To test whether each local minimum in the aggregated energy landscape is a characteristic feature of the observed data, we conducted non-parametric permutation tests on the mean energy values. For each resting-state network and each frequency band, we conducted 1,000 permutations. In each permutation, we randomly shuffled the pMEM parameters h and J (between ROIs and ROI pairs, respectively) that were fitted to individual participants. We then calculated an averaged energy landscape across all participants based on the shuffled parameters. This gave us a sampling distribution of the energy of each network state, under the null hypothesis that the energy values are not related to the observed oscillatory activities or observed pairwise regional co-occurrence. For each local minimum of the aggregated landscape, the level of significance (p value) of that local minimum’s energy was estimated by the fraction of the permutation samples that were higher than the mean energy Ē(sk) of that network state in the empirical data without shuffling. To account for the multiple statistical tests that were performed for all the local minima of each network, we evaluated the results using a Bonferroni-corrected threshold (p < 0.05) for significance.

Because the shape of an energy landscape was partly determined by the global minimum (Ezaki et al., 2018; Watanabe et al., 2014a), for each participant, we calculated the energy difference between a significant local minimum and the global energy minimum (i.e., the state with the lowest energy value on the landscape). We then compared this relative energy of the within-network local minima between JME patients and healthy controls. From the networks with significant alternations of relative energy values in JME patients, we constructed a disconnectivity graph to describe clusters of local minima and the relationships between them (Becker & Karplus, 1997), where a cluster represents a group of local minima with high probabilities of subsequent occurrences (Ezaki et al., 2017).

Basin duration at significant minima

On the aggregated energy landscape, the energy basin for each significant local minimum was identified using an existing method (Watanabe et al., 2014a). We started at an arbitrary network state and moved downhill on the energy landscape to one of its neighboring states with the lowest energy, until a local minimum was reached. The starting state is then assigned to the basin of the resulting local minimum. We repeated this procedure for all network states as the starting state.

We used the fitted pMEM as a generative model to simulate the dynamical changes in each resting-state network, and estimated the duration of the basin of each local minimum in the simulated dynamics. Similar to previous studies, we employed the Metropolis-Hastings algorithm to simulate time courses of network activity (Hastings, 1970). Each simulation started with a random network state sk. On each time step, one of the current state’s N neighboring state sk was selected with a probability of 1/N as the potential target of state transition, and the state transition occurred with a probability of exp[E(sk) − E(sk)] when E(sk) > E(sk) or 1 otherwise. For each participant, each network, and each frequency band, we simulated 20,000 time steps, and discarded the first 1,000 time steps to minimize the effect of initial condition. From the remaining 19,000 time steps, we calculated the proportion of duration of the network states that belongs to each energy basin.

Classification of individual patients based on energy values

To investigate the predictive power of pMEM energy measures, we used a support vector machine (SVM) classifier with a radial basis function kernel and a leave-one-out cross-validation procedure to classify individual JME patients and controls. The trade-off between errors of the SVM on training data and margin maximization was set to 1. For each resting-state network and each frequency band, the feature space for classification included the energy values of all the significant local minima. In each cross-validation fold, one participant was first removed and the remaining participants’ data were used as a training set to train the classifier. To avoid overfitting, the feature space (i.e., the local minima) was identified from the aggregated energy landscape constructed from the participants in the training set. The participant left out was then classified into one of the two groups (patients or controls). Classification performance was evaluated by the proportion of correctly classified participants over all cross-validations.

We used permutation tests to evaluate the classification results. The significance of each classification was determined by comparing the observed classification accuracy with its null distribution under the assumption of no difference between patients and controls. The null distribution was generated by 1,000 random permutations of leave-one-out classification results, with group labels shuffled in each permutation. We obtained a permutation p value by calculating the fraction of the permuted samples exceeding the observed classification accuracy.

Software and Data Accessibility

The scripts for analysis used in the current study are open-source and freely available online (https://github.com/dokato/energy_landscape). Raw MRI and MEG data are not publicly available due to informed consent restrictions concerning confidential patient information.

RESULTS

A summary of participant demographics and clinical characteristics is given in Table 1. The JME and control groups were well matched for age (F(1, 51) = 0.13, p = 0.72) and gender (p = 0.31, χ2 test). For each participant, we performed source localization of preprocessed MEG resting-state data and estimated oscillatory activity (i.e., Hilbert envelope) in each of the 90 ROIs from the AAL atlas, separately in the theta (4–8 Hz), alpha (8–12 Hz), beta (13–30 Hz), and low-gamma (35–60 Hz) bands. We focused our analysis on the differences between JME patients and controls in three resting-state networks (Figure 2): the FPN, the DMN, and the sensorimotor network (SMN).

Fitting of pairwise maximum entropy models (pMEM) to MEG oscillatory activity

We thresholded an ROI’s oscillatory amplitude at each time point t to assign the binary states of “high” (+1) or “low” (−1) activity. The state of a network at time t was then represented by a binary vector, consisting of the binarized activity of all the ROIs in the network. We fitted a pMEM to the series of binarized network oscillatory activities, separately for each participant, each resting-state network, and each frequency band (Equation 3, and see Methods for details). For a network of N ROIs, there are a total of 2N possible states. The pMEM provides a statistical model of the occurrence probabilities of the 2N network states, while it satisfies the empirical constraints of mean regional activities at each ROI and pairwise co-occurrence between each pair of ROIs within the network.

To evaluate the model fit, we compared the predicted and observed occurrence probabilities of the 2N possible network states, averaged across the participants in each group. There was a good agreement between the model predictions and observed data across networks in the JME (R2 > 0.90 in all networks and frequency bands, based on a log-log regression, Figure 3) and control groups (R2 > 0.89). We further used an accuracy index to quantify the goodness of fit of the pMEM (Equation 7). The accuracy index was calculated as the percentage of improvement of the pMEM fit to the empirical data compared with a null model, which assumed no pairwise co-occurrence between ROIs (i.e., an independent maximum entropy model). The pMEM achieved high accuracy indeces in both JME patients and controls (Figure 3). A Mann-Whitney U-test on accuracy indeces showed no significant main effect of group (JME vs. controls: U = 266.0, p = 0.19), suggesting the robustness of the pMEM on MEG oscillatory activities in both patients and controls. As determined by nonparametric repeated-measures Friedman test, there were main effects of the networks (χ2 = 87.5, p < 0.00001) and the frequency bands (χ2 = 46.27, p < 0.00001), suggesting that the distinct properties of the networks and information carried by the frequency bands affected the goodness of fit.

Inferences from pMEM energy landscape

The fitted pMEM yielded an energy value for each network state (Equation 4). We used energy values from the pMEM to depict an energy landscape of the network. The energy landscape is a graph representation of energy values from all possible network states (Figure 1F). We defined two network states being adjacent if there is one and only one ROI whose binarized activity (i.e. +1 or −1) is the opposite between the two states. According to the pMEM (Equations 3 and 4), network states with a higher energy would occur less frequently than those with a lower energy. As a result, transitions from high to low energy states would be more likely to occur than that from low to high energy states. Here, we examined the differences in three quantitative measures of energy landscape between patients with JME and controls: (1) the number of energy minima, (2) the relative energy values at the local minima, and (3) the generative basin duration at significant minima.

Number of energy minima

We located local minima on the energy landscape, defined as the network states with lower energy than all their adjacent states. Because a local minimum state would have a higher occurrence probability than all of its neighboring states, transitions of network states near an energy minimum is akin to a fixed point attractor in a deterministic dynamical system, and the number of energy minima quantifies the degree of multistability of a network.

We calculated the number of local minima for each participant (Figure 4) and compared it between groups, resting-state networks, and frequency bands with a repeated-measures ANOVA. Compared with controls, JME patients had significantly less local energy minima (F(1, 50) = 7.602, p = 0.008). Across all participants, there were significant main effects of the resting-state network (F(1.52, 76.25) = 99.89, p < 0.00001, Greenhouse corrected) and frequency band (F(2.83, 141.57) = 21.08, p < 0.00001). No significant network by group (F(1.52, 76.25) = 3.15, p = 0.07) or frequency band by group (F(2.83, 141.57) = 2.12, p = 0.11) interaction was observed. These results suggested that MEG oscillatory activities in JME patients had altered multistability in some networks and frequency bands.

Relative energy values of the local minima

To identify common energy minima at the group level, we averaged across all participants the energy value of each network state and identified the energy minima on the aggregated energy landscape. In all three resting-state networks (FPN, DMN, and SMN) and all frequency bands, permutation tests showed that the energy values of two network states, “all off” (i.e., all ROIs had low oscillatory activities [−1, −1, …, −1]) and “all on” (i.e., all ROIs had high oscillatory activities [+1, +1, …, +1]), did not differ significantly from those from a randomly shuffled energy landscape (p > 0.88, Bonferroni corrected). That is, the observed energy values at these two minima were not significantly sensitive to regional activation and pairwise coactivation in empirical data (see Methods for details). In addition, the “all off” state was also the global minimum of the energy landscape at both group and individual levels, which had the lowest energy value in all network states.

For each significant local minimum state that survived the permutation test, we calculated the relative energy difference between the local minimum and the “all off” state (i.e., the global minimum) for the individual participants. Then, we compared the obtained relative energy values between the JME and control groups. This subtraction step controlled for the individual variability in the occurrence probability of the global minimum state (Watanabe et al., 2014a).

In the FPN, the relative energy values at the local minima were significantly higher in JME patients than in controls in the theta band (Figure 5A, F(1, 50) = 18.90, p < 0.0001), beta band (Figure 5B, F(1, 50) = 15.43, p = 0.0002), and gamma band (Figure 5C, F(1, 50) = 7.2558, p = 0.009), but not in the alpha band (F(1, 50) = 0.80, p = 0.37). The aggregated energy landscapes in the beta and gamma bands contained the same set of 14 local minima. Post hoc tests showed that all the 14 local minima states had higher relative energy values in JME patients than controls in the beta band, and 5 of the 14 local minima states showed a significant group differences in the gamma band (p < 0.05, Šidák correction). The theta-band energy landscape contained six local minima states, which were a subset of the 14 local minima in the higher frequency bands, and all had higher relative energy values in JME patients than controls.

Figure 5. 

Relative energy values of the local minima in (A), theta FPN. (B), beta FPN. (C), gamma FPN, and (D), beta DMN. At the top of each panel, the disconnectivity graph showed the relative energy values of local minima from the aggregated energy landscape across all participants. The end of each branch on the disconnectivity graph represent a local minimum. The middle of each panel showed the network states of the corresponding local minima. White boxes denote high oscillatory activity (i.e., a binary value of +1) and gray box denote low oscillatory activity (i.e., a binary value of −1). The bottom of each panel showed the t-values from two sample t-tests (JME patients vs. controls) on the relative energy values of each local minimum. Asterisks indicate significant difference between JME patient and control groups (p < 0.05, FDR corrected).

Figure 5. 

Relative energy values of the local minima in (A), theta FPN. (B), beta FPN. (C), gamma FPN, and (D), beta DMN. At the top of each panel, the disconnectivity graph showed the relative energy values of local minima from the aggregated energy landscape across all participants. The end of each branch on the disconnectivity graph represent a local minimum. The middle of each panel showed the network states of the corresponding local minima. White boxes denote high oscillatory activity (i.e., a binary value of +1) and gray box denote low oscillatory activity (i.e., a binary value of −1). The bottom of each panel showed the t-values from two sample t-tests (JME patients vs. controls) on the relative energy values of each local minimum. Asterisks indicate significant difference between JME patient and control groups (p < 0.05, FDR corrected).

In the DMN, there were trends of higher relative energy values in the JME patients than controls in the beta band (F(1, 50) = 3.68, p = 0.06) and gamma band (F(1, 50) = 3.81, p = 0.06), and no significant difference in the theta band (F(1, 50) = 0.01, p = 0.92) or alpha band (F(1, 50) = 0.82, p = 0.37). One local minimum in the beta band, comprising co-activation in bilateral mPFC and ACC (Figure 5D), showed a group difference in post hoc tests at an uncorrected threshold (t(50) = 2.34, p = 0.03). In the SMN, there was no significant group difference in the relative energy values (theta: F(1, 50) = 1.26, p = 0.27; alpha: F(1, 50) = 0.06, p = 0.81; beta: F(1, 50) = 0.002, p = 0.97; gamma: F(1, 50) = 0.12, p = 0.73).

Overall, JME patients had higher relative energy values than controls in selective resting-state networks and frequency bands. This result indicates that some local minima on the aggregated energy landscape were less stable (i.e., having a higher relative energy level) in JME patients than controls.

Basin duration at significant minima

Each local minimum of an energy landscape is accompanied by a basin, which includes the local minimum itself and its neighboring states from which the local minimum is relatively easily reached (Ezaki et al., 2017). Therefore, the proportion of time for which each basin is visited gives a granular description of network dynamics. For each of the group-level significant minima on the aggregated energy landscape, we identified all the network states belonging to the same basin. For each participant, we then used the fitted pMEM to numerically simulate network dynamics, and calculated the proportion of time for which the simulated network activities visit each basin. The rationale to simulate basin durations is twofold. First, our simulation demonstrated the feasibility of the derived energy landscape to be used as a generative model of network dynamics. Second, because we removed MEG epochs strongly affected by artífacts, the source reconstructed data was not fully continuous in time, and hence basin duration estimated directly from the empirical data would be less accurate.

In the FPN, simulations showed that network dynamics in JME patients contained shorter basin duration at those significant local minima than controls in the theta (F(1, 50) = 42.72, p < 0.000001), beta (F(1, 50) = 10.49, p = 0.002), and gamma (F(1, 50) = 6.18, p = 0.016) bands, but not in the alpha band (F(1, 50) = 3.92, p = 0.053). There was no significant group difference in the basin duration in the DMN (theta: F(1, 50) = 0.015, p = 0.90; alpha: F(1, 50) = 2.67, p = 0.11; beta: F(1, 50) = 2.76, p = 0.10; gamma: F(1, 50) = 3.12, p = 0.08) or SMN (theta: F(1, 50) = 0.09, p = 0.76; alpha: F(1, 50) = 0.31, p = 0.58; beta: F(1, 50) = 1.59, p = 0.21; gamma: F(1, 50) = 1.25, p = 0.27).

Classification of individual patients

We used a leave-one-out cross-validation procedure for a binary classification of participant groups (JME patients and healthy controls), using the relative energy values of local minima as features. Consistent with the group comparisons (Figure 5), the relative energy values obtained from the fitted pMEM showed significant predictive power, with high classification accuracies from theta-band FPN (92.3%, p < 0.001, permutation test) and gamma-band FPN (67.3%, p = 0.012) (Figure 6). The classification based on the energy values from theta-band FPN achieved high specificity (89.3%) and sensitivity (94.8%). For the classification based on gamma-band FPN features, the specificity and sensitivity were 71.4% and 64.5%, respectively. The classification accuracy in the SMN, DMN, and other frequency bands of the FPN was not significant (p > 0.26, permutation test).

Figure 6. 

Support vector machine leave-one-out binary classification accuracy of JME patients versus controls. The energy values of the local minima were used as features for classifiers. Blue data points denote the mean classification accuracy. Black lines denote the 95% confidence level under the null hypothesis of no difference between the groups, based on 1,000 permutations of randomly shuffled labels of the data.

Figure 6. 

Support vector machine leave-one-out binary classification accuracy of JME patients versus controls. The energy values of the local minima were used as features for classifiers. Blue data points denote the mean classification accuracy. Black lines denote the 95% confidence level under the null hypothesis of no difference between the groups, based on 1,000 permutations of randomly shuffled labels of the data.

DISCUSSION

We proposed a pMEM approach to quantify the dynamics of MEG oscillatory activity and applied this method to derive energy landscape measures, quantifying the abnormal statistical characteristics of resting-state networks in JME patients. The number of within-network local minima from individual participant’s energy landscape indicates the degree of multistability from an attractor network perspective (Kelso, 2012) on MEG oscillatory power. The local minima are defined here, and should always be interpreted, in the context of a specific resting-state brain network. The energy values of minima on aggregated energy landscapes indicate the ease of transition from one stable state to another (Ezaki et al., 2017; Kang et al., 2019), and their effects on network dynamics was demonstrated in the simulation of basin duration. Furthermore, the activation profiles of local minima provided key anatomical insights into functional configurations of cortical networks that differ between JME patients and controls. Our approach described network abnormalities in multivariate data from a statistical account. This extended previous research on the temporal evolution of system dynamics leading to seizures, which measures chaoticity (Iasemidis et al., 1990) or entropy (Song et al., 2012) in single or combined channels.

In this study, we found that patients with JME showed altered pMEM-derived energy landscapes in selective resting-state networks and frequency bands (Figure 7). For the energy landscapes estimated at the individual level, JME patients exhibited fewer local minima than controls (Figure 4). For the aggregated energy landscapes estimated across participants, JME elevated relative energy values at the local minima of the FPN (theta, beta, and gamma bands) oscillatory activities (Figure 5). Our results confirmed the abnormalities of electrophysiological signals in JME (Aliberti, Grünewald, Panayiotopoulos, & Chroni, 1994), and provided new insights into JME pathophysiology affecting selective functional network configurations.

Figure 7. 

A schematic diagram of altered energy landscape of MEG oscillatory power in JME patients (left) compared with controls (right). In selective functional networks and frequency bands, JME patients exhibited less local energy minima and elevated energy values (e.g., in theta-band FPN), suggesting that resting-state networks exhibit changes in the degree of multistability and in the ease of state transitions.

Figure 7. 

A schematic diagram of altered energy landscape of MEG oscillatory power in JME patients (left) compared with controls (right). In selective functional networks and frequency bands, JME patients exhibited less local energy minima and elevated energy values (e.g., in theta-band FPN), suggesting that resting-state networks exhibit changes in the degree of multistability and in the ease of state transitions.

The fitted pMEM defined the energy values of all activity states of a network, from which an energy landscape of the network was depicted (Ezaki et al., 2017). Because a local minimum of the energy landscape refers to a network state with higher occurrence probability than its neighboring states, the fewer number of local minima and elevated energy values in JME suggested alterations in the multi-stable dynamics of the brain networks that may be prone to perturbation and ictogenesis, in line with the dynamical disease account for epilepsy (da Silva et al., 2003; Elger et al., 2000; C. Stam, 2005). The energy landscape further allowed to characterize clusters of energy minima and their hierarchies in terms of the disconnectivity graphs (Figure 5). In the FPN, the energy minima with bilateral high activation in the frontal or parietal regions were clustered separately and interleaved with lateralized energy minima (i.e., high activation in unilateral ROIs). This may indicate that network states with lateralized high activation represent transition statuses between frontal and parietal dominant states. In contrast, the DMN energy minima contained coactivation in bilateral ROIs, consistent with the evidence of strong interhemispheric and long-range connectivity in the DMN during awake states (Baker et al., 2014; Salvador et al., 2005).

Our results highlighted JME as a distributed network disorder involving frontal and parietal lobes (Fernandez et al., 2011; Niso et al., 2015; Wolf et al., 2015b). JME patients commonly exhibit impaired frontal cognitive functions (Piazzini, Turner, Vignoli, Canger, & Canevini, 2008), including working memory (Swartz, Halgren, Simpkins, & Syndulko, 1994), decision-making (Zamarian et al., 2013), response inhibition (S.-Y. Kim et al., 2007), and verbal fluency (O’Muircheartaigh et al., 2011). Demanding cognitive efforts during visuomotor coordination and decision-making can provoke myoclonic seizures in JME patients (Yacubian & Wolf, 2014), and the degree of cognitive dysfunctions were associated with fronto-parietal BOLD fMRI activity and connectivity (Vollmar et al., 2011). Cortical and subcortical pathology may underlie the cognitive phenotype in JME. Activities in the lateral parietal cortex and precuneus have a dominant role in initiating and sustaining characteristic spike-and-wave discharges in JME (Lee et al., 2014). MR spectroscopy imaging of JME patients has identified reduced N-Acetyl aspartate concentrations in the frontal lobe and the thalamus (Savic, Lekvall, Greitz, & Helms, 2000; Zhang, Li, Hong, & Zou, 2016), which, together with widespread cortical morphological abnormalities (Ronan et al., 2012), indicates dysfunctions in the corticothalamic loops in JME (Hattingen et al., 2014). Further research should extend our results to associate specific abnormal energy minima to JME patients’ cognitive and behavioral phenotypes.

We further demonstrated that the pMEM and energy landscapes can be used as a generative model to simulate the duration of the network activity in each energy basin (Figure 1) and as a predictive model for single-patient classification (Figure 6) beyond simple descriptive modeling (Shmueli, 2010): it allowed us to combine measures from multiple energy minima to make inferences at an individual level. Such analysis, as demonstrated in the current study, would be useful in clinical applications for identifying patients from controls, or for detecting changes in electrophysiological data prior to seizure onset in future studies (Song & Zhang, 2013). In addition, because classification-based analysis makes no assumption about data variances or distributions, it is a more stringent test than conventional statistical methods and provides accurate estimates of between-group differences (Kim & Oertzen, 2018). The normalized energies of the theta-band FPN minima achieved the best classification results (>90%), comparable with other studies (Goker et al., 2012) and consistent with our hypothesis of selective abnormalities of oscillatory activity in JME. Indeed, pathological theta oscillations were reported as a hallmark of idiopathic generalized epilepsy (Clemens, 2004), possibly owing to the involvement of the thalamus in initiating or facilitating theta oscillations through thalamocortical coherence (Sarnthein, Morel, Von Stein, & Jeanmonod, 2003).

The energy landscape measures for the SMN did not significantly differ between JME patients and controls. This result might seem counterintuitive, given that motor cortex hyperexcitability has been reported in JME (Badawy, Curatolo, Newton, Berkovic, & Macdonell, 2006). Nevertheless, previous research on resting-state functional connectivity also showed the lack of altered connectivity in the motor cortex in JME (Elshahabi et al., 2015; Li et al., 2015; Liao et al., 2011). Our results suggested that the network states (i.e., patterns of coactivation) in the SMN, comprising pre- and postcentral gyri as well as SMA, were not affected by JME during rest. However, this result does not rule out the possibility of network dysfunction in the motor circuit under stimulation or perturbation (Vollmar et al., 2011).

Our study provides new methods for studying the dynamics of MEG oscillatory activity. We showed that MEG oscillatory activity in resting-state networks was accurately described by the pMEM (Figure 3) and that the model fits were comparable between JME patients and controls. The pMEM was originally developed in the field of statistical mechanics and has been applied to population of spiking neurons (Schneidman et al., 2006; Yeh et al., 2010). More recently, it has been applied to quantify the dynamics of BOLD fMRI data (Ashourvan et al., 2017; Ezaki et al., 2018; Gu et al., 2018; Watanabe et al., 2013, 2014a). However, achieving satisfactory pMEM fitting requires a large number of data samples (Ezaki et al., 2017; Macke, Murray, & Latham, 2012). Because of the low temporal resolution of the BOLD signal, the applications of the pMEM to fMRI signals often need long scanning time that may be unrealistic for clinical populations, or to concatenate data across participants that limits the possibility of individual-level inferences (Ashourvan et al., 2017; Watanabe et al., 2014a). Here, we highlighted the feasibility and benefits of fitting the pMEM to MEG oscillatory power, which provided anatomically specific and frequency-dependent results. Capitalizing on the high sampling rate of MEG, we showed that one can make inferences on energy landscapes at the individual level from a short recording session that was well tolerated by patients. Future studies could use longer recording sessions to systematically examine the effect of data length on pMEM fitting to MEG data.

Other methods are available to describe transient network dynamics. Microstate analysis from scalp EEG has identified successive short time periods during which the configuration of the scalp potential field remains semistable (Baker et al., 2014), and the spatial patterns of EEG microstates have been mapped onto distinct mental states (Brodbeck et al., 2012; Michel & Koenig, 2018). Recent studies using hidden Markov models (HMM) characterized whole-brain spontaneous activity and identified hidden states with spatiotemporal patterns at durations of 100–200 ms (Quinn et al., 2018). Both microstate and HMM analyses are based on time-windowed approaches and provide abstractions of the interactions within large-scale networks. In the current study, we defined the state of a network as an instantaneous snapshot of regional activities, and the pMEM provided a probabilistic model of the network states with minimum assumptions.

There are several limitations of this study. First, to quantify network dynamics as the occurrence probability of a finite number of network states, the oscillatory power in each ROI needed to be binarized (i.e., high vs. low activity), similar to other functional connectivity studies (Liao et al., 2010). The binarization procedure for applying pMEM in neuroscience differs between data modalities. For single-unit recording and local field potentials (LFPs) (Tang et al., 2008; Yeh et al., 2010; Yu et al., 2011), a threshold based on signal variance was applied to continuous data to identify active states (spikes in single units or negative deflections in LFPs). For resting-state fMRI data, a threshold based on the mean of BOLD responses was used (Kang et al., 2019; Watanabe et al., 2014b). Unlike spiking trains or LFPs that have a clear definition of neuronal activity status, MEG oscillatory power reflects the level of synchronized activity in macroscopic neural populations, which, as a continuous measure, does not impose an a priori threshold for active/inactive binarization. The current study used the median of the oscillatory power envelope from each ROI as the threshold to binarize MEG source reconstructed data. The use of a median split is robust to signal outliers. Furthermore, our approach allows a common statistical criterion adaptive across regions and participants, appropriate for a potentially heterogeneous ensemble (Deco et al., 2012). Future research could consider more complex quantification scheme such as ternary quantization that reduces oscillatory power to ternary values (Zhu, Han, Mao, & Dally, 2016).

Second, the model fitting procedure for pMEM is computationally intensive. Currently, it is practically possible to exactly fit the pMEM to a network of approximately 15 ROIs, because the number of network states increases exponentially with more ROIs. As a result, the current study focused on the dynamics within well-established large-scale resting-state networks, rather than a whole-brain network comprising all the regions. Other approximate model fitting procedures may allow us to extend our approach to larger networks with more ROIs (Ezaki et al., 2017), which is beyond the scope of the current study. To facilitate future research, we have made our analysis scripts open source and freely available (https://github.com/dokato/energy_landscape).

Third, the current study chose, a priori, the AAL template for cortical parcellation. The AAL atlas is based on anatomical landmarks (Rolls, Joliot, & Tzourio-Mazoyer, 2015; Tzourio-Mazoyer et al., 2002) and commonly used in MEG resting-state analysis (Hillebrand et al., 2016; Papanicolaou et al., 2006; Routley, Singh, Hamandi, & Muthukumaraswamy, 2017). Previous studies have defined resting-state networks, including the ones used in our study, with the ROIs from the AAL atlas (Rosazza & Minati, 2011; Tewarie et al., 2013a). It is worth noting that there is an abundant group of atlases for cortical parcellation with various levels of granularity (Desikan et al., 2006; Destrieux, Fischl, Dale, & Halgren, 2010; Gordon et al., 2014; Klein & Tourville, 2012), and energy landscape measures from a network may change with different ROI definitions from an alternative atlas. Future research employing the pMEM for MEG needs to make informed decisions on the choice of parcellation scheme based on specific research questions and intended networks.

Fourth, the sample size in the current study is sufficient for comparing and classifying between JME patient and control groups. However, JME patients are often on a phenotypic spectrum, with variations among seizure frequencies, epileptic traits, and treatment response (Baykan & Wolf, 2017). A larger clinical cohort with comprehensive neuropsychological assessments is necessary to investigate whether our energy landscape approach is sensitive to the quantitative spectrum of JME. Moreover, the scanning protocol of this experiment did not enable continuous movement tracking. As a result, we could not directly compare the exact level of movements between groups. Nevertheless, even if the residual movement artifacts in MEG data did differ between patients and controls, they would affect multiple networks and hence could not readily explain the network-specific group differences in energy landscape measures.

In conclusion, by fitting a pMEM to MEG oscillatory activity, we showed that JME patients exhibited atypical energy landscapes in selective brain networks and frequency bands, with a smaller number of local minima of the energy and elevated energy levels leading to altered multistable network dynamics. We further demonstrated that the pMEM and energy landscape offered generative and predictive power for discriminating between JME patients and controls. These results have the potential to be exploited in future diagnostic and pharmacological studies for a mechanistic understanding of ictogenesis in JME.

AUTHOR CONTRIBUTIONS

Dominik Krzemiński: Conceptualization; Formal analysis; Investigation; Methodology; Software; Visualization; Writing - Original Draft; Writing - Review & Editing. Naoki Masuda: Conceptualization; Methodology; Software; Supervision; Writing - Review & Editing. Khalid Hamandi: Data curation; Writing - Review & Editing. Krish D. Singh: Methodology; Supervision; Writing - Review & Editing. Bethany Routley: Data curation; Writing - Review & Editing. Jiaxiang Zhang: Conceptualization; Formal analysis; Methodology; Project administration; Supervision; Validation; Writing - Original Draft; Writing - Review & Editing.

FUNDING INFORMATION

Krish D. Singh, Medical Research Council (http://dx.doi.org/10.13039/501100000265), Award ID: MR/K005464/1. Dominik Krzemiński, Engineering and Physical Sciences Research Council (http://dx.doi.org/10.13039/501100000266), Award ID: EP/N509449/1. Jiaxiang Zhang, European Research Council, Award ID: 716321. Bethany Routley, Medical Research Council (http://dx.doi.org/10.13039/501100000265), Award ID: MR/K501086/1. Khalid Hamandi, Health Care Research Wales.

TECHNICAL TERMS

     
  • Juvenile myoclonic epilepsy (JME):

    The most common form of generalized epilepsy, characterized by myoclonic jerks, generalized tonic-clonic seizures, and absence seizures.

  •  
  • Pairwise maximum entropy model (pMEM):

    A maximum entropy model that takes into account both average activity and pairwise interactions.

  •  
  • Network state:

    A binary vector that denotes high (+1) and low (−1) oscillatory activity in each region of a network.

  •  
  • Energy landscape:

    A graphical representation of the mapping from brain states to their energy values defined by the pairwise maximum entropy model.

  •  
  • Maximum entropy model (MEM):

    A statistical model of data with the highest entropy among all those that satisfy empirical constraints.

  •  
  • Local energy minima:

    The networks state with a lower energy value than all its neighbors.

  •  
  • Energy basin:

    A group of network states that are the most closely associated with a common local minimum.

REFERENCES

Aliberti
,
V.
,
Grünewald
,
R. A.
,
Panayiotopoulos
,
C. P.
, &
Chroni
,
E.
(
1994
).
Focal electroencephalographic abnormalities in juvenile myoclonic epilepsy
.
Epilepsia
,
35
(
2
),
297
301
. https://doi.org/10.1111/j.1528-1157.1994.tb02433.x
Ashourvan
,
A.
,
Gu
,
S.
,
Mattar
,
M. G.
,
Vettel
,
J. M.
, &
Bassett
,
D. S.
(
2017
).
The energy landscape underpinning module dynamics in the human brain connectome
.
NeuroImage
,
157
,
364
380
.
Babloyantz
,
A.
, &
Destexhe
,
A.
(
1986
).
Low-dimensional chaos in an instance of epilepsy
.
Proceedings of the National Academy of Sciences
,
83
(
10
),
3513
3517
. https://doi.org/10.1073/pnas.83.10.3513
Badawy
,
R.
,
Curatolo
,
J. M.
,
Newton
,
M.
,
Berkovic
,
S. F.
, &
Macdonell
,
R. A.
(
2006
).
Sleep deprivation increases cortical excitability in epilepsy
.
Neurology
,
67
(
6
),
1018
1022
. https://doi.org/10.1212/01.wnl.0000237392.64230.f7
Baker
,
A. P.
,
Brookes
,
M. J.
,
Rezek
,
I. A.
,
Smith
,
S. M.
,
Behrens
,
T.
,
Probert Smith
,
P. J.
, &
Woolrich
,
M.
(
2014
).
Fast transient networks in spontaneous human brain activity
.
eLife
,
3
,
e01867
. https://doi.org/10.7554/eLife.01867
Baykan
,
B.
, &
Wolf
,
P.
(
2017
).
Juvenile myoclonic epilepsy as a spectrum disorder: A focused review
.
Seizure
,
49
,
36
41
. https://doi.org/10.1016/j.seizure.2017.05.011
Becker
,
O. M.
, &
Karplus
,
M.
(
1997
).
The topology of multidimensional potential energy surfaces: Theory and application to peptide structure and kinetics
.
The Journal of Chemical Physics
,
106
(
4
),
1495
1517
. https://doi.org/10.1063/1.473299
Berkovic
,
S. F.
,
Howell
,
R. A.
,
Hay
,
D. A.
, &
Hopper
,
J. L.
(
1998
).
Epilepsies in twins: Genetics of the major epilepsy syndromes
.
Annals of Neurology
,
43
(
4
),
435
445
.
Betting
,
L. E.
,
Mory
,
S. B.
,
Li
,
L. M.
,
Lopes-Cendes
,
I.
,
Guerreiro
,
M. M.
,
Guerreiro
,
C. A.
, &
Cendes
,
F.
(
2006
).
Voxel-based morphometry in patients with idiopathic generalized epilepsies
.
NeuroImage
,
32
(
2
),
498
502
. https://doi.org/10.1016/j.neuroimage.2006.04.174
Bialek
,
W.
(
2017
).
Perspectives on theory at the interface of physics and biology
.
Reports on Progress in Physics
,
81
(
1
),
012601
.
Brodbeck
,
V.
,
Kuhn
,
A.
,
von Wegner
,
F.
,
Morzelewski
,
A.
,
Tagliazucchi
,
E.
,
Borisov
,
S.
, …
Laufs
,
H.
(
2012
).
EEG microstates of wakefulness and NREM sleep
.
NeuroImage
,
62
(
3
),
2129
2139
. https://doi.org/10.1016/j.neuroimage.2012.05.060
Brookes
,
M. J.
,
Woolrich
,
M.
,
Luckhoo
,
H.
,
Price
,
D.
,
Hale
,
J. R.
,
Stephenson
,
M. C.
, …
Morris
,
P. G
.
(
2011
).
Investigating the electrophysiological basis of resting state networks using magnetoencephalography
.
Proceedings of the National Academy of Sciences
,
108
(
40
),
16783
16788
. https://doi.org/10.1073/pnas.1112685108
Caeyenberghs
,
K.
,
Powell
,
H.
,
Thomas
,
R.
,
Brindley
,
L.
,
Church
,
C.
,
Evans
,
J.
, …
Hamandi
,
K.
(
2015
).
Hyperconnectivity in juvenile myoclonic epilepsy: A network analysis
.
NeuroImage: Clinical
,
7
,
98
104
. https://doi.org/10.1016/j.nicl.2014.11.018
Camfield
,
C. S.
,
Striano
,
P.
, &
Camfield
,
P. R.
(
2013a
).
Epidemiology of juvenile myoclonic epilepsy
.
Epilepsy & Behavior
,
28
,
S15
S17
.
Chowdhury
,
F. A.
,
Woldman
,
W.
,
FitzGerald
,
T. H. B.
,
Elwes
,
R. D. C.
,
Nashef
,
L.
,
Terry
,
J. R.
, &
Richardson
,
M. P.
(
2014
).
Revealing a brain network endophenotype in families with idiopathic generalised epilepsy
.
PLOS ONE
,
9
(
10
),
1
8
. https://doi.org/10.1371/journal.pone.0110136
Cirillo
,
E. N. M.
, &
Lebowitz
,
J. L.
(
1998
).
Metastability in the two-dimensional Ising model with free boundary conditions
.
Journal of Statistical Physics
,
90
(
1
),
211
226
. https://doi.org/10.1023/A:1023255802455
Clemens
,
B.
(
2004
).
Pathological theta oscillations in idiopathic generalised epilepsy
.
Clinical Neurophysiology
,
115
(
6
),
1436
1441
. https://doi.org/10.1016/j.clinph.2004.01.018
Clemens
,
B.
,
Puskás
,
S.
,
Besenyei
,
M.
,
Spisák
,
T.
,
Opposits
,
G.
,
Hollódy
,
K.
, …
Emri
,
M.
(
2013
).
Neurophysiology of juvenile myoclonic epilepsy: EEG-based network and graph analysis of the interictal and immediate preictal states
.
Epilepsy Research
,
106
(
3
),
357
369
. https://doi.org/10.1016/j.eplepsyres.2013.06.017
da Silva
,
F. H. L.
,
Blanes
,
W.
,
Kalitzin
,
S. N.
,
Parra
,
J.
,
Suffczynski
,
P.
, &
Velis
,
D. N.
(
2003
).
Dynamical diseases of brain systems: different routes to epileptic seizures
.
IEEE Transactions on Biomedical Engineering
,
50
(
5
),
540
548
. https://doi.org/10.1109/TBME.2003.810703
Deco
,
G.
,
Senden
,
M.
, &
Jirsa
,
V.
(
2012
).
How anatomy shapes dynamics: a semi-analytical study of the brain at rest by a simple spin model
.
Frontiers in Computational Neuroscience
,
6
,
68
. https://doi.org/10.3389/fncom.2012.00068
Delgado-Escueta
,
A. V.
, &
Enrile-Bacsal
,
F.
(
1984
).
Juvenile myoclonic epilepsy of janz
.
Neurology
,
34
(
3
),
285
285
. https://doi.org/10.1212/WNL.34.3.285
Desikan
,
R. S.
,
Ségonne
,
F.
,
Fischl
,
B.
,
Quinn
,
B. T.
,
Dickerson
,
B. C.
,
Blacker
,
D.
, …
Killiany
,
R. J.
(
2006
).
An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest
.
NeuroImage
,
31
(
3
),
968
980
. https://doi.org/10.1016/j.neuroimage.2006.01.021
Destrieux
,
C.
,
Fischl
,
B.
,
Dale
,
A.
, &
Halgren
,
E.
(
2010
).
Automatic parcellation of human cortical gyri and sulci using standard anatomical nomenclature
.
NeuroImage
,
53
(
1
),
1
15
. https://doi.org/10.1016/j.neuroimage.2010.06.010
Elger
,
C. E.
,
Widman
,
G.
,
Andrzejak
,
R.
,
Arnhold
,
J.
,
David
,
P.
, &
Lehnertz
,
K.
(
2000
).
Nonlinear EEG analysis and its potential role in epileptology
.
Epilepsia
,
41
(
s3
),
S34
S38
. https://doi.org/10.1111/j.1528-1157.2000.tb01532.x
Elshahabi
,
A.
,
Klamer
,
S.
,
Sahib
,
A. K.
,
Lerche
,
H.
,
Braun
,
C.
, &
Focke
,
N. K.
(
2015
).
Magnetoencephalography reveals a widespread increase in network connectivity in idiopathic/genetic generalized epilepsy
.
PLOS ONE
,
10
(
9
),
1
16
. https://doi.org/10.1371/journal.pone.0138119
Engel
,
J. P. T.
(
1997
).
Epilepsy: A comprehensive textbook
.
Lippincott-Raven
.
Ezaki
,
T.
,
Sakaki
,
M.
,
Watanabe
,
T.
, &
Masuda
,
N.
(
2018
).
Age-related changes in the ease of dynamical transitions in human brain activity
.
Human Brain Mapping
,
39
(
6
),
2673
2688
. https://doi.org/10.1002/hbm.24033
Ezaki
,
T.
,
Watanabe
,
T.
,
Ohzeki
,
M.
, &
Masuda
,
N.
(
2017
).
Energy landscape analysis of neuroimaging data
.
Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences
,
375
(
2096
),
20160287
. https://doi.org/10.1098/rsta.2016.0287
Fernandez
,
S.
,
Donaire
,
A.
,
Maestro
,
I.
,
Seres
,
E.
,
Setoain
,
X.
,
Bargalló
,
N.
, …
Carreño
,
M.
(
2011
).
Functional neuroimaging in startle epilepsy: Involvement of a mesial frontoparietal network
.
Epilepsia
,
52
(
9
),
1725
1732
. https://doi.org/10.1111/j.1528-1167.2011.03172.x
Goker
,
I.
,
Osman
,
O.
,
Ozekes
,
S.
,
Baslo
,
M. B.
,
Ertas
,
M.
, &
Ulgen
,
Y.
(
2012
).
Classification of juvenile myoclonic epilepsy data acquired through scanning electromyography with machine learning algorithms
.
Journal of Medical Systems
,
36
(
5
),
2705
2711
. https://doi.org/10.1007/s10916-011-9746-6
Gordon
,
E. M.
,
Laumann
,
T. O.
,
Adeyemo
,
B.
,
Huckins
,
J. F.
,
Kelley
,
W. M.
, &
Petersen
,
S. E.
(
2014
).
Generation and evaluation of a cortical area parcellation from resting-state correlations
.
Cerebral Cortex
,
26
(
1
),
288
303
. https://doi.org/10.1093/cercor/bhu239
Gotman
,
J.
,
Grova
,
C.
,
Bagshaw
,
A.
,
Kobayashi
,
E.
,
Aghakhani
,
Y.
, &
Dubeau
,
F.
(
2005
).
Generalized epileptic discharges show thalamocortical activation and suspension of the default state of the brain
.
Proceedings of the National Academy of Sciences
,
102
(
42
),
15236
15240
. https://doi.org/10.1073/pnas.0504935102
Gu
,
S.
,
Cieslak
,
M.
,
Baird
,
B.
,
Muldoon
,
S. F.
,
Grafton
,
S. T.
,
Pasqualetti
,
F.
, &
Bassett
,
D. S.
(
2018
).
The energy landscape of neurophysiological activity implicit in brain network structure
.
Scientific Reports
,
8
(
1
),
2507
. https://doi.org/10.1038/s41598-018-20123-8
Hall
,
E. L.
,
Woolrich
,
M. W.
,
Thomaz
,
C. E.
,
Morris
,
P. G.
, &
Brookes
,
M. J.
(
2013
).
Using variance information in magnetoencephalography measures of functional connectivity
.
NeuroImage
,
67
,
203
212
. https://doi.org/10.1016/j.neuroimage.2012.11.011
Hamandi
,
K.
,
Salek-Haddadi
,
A.
,
Laufs
,
H.
,
Liston
,
A.
,
Friston
,
K.
,
Fish
,
D. R.
, …
Lemieux
,
L.
(
2006
).
EEG-fMRI of idiopathic and secondarily generalized epilepsies
.
NeuroImage
,
31
(
4
),
1700
1710
. https://doi.org/10.1016/j.neuroimage.2006.02.016
Hastings
,
W. K.
(
1970
).
Monte Carlo sampling methods using Markov chains and their applications
.
Biometrika
,
57
(
1
),
97
109
. https://doi.org/10.1093/biomet/57.1.97
Hattingen
,
E.
,
Lückerath
,
C.
,
Pellikan
,
S.
,
Vronski
,
D.
,
Roth
,
C.
,
Knake
,
S.
, …
Pilatus
,
U
.
(
2014
).
Frontal and thalamic changes of gaba concentration indicate dysfunction of thalamofrontal networks in juvenile myoclonic epilepsy
.
Epilepsia
,
55
(
7
),
1030
1037
. https://doi.org/10.1111/epi.12656
Hillebrand
,
A.
,
Tewarie
,
P.
,
van Dellen
,
E.
,
Yu
,
M.
,
Carbo
,
E. W. S.
,
Douw
,
L.
, …
Stam
,
C. J.
(
2016
).
Direction of information flow in large-scale resting-state networks is frequency-dependent
.
Proceedings of the National Academy of Sciences
,
113
(
14
),
3867
3872
. https://doi.org/10.1073/pnas.1515657113
Hipp
,
J. F.
,
Hawellek
,
D. J.
,
Corbetta
,
M.
,
Siegel
,
M.
, &
Engel
,
A. K.
(
2012
).
Large-scale cortical correlation structure of spontaneous oscillatory activity
.
Nature Neuroscience
,
15
,
884
. https://doi.org/10.1038/nn.3101
Horwitz
,
B.
(
2003
).
The elusive concept of brain connectivity
.
NeuroImage
,
19
(
2
),
466
470
. https://doi.org/10.1016/S1053-8119(03)00112-5
Iasemidis
,
L. D.
,
Chris Sackellares
,
J.
,
Zaveri
,
H. P.
, &
Williams
,
W. J.
(
1990
).
Phase space topography and the Lyapunov exponent of electrocorticograms in partial seizures
.
Brain Topography
,
2
(
3
),
187
201
. https://doi.org/10.1007/BF01140588
Jaynes
,
E. T.
(
1957
).
Information theory and statistical mechanics
.
Physical Review. Series II
,
106
(
4
),
620
630
. https://doi.org/10.1103/PhysRev.106.620
Kang
,
J.
,
Pae
,
C.
, &
Park
,
H.-J.
(
2019
).
Graph-theoretical analysis for energy landscape reveals the organization of state transitions in the resting-state human cerebral cortex
.
PLoS One
,
14
(
9
),
1
25
. https://doi.org/10.1371/journal.pone.0222161
Kannathal
,
N.
,
Choo
,
M. L.
,
Acharya
,
U. R.
, &
Sadasivan
,
P.
(
2005
).
Entropies for detection of epilepsy in EEG
.
Computer Methods and Programs in Biomedicine
,
80
(
3
),
187
194
. https://doi.org/10.1016/j.cmpb.2005.06.012
Kelso
,
J. A. S.
(
2012
).
Multistability and metastability: Understanding dynamic coordination in the brain
.
Philosophical transactions of the Royal Society of London. Series B, Biological sciences
,
367
(
1591
),
906
918
. https://doi.org/10.1098/rstb.2011.0351
Kim
,
B.
, &
Oertzen
,
T. v.
(
2018
).
Classifiers as a model-free group comparison test
.
Behavior Research Methods
,
50
(
1
),
416
426
. https://doi.org/10.3758/s13428-017-0880-z
Kim
,
J. H.
,
Lee
,
J. K.
,
Koh
,
S.-B.
,
Lee
,
S.-A.
,
Lee
,
J.-M.
,
Kim
,
S. I.
, &
Kang
,
J. K.
(
2007
).
Regional grey matter abnormalities in juvenile myoclonic epilepsy: A voxel-based morphometry study
.
NeuroImage
,
37
(
4
),
1132
1137
. https://doi.org/10.1016/j.neuroimage.2007.06.025
Kim
,
S.-Y.
,
Hwang
,
Y.-H.
,
Lee
,
H.-W.
,
Suh
,
C.-K.
,
Kwon
,
S.-H.
, &
Park
,
S.-P.
(
2007
).
Cognitive impairment in juvenile myoclonic epilepsy
.
Journal of Clinical Neurology (Seoul, Korea)
,
3
(
2
),
86
92
. https://doi.org/10.3988/jcn.2007.3.2.86
Klein
,
A.
, &
Tourville
,
J.
(
2012
).
101 labeled brain images and a consistent human cortical labeling protocol
.
Frontiers in Neuroscience
,
6
,
171
. https://doi.org/10.3389/fnins.2012.00171
Lee
,
C.
,
Kim
,
S.-M.
,
Jung
,
Y.-J.
,
Im
,
C.-H.
,
Kim
,
D. W.
, &
Jung
,
K.-Y.
(
2014
).
Causal influence of epileptic network during spike-and-wave discharge in juvenile myoclonic epilepsy
.
Epilepsy Research
,
108
(
2
),
257
266
. https://doi.org/10.1016/j.eplepsyres.2013.11.005
Li
,
Q.
,
Cao
,
W.
,
Liao
,
X.
,
Chen
,
Z.
,
Yang
,
T.
,
Gong
,
Q.
, …
Yao
,
D.
(
2015
).
Altered resting state functional network connectivity in children absence epilepsy
.
Journal of the Neurological Sciences
,
354
(
1
),
79
85
. https://doi.org/10.1016/j.jns.2015.04.054
Liao
,
W.
,
Zhang
,
Z.
,
Pan
,
Z.
,
Mantini
,
D.
,
Ding
,
J.
,
Duan
,
X.
, …
Chen
,
H.
(
2010
).
Altered functional connectivity and small-world in mesial temporal lobe epilepsy
.
PLoS One
,
5
(
1
),
1
11
. https://doi.org/10.1371/journal.pone.0008525
Liao
,
W.
,
Zhang
,
Z.
,
Pan
,
Z.
,
Mantini
,
D.
,
Ding
,
J.
,
Duan
,
X.
, …
Chen
,
H.
(
2011
).
Default mode network abnormalities in mesial temporal lobe epilepsy: a study combining fMRI and DTI
.
Human Brain Mapping
,
32
(
6
),
883
895
. https://doi.org/10.1002/hbm.21076
Macke
,
J. H.
,
Murray
,
I.
, &
Latham
,
P. E.
(
2012
).
How biased are maximum entropy models?
Advances in Neural Information Processing Systems
,
1
9
.
McGill
,
M. L.
,
Devinsky
,
O.
,
Kelly
,
C.
,
Milham
,
M.
,
Castellanos
,
F. X.
,
Quinn
,
B. T.
, …
Thesen
,
T.
(
2012
).
Default mode network abnormalities in idiopathic generalized epilepsy
.
Epilepsy & Behavior
,
23
(
3
),
353
359
. https://doi.org/10.1016/j.yebeh.2012.01.013
Michel
,
C. M.
, &
Koenig
,
T.
(
2018
).
EEG microstates as a tool for studying the temporal dynamics of whole-brain neuronal networks: a review
.
NeuroImage
,
180
,
577
593
. https://doi.org/10.1016/j.neuroimage.2017.11.062
Muthukumaraswamy
,
S. D.
,
Carhart-Harris
,
R. L.
,
Moran
,
R. J.
,
Brookes
,
M. J.
,
Williams
,
T. M.
,
Errtizoe
,
D.
, …
Nutt
,
D. J.
(
2013
).
Broadband cortical desynchronization underlies the human psychedelic state
.
Journal of Neuroscience
,
33
(
38
),
15171
15183
. https://doi.org/10.1523/JNEUROSCI.2063-13.2013
Niedermeyer
,
E.
(
2005
).
Electroencephalography
.
Lippincott Williams & Wilkins
.
Niso
,
G.
,
Carrasco
,
S.
,
Gudín
,
M.
,
Maestú
,
F.
,
del Pozo
,
F.
, &
Pereda
,
E.
(
2015
).
What graph theory actually tells us about resting state interictal MEG epileptic activity
.
NeuroImage: Clinical
,
8
,
503
515
. https://doi.org/10.1016/j.nicl.2015.05.008
Nolte
,
G.
(
2003
).
The magnetic lead field theorem in the quasi-static approximation and its use for magnetoencephalography forward calculation in realistic volume conductors
.
Physics in Medicine and Biology
,
48
(
22
),
3637
3652
.
O’Muircheartaigh
,
J.
,
Vollmar
,
C.
,
Barker
,
G.
,
Kumari
,
V.
,
Symms
,
M.
,
Thompson
,
P.
, …
Richardson
,
M.
(
2011
).
Focal structural changes and cognitive dysfunction in juvenile myoclonic epilepsy
.
Neurology
,
76
(
1
),
34
40
. https://doi.org/10.1212/WNL.0b013e318203e93d
Papanicolaou
,
A.
,
Pazo-Alvarez
,
P.
,
Castillo
,
E.
,
Billingsley-Marshall
,
R.
,
Breier
,
J.
,
Swank
,
P.
, …
Passaro
,
A.
(
2006
).
Functional neuroimaging with MEG: Normative language profiles
.
NeuroImage
,
33
(
1
),
326
342
. https://doi.org/10.1016/j.neuroimage.2006.06.020
Piazzini
,
A.
,
Turner
,
K.
,
Vignoli
,
A.
,
Canger
,
R.
, &
Canevini
,
M. P.
(
2008
).
Frontal cognitive dysfunction in juvenile myoclonic epilepsy
.
Epilepsia
,
49
(
4
),
657
662
. https://doi.org/10.1111/j.1528-1167.2007.01482.x
Quinn
,
A. J.
,
Vidaurre
,
D.
,
Abeysuriya
,
R.
,
Becker
,
R.
,
Nobre
,
A. C.
, &
Woolrich
,
M. W.
(
2018
).
Task-evoked dynamic network analysis through hidden markov modeling
.
Frontiers in Neuroscience
,
12
,
603
. https://doi.org/10.3389/fnins.2018.00603
Rolls
,
E. T.
,
Joliot
,
M.
, &
Tzourio-Mazoyer
,
N.
(
2015
).
Implementation of a new parcellation of the orbitofrontal cortex in the automated anatomical labeling atlas
.
NeuroImage
,
122
,
1
5
. https://doi.org/10.1016/j.neuroimage.2015.07.075
Ronan
,
L.
,
Alhusaini
,
S.
,
Scanlon
,
C.
,
Doherty
,
C. P.
,
Delanty
,
N.
, &
Fitzsimons
,
M.
(
2012
).
Widespread cortical morphologic changes in juvenile myoclonic epilepsy: Evidence from structural MRI
.
Epilepsia
,
53
(
4
),
651
658
. https://doi.org/10.1111/j.1528-1167.2012.03413.x
Rosazza
,
C.
, &
Minati
,
L.
(
2011
).
Resting-state brain networks: literature review and clinical applications
.
Neurological Sciences
,
32
(
5
),
773
785
. https://doi.org/10.1007/s10072-011-0636-y
Routley
,
B. C.
,
Singh
,
K. D.
,
Hamandi
,
K.
, &
Muthukumaraswamy
,
S. D.
(
2017
).
The effects of AMPA receptor blockade on resting magnetoencephalography recordings
.
Journal of Psychopharmacology
,
31
(
12
),
1527
1536
. https://doi.org/10.1177/0269881117736915
Salvador
,
R.
,
Suckling
,
J.
,
Coleman
,
M. R.
,
Pickard
,
J. D.
,
Menon
,
D.
, &
Bullmore
,
E.
(
2005
).
Neurophysiological architecture of functional magnetic resonance images of human brain
.
Cerebral Cortex
,
15
(
9
),
1332
1342
. https://doi.org/10.1093/cercor/bhi016
Sarnthein
,
J.
,
Morel
,
A.
,
Von Stein
,
A.
, &
Jeanmonod
,
D.
(
2003
).
Thalamic theta field potentials and EEG: high thalamocortical coherence in patients with neurogenic pain, epilepsy and movement disorders
.
Thalamus & Related Systems
,
2
(
3
),
231
238
.
Savic
,
I.
,
Lekvall
,
A.
,
Greitz
,
D.
, &
Helms
,
G.
(
2000
).
Mr spectroscopy shows reduced frontal lobe concentrations of n-acetyl aspartate in patients with juvenile myoclonic epilepsy
.
Epilepsia
,
41
(
3
),
290
296
. https://doi.org/10.1111/j.1528-1157.2000.tb00158.x
Schneidman
,
E.
,
Berry
,
M. J.
,
Segev
,
R.
, &
Bialek
,
W.
(
2006
).
Weak pairwise correlations imply strongly correlated network states in a neural population
.
Nature
,
440
(
7087
),
1007
1012
. https://doi.org/10.1038/nature04701
Shmueli
,
G.
(
2010
).
To explain or to predict?
Statistical Science
,
25
(
3
),
289
310
. https://doi.org/10.1214/10-STS330
Song
,
Y.
,
Crowcroft
,
J.
, &
Zhang
,
J.
(
2012
).
Automatic epileptic seizure detection in EEGs based on optimized sample entropy and extreme learning machine
.
Journal of Neuroscience Methods
,
210
(
2
),
132
146
. https://doi.org/10.1016/j.jneumeth.2012.07.003
Song
,
Y.
, &
Zhang
,
J.
(
2013
).
Automatic recognition of epileptic EEG patterns via extreme learning machine and multiresolution feature extraction
.
Expert Systems with Applications
,
40
(
14
),
5477
5489
. https://doi.org/10.1016/j.eswa.2013.04.025
Stam
,
C.
(
2005
).
Nonlinear dynamical analysis of EEG and MEG: Review of an emerging field
.
Clinical Neurophysiology
,
116
(
10
),
2266
2301
. https://doi.org/10.1016/j.clinph.2005.06.011
Stam
,
C. J.
, &
Straaten
,
E. C. W. v.
(
2012
).
The organization of physiological brain networks
.
Clinical Neurophysiology
,
123
(
6
),
1067
1087
. https://doi.org/10.1016/j.clinph.2012.01.011
Swartz
,
B.
,
Halgren
,
E.
,
Simpkins
,
F.
, &
Syndulko
,
K.
(
1994
).
Primary memory in patients with frontal and primary generalized epilepsy
.
Journal of Epilepsy
,
7
(
3
),
232
241
. https://doi.org/10.1016/0896-6974(94)90034-5
Tang
,
A.
,
Jackson
,
D.
,
Hobbs
,
J.
,
Chen
,
W.
,
Smith
,
J. L.
,
Patel
,
H.
, …
Beggs
,
J. M.
(
2008
).
A maximum entropy model applied to spatial and temporal correlations from cortical networks in vitro
.
Journal of Neuroscience
,
28
(
2
),
505
518
. https://doi.org/10.1523/JNEUROSCI.3359-07.2008
Tewarie
,
P.
,
Schoonheim
,
M. M.
,
Stam
,
C. J.
,
van der Meer
,
M. L.
,
van Dijk
,
B. W.
,
Barkhof
,
F.
, …
Hillebrand
,
A.
(
2013a
).
Cognitive and clinical dysfunction, altered MEG resting-state networks and thalamic atrophy in multiple sclerosis
.
PLoS One
,
8
(
7
),
1
11
. https://doi.org/10.1371/journal.pone.0069318
Tkacik
,
G.
,
Schneidman
,
E.
,
Berry
,
I.
,
Michael
,
J.
, &
Bialek
,
W.
(
2006
).
Ising models for networks of real neurons
.
arXiv preprint q-bio/0611072
.
Trenité
,
D. G. A. K.-N.
,
Schmitz
,
B.
,
Janz
,
D.
,
Delgado-Escueta
,
A. V.
,
Thomas
,
P.
,
Hirsch
,
E.
, …
Genton
,
P.
(
2013
).
Consensus on diagnosis and management of JME: From founder’s observations to current trends
.
Epilepsy & Behavior
,
28
,
S87
S90
. https://doi.org/10.1016/j.yebeh.2012.11.051
Tzourio-Mazoyer
,
N.
,
Landeau
,
B.
,
Papathanassiou
,
D.
,
Crivello
,
F.
,
Etard
,
O.
,
Delcroix
,
N.
, …
Joliot
,
M.
(
2002
).
Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain
.
NeuroImage
,
15
(
1
),
273
289
. https://doi.org/10.1006/nimg.2001.0978
Vollmar
,
C.
,
O’Muircheartaigh
,
J.
,
Symms
,
M.
,
Barker
,
G.
,
Thompson
,
P.
,
Kumari
,
V.
, …
Koepp
,
M.
(
2012
).
Altered microstructural connectivity in juvenile myoclonic epilepsy
.
Neurology
,
78
(
20
),
1555
1559
. https://doi.org/10.1212/WNL.0b013e3182563b44
Vollmar
,
C.
,
O’Muircheartaigh
,
J.
,
Barker
,
G. J.
,
Symms
,
M. R.
,
Thompson
,
P.
,
Kumari
,
V.
, …
Koepp
,
M. J.
(
2011
).
Motor system hyperconnectivity in juvenile myoclonic epilepsy: a cognitive functional magnetic resonance imaging study
.
Brain
,
134
(
6
),
1710
1719
. https://doi.org/10.1093/brain/awr098
Vrba
,
J.
, &
Robinson
,
S. E.
(
2001
).
Signal processing in magnetoencephalography
.
Methods
,
25
(
2
),
249
271
. https://doi.org/10.1006/meth.2001.1238
Watanabe
,
T.
,
Hirose
,
S.
,
Wada
,
H.
,
Imai
,
Y.
,
Machida
,
T.
,
Shirouzu
,
I.
, …
Masuda
,
N.
(
2013
).
A pairwise maximum entropy model accurately describes resting-state human brain networks
.
Nature Communications
,
4
,
1370
. https://doi.org/10.1038/ncomms2388
Watanabe
,
T.
,
Hirose
,
S.
,
Wada
,
H.
,
Imai
,
Y.
,
Machida
,
T.
,
Shirouzu
,
I.
, …
Rees
,
G.
(
2014a
).
Energy landscapes of resting-state brain networks
.
Frontiers in Neuroinformatics
,
8
,
12
. https://doi.org/10.3389/fninf.2014.00012
Watanabe
,
T.
,
Masuda
,
N.
,
Megumi
,
F.
,
Kanai
,
R.
, &
Rees
,
G.
(
2014b
).
Energy landscape and dynamics of brain activity during human bistable perception
.
Nature Communications
,
5
(
1
),
4765
. https://doi.org/10.1038/ncomms5765
Wolf
,
P.
, &
Beniczky
,
S.
(
2014
).
Understanding ictogenesis in generalized epilepsies
.
Expert Review of Neurotherapeutics
,
14
(
7
),
787
798
. https://doi.org/10.1586/14737175.2014.925803
Wolf
,
P.
,
Yacubian
,
E. M. T.
,
Avanzini
,
G.
,
Sander
,
T.
,
Schmitz
,
B.
,
Wandschneider
,
B.
, &
Koepp
,
M.
(
2015a
).
Juvenile myoclonic epilepsy: A system disorder of the brain
.
Epilepsy Research
,
114
,
2
12
. https://doi.org/10.1016/j.eplepsyres.2015.04.008
Wolf
,
P.
,
Yacubian
,
E. M. T.
,
Avanzini
,
G.
,
Sander
,
T.
,
Schmitz
,
B.
,
Wandschneider
,
B.
, &
Koepp
,
M.
(
2015b
).
Juvenile myoclonic epilepsy: A system disorder of the brain
.
Epilepsy Research
,
114
,
2
12
. https://doi.org/10.1016/j.eplepsyres.2015.04.008
Yacubian
,
E. M.
, &
Wolf
,
P.
(
2014
).
Praxis induction. definition, relation to epilepsy syndromes, nosological and prognostic significance. a focused review
.
Seizure
,
23
(
4
),
247
251
. https://doi.org/10.1016/j.seizure.2014.01.011
Yeh
,
F.-C.
,
Tang
,
A.
,
Hobbs
,
J. P.
,
Hottowy
,
P.
,
Dabrowski
,
W.
,
Sher
,
A.
, …
Beggs
,
J. M.
(
2010
).
Maximum entropy approaches to living neural networks
.
Entropy
,
12
(
1
),
89
106
. https://doi.org/10.3390/e12010089
Yu
,
S.
,
Yang
,
H.
,
Nakahara
,
H.
,
Santos
,
G. S.
,
Nikolić
,
D.
, &
Plenz
,
D.
(
2011
).
Higher-order interactions characterized in cortical activity
.
Journal of Neuroscience
,
31
(
48
),
17514
17526
. https://doi.org/10.1523/JNEUROSCI.3127-11.2011
Zamarian
,
L.
,
Höfler
,
J.
,
Kuchukhidze
,
G.
,
Delazer
,
M.
,
Bonatti
,
E.
,
Kemmler
,
G.
, &
Trinka
,
E.
(
2013
).
Decision making in juvenile myoclonic epilepsy
.
Journal of Neurology
,
260
(
3
),
839
846
. https://doi.org/10.1007/s00415-012-6715-z
Zhang
,
L.
,
Li
,
H.
,
Hong
,
P.
, &
Zou
,
X.
(
2016
).
Proton magnetic resonance spectroscopy in juvenile myoclonic epilepsy: A systematic review and meta-analysis
.
Epilepsy Research
,
121
,
33
38
. https://doi.org/10.1016/j.eplepsyres.2016.01.004
Zhu
,
C.
,
Han
,
S.
,
Mao
,
H.
, &
Dally
,
W. J.
(
2016
).
Trained ternary quantization
.
arXiv preprint arXiv:1612.01064
.

Competing Interests

Competing Interests: The authors have declared that no competing interests exist.

Author notes

Handling Editor: Daniele Marinazzo

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 the original work is properly cited. For a full description of the license, please visit https://creativecommons.org/licenses/by/4.0/legalcode.