Energy landscape of resting magnetoencephalography reveals fronto-parietal network impairments in epilepsy

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.

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.

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.

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).

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.

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).

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.

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 r_i(t), with r_i(t) = +1 for high oscillatory activity and r_i(t) = −1 for low oscillatory activity. The activity pattern of a N-dimensional binary vector s(t) = [r₁(t), r₂(t), …, r_N(t)] represents the state of the network at time t.

We used a gradient ascent algorithm to iteratively update h and J, until 〈r_i〉_mod and 〈r_ir_j〉_mod match 〈r_i〉_emp and 〈r_ir_j〉_emp from the observed data, with a stop criterion of 5 × 10⁶ steps. In each iteration, the updates of the parameters were given by $hinew$ = $hiold$ + ϵ(〈r_i〉_emp − 〈r_i〉_mod) and $Jijnew$ = $Jijold$ + ϵ(〈r_ir_j〉_emp − 〈r_ir_j〉_mod). The learning rate ϵ was set to 10⁻⁸.

The pMEM parameters h and J determine the energy E(s_k) of each network state s_k(k = 1, …, 2^N), 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(s_k) 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(s_k) indicates the model prediction of the inverse appearance probability of the state s_k under the empirical constraints of regional activity (parameter h) and regional interactions (parameter J). For instance, if E(s_i) < E(s_j), the pMEM predicts that the network activity pattern is more likely to be at the state s_i than s_j.

For each resting-state network and each frequency band, we depicted an energy landscape as a graph of the energy function across the 2^N possible network states s_k, 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(s_k) 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.

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.

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 2^N 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).

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 Ē(s_k) of each network state s_k 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 Ē(s_k) 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).

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 s_k. On each time step, one of the current state’s N neighboring state s_k′ 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(s_k) − E(s_k′)] when E(s_k′) > E(s_k) 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.

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.

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.

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, χ² 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).

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 2^N possible states. The pMEM provides a statistical model of the occurrence probabilities of the 2^N 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 2^N 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 (R² > 0.90 in all networks and frequency bands, based on a log-log regression, Figure 3) and control groups (R² > 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 (χ² = 87.5, p < 0.00001) and the frequency bands (χ² = 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.

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.

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.

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.

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.

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).

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).

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.

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.

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.

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.

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.