MEG energy landscape abnormalities in juvenile myoclonic epilepsy

Juvenile myoclonic epilepsy (JME) is a network disorder affecting brain activity and connectivity. However, it is unclear whether JME leads to widespread abnormalities in the network dynamics across different functional networks. Here, we used a pairwise maximum entropy model (pMEM) and energy landscape analysis to characterize network dynamics in MEG resting-state data and its abnormalities in JME. We fitted the pMEM to the MEG oscillatory power in three functional networks: the default mode network (DMN), the frontoparietal network (FPN) and the sensorimotor network (SMN). The pMEM provided an accurate fit to the MEG data in both patient and control groups. We then used pMEM-derived energy values to depict an energy landscape of each network, with a higher energy state corresponding to a lower occurrence probability. JME patients exhibited a lower number of local energy minima than controls, and had elevated energy values in the theta, beta and gamma-band of FPN oscillatory activity as well as the beta-band DMN activity, but not in the SMN. Our findings suggested that JME patients had impaired multi-stability in selective functional networks and frequency bands in the frontoparietal cortices.


Introduction
Juvenile myoclonic epilepsy (JME) has been recognised as a network disorder affecting brain activity and connectivity that leads to cognitive impairments (Wolf et al., 2015). Electrophysiological data suggests that JME has an impact on multiple functional networks, including the fronto-parietal network (FPN), the default mode network (DMN), and the sensorimotor network (SMN) (Clemens et al., 2013). However, it is yet unclear whether JME patients had widespread or selective abnormalities in the network dynamics across functional networks, and whether the abnormalities are frequency specific.
We addressed these questions by proposing a pairwise maximum entropy model (pMEM) approach (Yeh et al., 2010) to analyse MEG resting-sate oscillatory activity. PMEM is a statistical model of the occurrence probability of network states, which has previously been applied to the collective behaviour of spiking neural networks (Schneidman et al., 2006) and BOLD fMRI responses (Ashourvan et al., 2017).
Here, We extended this theoretical framework to MEG oscillatory activity. Based on the fitted pMEM to individual participants, we depicted an energy landscape for each of the networks at theta (4-7 Hz), alpha (8-13 Hz), beta (15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25) and gamma (30-60 Hz) bands. The energy landscape is a graph representation of energy values from all network states, with the graph's connectivity described by an adjacency matrix of regional activation states (Ezaki et al., 2017). We observed that pMEM provided a good fit to the statistical regularities of functional networks in both JME and control groups. Furthermore, JME patients showed reduced numbers of local energy minima and elevated energy values in the theta, beta and gamma-band FPN activity and beta-band DMN activity, but not in the SMN. These findings suggested anatomicallyand frequency-specific network abnormalities in JME.

Methods
Participants 26 JME patients were recruited from a specialist clinic for epilepsy at University Hospital of Wales in Cardiff. 26 agematched healthy control participants were recruited from a regional volunteer panel. 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.

MEG source localization of oscillatory activity
5-minutes whole-head resting-state MEG recordings were made using a 275-channel CTF radial gradiometer system at a sampling rate of 600 Hz. Whole-brain T1-weighted MRI data were acquired using a GE HDx 3T MRI scanner.
For each participant, the structural MRI scan was coregistered to the MEG sensor space. The pre-processed MEG 845 This work is licensed under the Creative Commons Attribution 3.0 Unported License. To view a copy of this license, visit http://creativecommons.org/licenses/by/3.0 data was band-pass filtered into four frequency bands: theta 4-8 Hz, alpha 8-12 Hz, beta 13-30 Hz, and low-gamma 35-60 Hz. For each frequency band, we computed the inverse source reconstruction using LCMV beamformer. The atlasbased source reconstruction was used to derive virtual sensors for every voxel in each of the 90 regions of the Automated Anatomical Label (AAL) atlas.
We focused our analysis on three networks: FPN, DMN and SMN. Each network comprised of bilateral regions of interest (ROIs) from the AAL atlas ( Figure 1). For each ROI, its representative time course was obtained from the voxel in that ROI with the highest temporal standard-deviation. 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 in each frequency band.

PMEM of MEG oscillatory activity
We fitted a pMEM to individual participant's MEG data, separately for each network and each frequency band. According to the principle of maximum entropy, among all probabilistic models, one should choose the one with the largest uncertainty (i.e., entropy), because it makes minimum assumption of additional information that would otherwise lower the uncertainty (Yeh et al., 2010).
Consider a network consisting of N ROIs. 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 network state at time t is thus given by a N-dimensional binary vector s(t) = [r 1 (t), r 2 (t), ..., r N (t)].
The N-ROI network has a total of 2 N possible states s k (k = 1, ..., 2 N ). From the binarized MEG oscillatory activity, we calculated the probability of occurrence of each network state, denoted by P emp (s k ). We further estimated the empirical average activation rate for each ROI r i emp and all pairwise correlations between any two ROIs r i r j emp : where T denotes the total number of time points in the data.
The fitting procedure aimed to identify a pMEM model that preserves the constraints in Equation 1 and reproduces the empirical state probability P emp (s k ) with the maximum entropy.
It is known that the pMEM follows the Boltzman distribution: and E(s k ) represents the energy of the network state s k , by: where r i (s k ) refers to the binary value of r i for the network state s k . h and J are the model parameters to be estimated from the data. h = [h 1 , h 2 , ..., h N ] represents the bias of high oscillatory activity in each ROI, and J = [J 1,1 , J 1,2 , ..., J N,N ] represents the coupling strength between two ROIs. The model predictions of the average activation rate r i mod and pairwise correlations r i r j mod are given by: We used a gradient ascent algorithm to iteratively update h and J, until r i mod and r i r j mod from the pMEM match the constraints r i emp and r i r j emp . In each iteration, the updates of the parameters were given by h new i = h old i + ε( r i emp − r i mod ) and J new i j = J old i j +ε( r i r j emp − r i r j mod ). The learning rate ε was set to 10 −8 .

Energy landscape of resting-state networks
The pMEM defines the energy E(s k ) of every network state s k , and its value indicates the model prediction of the inverse appearance probability of s k . 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 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 , characterising state probabilities and state transitions from the perspective of attractor dynamics (Watanabe et al., 2014). 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 connected, if and only if one ROI of the network had different binarized oscillatory activity (high vs. low).

Quantitative measures of energy landscape
We used two measures to understand the differences in the energy landscape between JME patients and healthy controls.
1. The number of local energy minima A local energy minimum was defined as the network state with a lower energy value than all its adjacent neighbouring states (Watanabe et al., 2014). Because lower energy corresponds to a higher probability of occurrence, network states of local minima can be conceptualised as attractors in attractor dynamics.

The relative energy of the local minima
We calculated the mean energy of each state s k across all participants.
The mean energy was used to depict an aggregated landscape, allowing to identify common energy minima shared between JME patients and controls. Because the shape of a energy landscape was partly determined by the global minimum (Watanabe et al., 2014), for each participant, we calculated the difference between a significant local minimum and the global energy minimum (i.e., the lowest energy value). We then compared this relative energy of the local minima between JME patients and controls. From the networks with significant alternations of relative energy values in JME patients, we constructed a disconnectivity graph to describe clusters of energy minima and the relationships between them (Ezaki et al., 2017), where the clusters represent groups of local minima with higher probabilities of co-occurrence.

Fitting of pMEM to MEG oscillatory activity
The empirical occurrence probabilities of possible network states closely matched the model predictions from fitted pMEM (R 2 > 0.8 in all networks). We further used an accuracy index to quantify the goodness of fit, which was calculated as the percentage of improvement of the pMEM fit to the empirical data compared with a null model. The null model assumed no pairwise correlation between ROIs (i.e., an independent maximum entropy model). The accuracy indexes were high (> 80%) and did not differ significantly between groups (F(1, 50) = 2.71, p = 0.11), suggesting the robustness of pMEM in both patients and controls.
Inferences from pMEM energy landscape 1. The number of energy minima We located local minima on the pMEM-derived energy landscape. Because a local minimum state has a higher occurrence probability than all of its neighbouring states, transitions of network states near an energy minimum is akin to a fixed point attractor in a dynamic system, and the number of energy minima quantifies the degree of multi-stability of a network. JME patients had significantly less local energy minima than controls (F(1, 50) = 5.27, p = 0.026). No significant network by group (F(1.86, 93.27) = 2.06, p = 0.14) or frequency band by group (F(2.82, 144.17) = 1.34, p = 0.26) interaction was observed.

The relative energy values of the local minima
In the FPN (Figure 2), the relative energy values at the lo-cal minima were significantly higher in JME patients than in controls in the theta-band (F(1, 50) = 18.90, p < 0.0001), beta-band (F(1, 50) = 15.43, p = 0.0002), and gamma band (F(1, 50) = 7.2558, p = 0.009), but not in the alpha band. In the DMN, the relative energy values were significantly higher in JME patients than in controls in the beta-band (F(1, 50) = 13.72, p = 0.0005). In the SMN, there was no significant group difference in the relative energy values.

Discussion
We found that patients with JME showed altered pMEMderived energy landscapes in selective resting-state networks and frequency bands. For the energy landscapes estimated at the individual level, JME patients exhibited lower numbers of local minima than controls. 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) and DMN (beta band) oscillatory activities, but not the SMN. Our results confirmed the abnormalities of electrophysiological signals in JME, and provided new insights into JME pathophysiology affecting selective functional network configurations.
The fewer number of local minima and elevated energy values in JME suggested abnormalities in the multi-stable dynamics properties of the brain networks that may be prone to perturbation and ictogenesis. The energy landscape further allowed to characterise clusters of energy minima and their hierarchies from the disconnectivity graphs. 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, all the DMN energy minima contained co-activation in bilateral ROIs, consistent with the evidence of strong interhemispheric and longrange connectivity in the DMN during awake.
Interestingly, the MEG energy landscape measures of the SMN did not differ significantly between JME patients and controls. This might seem counter-intuitive, giving that motor cortex hyperexcitability has been reported in JME (Badawy et al., 2006). Nevertheless, previous research on resting-state functional connectivity also showed the lack of altered connectivity in the motor cortex in JME (Liao et al., 2011). Our results suggested that the network states (i.e., patterns of co-activation) in the SMN were not affected by JME during rest, this does not rule out the possibility of network dysfunction in the motor circuit under stimulation or perturbation.
Our study provides new methods for studying the dynamics of MEG oscillatory activity. The pMEM has previously been applied to quantify the dynamics of BOLD fMRI data (Watanabe et al., 2014). However, achieving satisfactory pMEM fitting requires a large number of data samples. Because of the low temporal resolution of the BOLD signal, the applications of pMEM for fMRI often need long scanning time that Figure 2: Activity patterns of significant conditions: theta, beta and gamma FPN and beta DMN. Top panel shows disconnectivity graph of an energy barrier between pairs of local minima basins. Middle panel depicts anatomical location of the activity pattern displayed by each local minimum. Grey box denotes low oscillatory activity (-1) and white box denotes high activity (+1). Bottom panel represents T-value of pairwise t-test comparison between normalised energy values from JME patients and controls.
Asterisk above the bar means that the comparison was significant with p < 0.05 (Sidak corrected).
may not be practical for clinical populations, or to concatenate data across participants that limits the possibility of individuallevel inferences. Here, we demonstrated that the high sampling rate afforded by source-localised MEG data suited well for the pMEM analysis, providing anatomically-specific and frequency-dependent results. It is possible to apply our approach to investigate rapid changes of network energy landscapes during active tasks.
In conclusion, JME patients exhibited abnormal energy landscapes of MEG oscillatory activity in selective brain networks and frequency bands, with less local energy minimal and elevated energy levels leading to compromised multistable network dynamics. These results have the potential to be exploited in future diagnostic and pharmacological studies for a mechanistic understanding of ictogenesis in JME.