Accumulated source imaging of brain activity with both low and high-frequency neuromagnetic signals

Recent studies have revealed the importance of high-frequency brain signals (>70 Hz). One challenge of high-frequency signal analysis is that the size of time-frequency representation of high-frequency brain signals could be larger than 1 terabytes (TB), which is beyond the upper limits of a typical computer workstation's memory (<196 GB). The aim of the present study is to develop a new method to provide greater sensitivity in detecting high-frequency magnetoencephalography (MEG) signals in a single automated and versatile interface, rather than the more traditional, time-intensive visual inspection methods, which may take up to several days. To address the aim, we developed a new method, accumulated source imaging, defined as the volumetric summation of source activity over a period of time. This method analyzes signals in both low- (1~70 Hz) and high-frequency (70~200 Hz) ranges at source levels. To extract meaningful information from MEG signals at sensor space, the signals were decomposed to channel-cross-channel matrix (CxC) representing the spatiotemporal patterns of every possible sensor-pair. A new algorithm was developed and tested by calculating the optimal CxC and source location-orientation weights for volumetric source imaging, thereby minimizing multi-source interference and reducing computational cost. The new method was implemented in C/C++ and tested with MEG data recorded from clinical epilepsy patients. The results of experimental data demonstrated that accumulated source imaging could effectively summarize and visualize MEG recordings within 12.7 h by using approximately 10 GB of computer memory. In contrast to the conventional method of visually identifying multi-frequency epileptic activities that traditionally took 2–3 days and used 1–2 TB storage, the new approach can quantify epileptic abnormalities in both low- and high-frequency ranges at source levels, using much less time and computer memory.

Introduction

Recent studies have revealed the significance of high-frequency brain signals – such as high-frequency oscillations (HFOs, 90–200 Hz), ripples (80–250 Hz) and fast ripples (250–500 Hz) relative to the conventional lower frequency brain signals (<70 Hz) (Pulvermuller et al., 1997; Guggisberg et al., 2007; Gotman, 2010; Worrell et al., 2012). One of the important motivations behind the study of high-frequency brain signals is their potential clinical applications. HFOs may be important biomarkers of epileptogenicity, a revolutionary finding revealed in recent years (Xiang et al., 2004, 2009a, 2010). Clinical data have revealed that removal of HFO-generating areas lead to improved surgical outcomes (Haegelen et al., 2013). In addition, by using HFOs, it is possible to substantially reduce the extent of cortical resections in epilepsy surgery procedures without compromising seizure control (Weiss et al., 2013). Furthermore, HFOs also play a very important role in many brain disorders (Uhlhaas et al., 2011). For example, schizophrenia is associated with abnormal amplitude and synchrony of high frequency activities (Uhlhaas and Singer, 2013). Of note, the study of high-frequency brain signals may shed light on some of the fundamental mechanisms of neuronal functions and brain disorders.

Numerous challenges exist in the study of high-frequency brain signals with magnetoencephalography (MEG) and electroencephalography (EEG) (Xiang et al., 2004, 2010, 2013; Dalal et al., 2008; Papadelis et al., 2009; Chen et al., 2010; Gotman, 2010; Gummadavelli et al., 2013). First, the size of high sampling rate data can be over 12 terabytes (TB) (Blanco et al., 2011). The size of high sampling rate data can cause a substantial amount of data, posing a challenge for data transfer, storage, archiving, sharing and analysis (Van Essen et al., 2012; Worrell et al., 2012; Zafeiriou and Vargiami, 2012; Zijlmans et al., 2012b). Given the massive amounts of high-sampling rate MEG/EEG data that are collected from patients and research subjects, it is impractical to rely on a visual review of HFOs (Haegelen et al., 2013; Tort et al., 2013; Xiang et al., 2013). Second, in clinical practice, MEG/EEG data are typically analyzed with other neuroimaging data such as invasive recordings, magnetic resonance imaging (MRI) and functional MRI (fMRI). The considerable volume of multi-modal neuroimaging data produced across different communities has posed a daunting challenge to the traditional methods of data sharing, data archiving, data processing, and data interpreting (Van Essen et al., 2012; Worrell et al., 2012; Zafeiriou and Vargiami, 2012; Zijlmans et al., 2012b). Though the multi-modal data enhance our collective understanding of the structure and function of the brain, it is a challenge to handle these varied and heterogeneous datasets. Even with modern computational innovations, there remain technical challenges in data transfer, storage, and analysis of large data sets of more than 12 TB (Brinkmann et al., 2009; Le Van Quyen et al., 2010). Third, the best way to clinically utilize analysis from high-frequency brain signals remains a challenge. While it has been demonstrated that the brain generates signals in wide frequency ranges, there are currently no established criteria for distinguishing physiologic high-frequency signals from pathologic neuromagnetic signals (Worrell et al., 2012; Zijlmans et al., 2012b; Haegelen et al., 2013; Matsumoto et al., 2013; Pail et al., 2013; Srejic et al., 2013; Tort et al., 2013). Although multiple studies with invasive recordings have shown the feasibility and potential clinical importance of detecting HFOs (Jirsch et al., 2006; Engel et al., 2009; Jacobs et al., 2009, 2012; Levesque et al., 2011; Andrade-Valenca et al., 2012; Dumpelmann et al., 2012; Zijlmans et al., 2012b), there is no noninvasive method which can be used for clinical purposes. One remaining important clinical question is whether a noninvasive method can extract and visualize meaningful HFOs from the brain for research and clinical purposes.

This study aimed to resolve the aforementioned challenges associated with large scale high-frequency signal processing by developing novel analysis methodologies and workflows for the MEG data. Since the computer memory limits for a 32 bit and 64 bit operating system are 4 GB and 192 GB (Windows 7, respectively) (http://msdn.microsoft.com/en-us/library/windows/desktop/aa366778(v=vs.85).aspx) and the size of high-frequency brain signals are usually larger than 12 TB (Blanco et al., 2011), one methodological question this study would like to address is whether new algorithms could minimize the use of computer memory and storage. To solve the challenges of analyzing more than 12 TB of both high and low frequency MEG data, we mathematically and experimentally developed a systematic approach to extract meaningful frequency specific and spatiotemporal information from MEG data. Accumulated spectrograms, a technique which maximizes the signal power of the frequency of interest while simultaneously minimizing other frequency contents, provides a novel method of quantifying and visualizing the frequency signatures of brain activity in both low- and high-frequency ranges. Accumulated source imaging, which volumetrically reconstructs source activity in multiple frequency ranges, provides source images for clinicians to analyze epileptic activity at source levels. The central hypothesis of our research is that neuromagnetic brain signals in both low and high frequency ranges could be localized and visualized with accumulated source imaging. The new algorithm calculated optimal channel-cross-channel (CxC) matrices and source location-orientation weights for volumetric source imaging, minimizing multi-source interference and reducing computational cost. To demonstrate the advancements of the new methods in research and clinical settings, MEG data from subjects were obtained, analyzed, and demonstrated in 2D and 3D environments.

Materials and Methods

Detection of Low- and High-Frequency MEG Signals at Sensor Levels

Multi-channel MEG data had to be digitized at a high sampling rate because the sampling rate must be at least two times higher than the frequency edge of interest. For the analysis of low-frequency signals, MEG data could be resampled to minimize the use of memory and to improve the computational efficiency. Resampling was done by decimating signals to extract the low frequency data. A low-pass anti-aliasing filter was applied before resampling. The high and low frequency pass-bands depended on the sampling rate and the frequency ranges of interest. In this study, two pass-bands of 1–70 Hz and 70–200 Hz were used. To compute the accumulated spectrogram, filtered MEG data were then segmented into small data segments. The length of the data segments depended on the time window of wavelet-transformation. In this study, we used a 5 s time-window and 600 frequency bins. There was no overlap between segments. Of note, the total length of recorded MEG data did not always match exactly with all of the segments. To solve this problem, data padding (typically, adding zero to make up enough data points for computing) was applied. If there were more than enough data points, the program also allowed for discarding of “extra” data points. The time duration of these segments depended on several factors including the available computer memory, storage spaces and research purposes. Once the time-frequency representations were computed, they were accumulated into one spectrum by adding them together. The “threshold” was used during data accumulating. There were two threshold values: a minimum threshold value and a maximum threshold value. If a time-frequency value was smaller than the minimum threshold value (e.g., background activity) or larger than the maximum threshold value (e.g., artifacts), the value was discarded. Accumulated spectrum is different from an averaged spectrum because the process of accumulating has several parameters: (1) accumulating has two thresholds; and (2) the accumulated data do not have to be averaged. Since the analysis of high-frequency components required high-sampling data, the re-sampling function was critical for low-frequency spectral analysis, which also minimized the use of computer memory. The workflows of data analyses at sensor levels are illustrated in Figure 1.

FIGURE 1

Figure 1. Workflow for computing accumulated spectrogram (left) and the basic principle of computing accumulated spectrogram (right). Since the analysis of high-frequency MEG signals requires high-sampling rate MEG data, MEG data are digitized in a high-frequency range. To improve the performance and optimize the use of computer memory for analyzing both low- and high-frequency MEG signals, the new method can re-sample MEG data dynamically according to the analysis frequency ranges. If the data points of the recorded data are smaller than the minimum data point of wavelet transform in frequency range, the “Data Padding” function can pad some data points so as to meet the requirements of wavelet transform. The “Thresholding” indicates that a spectral value can be rejected or accepted by the accumulated spectrogram according to a threshold value. MEG data recorded are waveforms, which are divided to segments (e.g., “Waveform 1,” “Waveform N”) to minimize the use of memory for wavelet transform (“Wavelet transform”). In the new method, wavelet transform transfers each segment of waveform data to a spectrum (e.g., “Spectrum 1,” “Spectrum N”). Of note, “N” indicates the total number of segments or spectra, which can be theoretically infinitely large. The “+” indicates the process of accumulation, which add all spectra together to produce an accumulated spectrum (“Accumulated Spectrum”). The left view of the sensor distribution of our MEG system is shown on the top right.

Morlet continuous wavelet transform was used for transforming time-domain data to frequency-domain data (see Figure 1). The Morlet wavelet was used because brain activity is nonstationary and the wavelet is better suited for nonstationary data (Ghuman et al., 2011). Wavelet transform can be described by the following equation:

$\begin{array}{cc}G(t,f)=\frac{1}{\sqrt{2\pi f}}{e}^{\left(\frac{-t^2}{2{\sigma }^2}\right)}{e}^{i2\pi ft}& (1)\end{array}$

In the above formula, t indicates time, f indicates frequency, and σ represents the standard deviation of the Gaussian curve in the time domain. To ensure stability of the wavelet transform, σ is typically larger than $\frac{5}{2\pi f}$. Since the wavelet convolution brings Gaussian temporal blurring with a standard deviation of σ, the effective number of independent samples is $\frac{N-1}{\sqrt{2\pi {(f_s\sigma )}^2}}$. The f_s represents the sampling frequency of the data and N represents the number of data points.

Since brain activation in a given time-window might occur in different frequency ranges and different frequencies might have different corresponding amplitudes, we used a different sigma value for each frequency to capture the time-frequency changes. Consequently, wavelet Equation (1) can be represented with an alternate representation in Equation (2) as follows:

$\begin{array}{cc}G(t,f)=C_{\sigma }{\pi }^{-\frac{1}{4}}{e}^{-\frac{1}{2}{t}^2}(e^{i\sigma t}-{\kappa }_{\sigma })& (2)\end{array}$

In the formula, t indicates time and f indicates frequency. Each wavelet transform has its own sigma value. Sigma is the scaling parameter that affects the width of the window. The sigma values are derived from the mother function in wavelet transform by computing the number of small waves for a time-frequency analyses (Ghuman et al., 2011). Sigma values could also be experimentally determined. κ_σ represents the admissibility and C_σ represents a normalized constant. σ represents the standard deviation of the Gaussian curve in the time domain. If signals appeared in the given sensitive time (a small sigma value) and sensitive frequency (a large sigma value) ranges, they would be enhanced.

An accumulated spectrum was defined as the time-frequency summation of a long-time or continuous recording which had a time period at least two times longer than that of the time window of the spectrum. The equation of computing accumulated spectra is given by:

$\begin{array}{cc}Atf(s,f)={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{f=1}^{F}G}}(t,f)& (3)\end{array}$

In Equation (3), Atf represents an accumulated spectrum; s indicates the time slice of the spectrum; f indicates frequency bands (or bins) of MEG data; T indicates total time points of MEG data and F indicate the total frequency bands. We defined s ≥ 1 and s ≤ T/2. From computer program point of view, the use of computer memory and storage space by Equation (3) depends on the s. Even though T could be infinitively increasing, the requirements for computer memory and storage remain the same. Consequently, the approach automatically avoided possible “overflow” or “out of space” problems in a long-time or continuous recording for capturing epileptic activity.

An accumulated spectrogram was computed by sequentially transforming each of the segments of waveform data to time-frequency representations using Morlet wavelet algorithm Equation (2) and then accumulating all the spectra together Equation (3). In this procedure, the different spectrograms of individual time segments were mathematically summed together to a single new overall spectrogram. An accumulated spectrogram can reveal brain activity in a consistent frequency range at multiple time windows. It can be considered as a “collective result” for a long-time recording. Figure 1 demonstrates the basic principles of computing an accumulated spectrogram. An accumulated spectrogram could reveal brain activity in a consistent frequency range while minimizing noise at random frequency ranges (Figure 1). Therefore, it could be considered to be a “collective result” of spatial- and frequency locked signals in multiple epochs of MEG data. To identify the frequency profile of the entire brain for a recording, we developed an accumulated global spectrogram. An accumulated global spectrogram was an averaged spectrogram of all accumulated spectrograms from the entire MEG sensor array. The accumulated global spectrogram was the “spatial summation” of the entire MEG sensor array' accumulated spectrograms. Since each sensor was positioned in a distinct location around the brain if there was a subject, an accumulated global spectrogram should represent the magnetic field of the entire brain. The mathematical principles have been described in previous reports (Rau et al., 2002). The neuromagnetic activity at each sensor was visualized with contour maps, which showed small spectrograms at the position of each MEG sensor. The equation of computing global spectrogram is given by:

$\begin{array}{cc}G(s,f)=\frac{1}{M}{\displaystyle \sum _{m=1}^{M}A}ft(s,f)& (4)\end{array}$

In Equation (4), G represents the global spectrogram; Atf represents an accumulated spectrum of one MEG sensor data; m indicates MEG sensor index and M indicates the total number of MEG sensors; s indicates the time slice of the spectrum; f indicates frequency bands (or bins) of MEG data. Since each sensor was positioned in a distinct location around the head (Figure 1), the global spectrogram is considered to be a “spatial summation” for each epoch of data (Xiang et al., 2009a).

Detection of Low- and High-Frequency MEG Signals at Source Levels

To detect low- and high-frequency neuromagnetic signals at source levels, two computing pipelines were developed. One computing pipeline generated multi-frequency datasets by processing MEG data with filter or wavelet transforms. MEG signals in multi-frequency datasets were in a set of frequency ranges. Of note, the frequency ranges depended on the research tasks and can be predefined. Another computing pipeline performed four tasks: (1) creating a three-dimensional source grid (3D grid), where each grid node represents a possible source; (2) conducting forward solution by calculating lead fields for each source (node) for the entire grid; (3) computing the lead field norm (or magnitude) and ranking the norm for each source for all sensors; (4) producing the node-beam lead field, performing single value decomposition (SVD) and calculating spatial filter weights. The node-beam lead field, which represents a form of sub-space solution, was completed by selecting a group of sensors which had a larger lead field norm. According to our tests, the optimal number of sensors for a node-beam lead field was in a range of 3 to M/3; here M indicates the total number of sensors of a whole cortex MEG system. For example, in our study, the total number of MEG sensor was 275. Thus, the suitable number of sensors that could be used for node-beam lead field was 3–91 (275/3). Of course, all sensors could be used for source scan. A small number of sensors was used in node-beam lead field because high-frequency brain signals were typically very weak and appeared only in a focal group of sensors. The node-beam sensors were also used to generate beam sensor MEG datasets so that the sensors in forward solution matched with measured magnetic signals. The final step was to compute source moments and to generate source data. Additional components were optional (red lines, which will be discussed in following sections). The main workflow for localizing both low- and high-frequency MEG signals is shown in Figure 2.

FIGURE 2

Figure 2. Workflow for computing accumulated source images. The workflow includes two main computing pipelines. One computing pipeline processes MEG data with filter or wavelet transforms so as to generate multi-frequency datasets. MEG signals in multi-frequency datasets are in a set of frequency ranges. Another computing pipeline works on several tasks, which included the creation of a three-dimensional source grid (3D grid), performing forward solution by calculating lead fields, ranking the norm for each source for all sensors, and performing SVD. The node-beam lead field is completed by selecting a group of sensors which have a larger lead field norm (or weights). Of note, each location in accumulated source imaging can have multiple parameters (e.g., “Frequency Index,” “Source Strength”). Some processes are optional (red lines) and additional parameters can also be added to the workflow.

Differing from the conventional volumetric source imaging or distributed source map, each grid node consisted of multiple data items including the strength and frequency of the source activity (Figure 2). Building on previous reports (Mosher and Leahy, 1998; Vrba and Robinson, 2001; De Gooijer-Van De Groep, 2013), the mathematic relationship between measured MEG data and source activity can be expressed as following equation:

$\begin{array}{cc}B=LQ+N& (5)\end{array}$

In Equation (4), B represents the MEG data; L represents the lead field, Q represents the source strength, and N represents the noise. For a given MEG dataset, B is known and L can be computed for each node with a forward solution. The forward solution in this study was computed according to Sarvas' formula for outside hemispherical conductors in Cartesian coordinates (Sarvas, 1987).

The determination of source strength and orientation of Q has been a challenge as discussed in many previous reports (Mosher et al., 1999; Huang et al., 2004; Robinson, 2004; De Munck and Bijma, 2009; Ou et al., 2009). According to our tests, the determination of MEG data in both low- and high-frequency ranges with conventional beamforming required considerable time and computing power to decompose MEG sensor data to subspaces because the data in both low- and high-frequency ranges had more data points as compared with the previous reports typically focusing on a single frequency range. However, for a given MEG data set in multiple frequency ranges in a limited time window (2 min in this study), the positions of sensor array and the 3D source grid were fixed; consequently, lead fields could be computed once and then used for both low and high-frequency ranges. Under these assumptions, we propose using SVD to decompose the lead field as following:

$\begin{array}{cc}L=US{V}^T& (6)\end{array}$

Where U ∈ R^mxm is an orthogonal (unitary in the complex case) matrix. The columns of U are the left singular vectors of L. V ∈ R^mxm is an orthogonal (unitary in the complex case) matrix. The columns of V are right singular vectors of L. S = diag(σ₁, σ₂,…σ_p) is an M × N diagonal matrix with p = (m, n) and σ₁, σ₂,…σ_p are the singular values of L. M indicates the number of sensors and N indicates the number of source orientations. For a single source, p = 3. The Moore-Penrose pseudo inverse of L is given by:

$\begin{array}{cc}{L}^{+}=V{S}^{+}{U}^T& (7)\end{array}$

Where S⁺ is a diagonal formed with the multiplicative inverses of the nonzero singular values of L placed on the diagonal. Assuming there was no noise (N = 0), the measured MEG data, B, can be described by the following equations:

$\begin{array}{cc}B=LQ=US{V}^{T}Q& (8)\end{array}$

$\begin{array}{cc}Q=B{L}^{-1}& (9)\end{array}$

By replacing L⁻¹ in Equation (9) with L in Equation (8), the estimated moment, $\overrightarrow{Q}$, can be computed with a SVD back substitution as described in the following equation:

$\begin{array}{cc}\overrightarrow{Q}=BV{S}^{+}{U}^T& (10)\end{array}$

Of note, L⁺, pseudo inverse of L, could be computed once and used for the analysis of data in all frequency ranges, which makes the computation of source strength and probability more efficient. In addition, once the $\overrightarrow{Q}$ is determined, virtual sensor spectrograms can be also computed with $\overrightarrow{Q}$ for each frequency range and time window.

$\begin{array}{cc}V(t,f)={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{f=1}^F|}}|\overrightarrow{Q}|{|}_2{(TF)}^{-1}& (11)\end{array}$

In Equation (11), V represents the computed virtual sensor spectral data. The t and T indicate time slice and total number of time windows, respectively. The f and F indicate frequency band and total number of frequency bands, respectively. Magnetic signals generated by $\overrightarrow{Q}$ can be computed with the follow equation:

$\begin{array}{cc}{X}_{cmp}=L\overrightarrow{Q}& (12)\end{array}$

where Xcmp represents computed magnetic signals at individual sensors from source $\overrightarrow{Q}$. We used Xmea to represent the measured magnetic signals at individual sensors, which were different from B in Equation (3), which represents MEG data in general.

Reliability Assessment of Source Activity at Low- and High-Frequency Ranges

To minimize the “ill-posed” inverse problem in MEG, the theory that a given MEG sensor data pattern may have an infinite number of possible “correct” answers (Hamalainen and Sarvas, 1987; Sarvas, 1987), we developed a channel-cross-channel (CxC) function to analyze the spatial pattern of MEG signals. Building on the use of covariance matrix for MEG beamforming in our previous studies (Kotecha et al., 2009; Gummadavelli et al., 2013), we applied a subtraction operation to all possible channel-pairs to generate a matrix which described the spatial gradient of magnetic signals among the sensors. Mathematically, each entry outside of the main diagonal in a CxC matrix represents the difference of a channel-pair. The diagonal entries represent the values of the corresponding sensors. To assess the reliability of source activity, the similarity of the measured MEG signal (Xmea) and the computed MEG signals (Xcmp) were statistically analyzed with the CxC matrix by computing the covariance and correlation factors with the following formulas:

$\begin{array}{cc}C\text{}\left(x_{mea},x_{cmp}\right)=\frac{{\displaystyle {\sum }_{i=1}^K(x_{mea}i-\overline{x_{mea}})(x_{cmp}i-\overline{x_{cmp}})}}{N-1}& (13)\end{array}$

$\begin{array}{cc}R\text{}\left(x_{mea},x_{cmp}\right)=\frac{C(x_{mea},x_{cmp})}{S{x}_{mea}S{x}_{cmp}}& (14)\end{array}$

Where C (x_mea, x_cmp) indicates the covariance and R (x_mea, x_cmp) indicates the correlation in the CxC matrices. The x_mea and x_cmp indicate signals in two channels which were paired for computing CxC. x_mea and x_cmp represent the mean of the signals in the measured and computed datasets, respectively. Sx_mea and Sx_cmp indicate the standard deviation of the signals in the two datasets, respectively. K indicates the number of sensors used for source estimation, which was smaller or equal to the total number of measuring sensors. To statistically determine the spatial correlations for each node in the 3D grid, t-values were computed for all sources.

$\begin{array}{cc}Tp=R\sqrt{\frac{K-2}{1-R^2}}& (15)\end{array}$

In Equation (15), Tp is the t-value of a source; R indicates the correlation of the measured and computed MEG signals for the source; K indicates the number of sensors related to the source.

A careful observation of Equation (13) could find that x_cmp is similar to the weights of the conventional beamforming because x_cmp represents signals from a predefined location and estimated source orientation. Similar to the conventional beamforming, the use of x_cmp could maximize signals from the source and minimize environmental noise and signals from other locations. For the analyses of multi-frequency signals, the location-orientation weights were computed from the optimal CxC matrix for each frequency. Thus, the source orientation was independent of frequency and only dependent on the orientation of the cortical normal vector. In other words, the solutions are approximations; the orientation portion was frequency independent.

Building on previous reports that the spectral signatures of low- and high-frequency signals at source levels can be measured with the combination of accumulated spectrogram and virtual sensors (Xiang et al., 2004, 2009a,b; Xiang and Xiao, 2009), the present study developed accumulated source imaging (Figure 2). With this technique, an accumulated source image was generated by accumulating all the source data computed for each location and each frequency band from the entire epoch of the MEG data. Of note, the computing of accumulated source images maintained spatial- and frequency-locked signals and minimized signals in random-space and frequency.

Magnetic Source Imaging with Multi-Parameters Per Location (MPPL)

This study moved one step further by developing magnetic source images with multi-values per location or MPPL. Specifically, each location has multi-parameters: (1) the first parameter describes the frequency range, which is represented with a frequency index for minimizing the use of computer memory and storage spaces; (2) the second parameter describes the strength of source activity; (3) the third parameter describes the reliability of the source; (4) the fourth parameter describes the Kurtosis or “peakedness” of source activity. The frequency index was directly obtained from the processed MEG data (the values of high-pass and low-pass filters or the frequency index in time-frequency representation). The strength of source activity was the source moment computed with Equation (10). The reliability could be computed with Equations (14) or (15). Building on previous report (Robinson et al., 2004), the Kurtosis was computed with following equation.

$\begin{array}{cc}K=\frac{{\displaystyle {\sum }_{t=1}^T{(q(t)-u)}^4}}{T{\sigma }_t^4}-3& (16)\end{array}$

Where T is the length of source data t in a time window, which has a mean of u and a standard deviation of σ. K represents the kurtosis values and is stored in parameter 4 in accumulated source imaging.

As shown in Figure 2, the analyses of MEG signals at both low- and high-frequency ranges generated more than one value for each location or each node of the 3D grid (e.g., strength, reliability and frequency of source activity). Notably, conventional magnetic source imaging, which encodes one value for one location or voxel, cannot represent the source data computed with the developed methods. The main differences between the new methods and existing methods are summarized in Table 1.

TABLE 1

Table 1. Differences between accumulated source imaging (ASI) and similar methods.

Source Localization with Accumulating

Accumulated source imaging was defined as the volumetric summation of source activity over a period of time which was at least two times longer than that of the time window of the source image. Of note, accumulated source imaging could have more than 1 time slices to reveal the fluctuation of source activity in space and time. Accumulated source imaging can be described as the following equation:

$\begin{array}{cc}Asi(r,s)={\displaystyle \sum _{t=1}^{t=n}Q}(r,t)& (17)\end{array}$

In Equation (17), Asi represents accumulated source strength at location r; s indicates the time slice; t indicates time point of MEG data; n indicates total time points of MEG data and Q indicate the source activity at source r and at time point t. We defined that s ≥ 1 and s ≤ n/2. From a computer program point of view, the use of computer memory and storage space by Equation (12) is dependent on the s for a fixed source imaging configuration (e.g., spatial resolution and dimension). Even though n could be infinitely increasing, the requirements for computer memory and storage remain the same. Consequently, the approach automatically avoided possible “overflow” or “out of space” problems in a long-time or continuous recording for capturing epileptic activity such as spikes. Since accumulated source imaging accumulates the results of source data, it is different from previous reports which compute a covariance matrix or kurtosis of sensor data for a long-time recording. Specifically, using a covariance matrix or kurtosis computed with sensor data for a long-time recording for source localized is based on the assumption that the source was stationary during the long-time recording. Our approach, on the other hand, did not make this assumption. Therefore, our approach has the capability to detect both stationary and nonstationary source activity.

MEG Experiments, MRI Scan and Intracranial Recordings

Participants

Ten healthy children (5 girls; 5 boys; age: 6–18 years; mean age: 12.8 years) were recruited for this study. Inclusion criteria were: (1) healthy without a history of neurological disorders or brain injuries; (2) age-appropriate functions including hearing, vision, and hand movement; (3) head movement during MEG recording was less than 5 mm. Ten pediatric patients (5 girls; 5 boys; age: 6–18 years; mean age: 12.7 years) with clinically diagnosed epilepsy were retrospectively studied. Patient inclusion criteria were: (1) clinically diagnosed epilepsy; (2) head movement during MEG recording was less than 5 mm; and (3) epileptic foci were confirmed with electrocorticography (ECoG) and/or neuroimaging data. Exclusion criteria were: (1) inability to remain still; and (2) presence of an implant such as a cochlear implant device, a pacemaker, or a neuro-stimulator containing electrical circuitry, generating magnetic signals, or having other metal that could produce visible magnetic noise (>6 pT) in the MEG data. Written consent, formally approved by the Institutional Review Board (IRB) at Cincinnati Children's Hospital Medical Center (CCHMC) and Nanjing Brain Hospital, was obtained from each healthy participant prior to testing. This study was approved by IRB at CCHMC.

MEG recordings

MEG signals were recorded in a magnetically shielded room (MSR) using a whole head CTF 275-Channel MEG system (VSM MedTech Systems Inc., Coquitlam, BC, Canada) in the MEG Center at CCHMC. Before data acquisition commenced, three electromagnetic coils were attached to the nasion, left and right pre-auricular points of each subject. These three coils were subsequently activated at different frequencies for measuring each subject's head position relative to the MEG sensors. Each subject lay comfortably in the supine position, his or her arms resting on either side, during the entire procedure. MEG data were recorded at a sampling rate of 4000 Hz. Continuous MEG recordings were completed for an epoch of 2 min. To ensure the reproducibility, at least two epochs were recorded for each subject. All MEG data were recorded with a noise cancellation of third order gradients and without on-line filtering. To identify system and environmental noise, we routinely recorded one MEG dataset without a subject immediately prior to the experiment.

MRI scan

Three-dimensional magnetic resonance imaging (MRI) was obtained using a 3-T Philips Achieva scanner (Philips Healthcare, Andover, MA). Three fiduciary marks were placed in identical locations to the positions of the 3 coils used in the MEG recordings with the aid of digital photographs to allow for an accurate co-registration of the 2 data sets. Subsequently, all anatomic landmarks were made identifiable in the MRIs.

Similar to previous reports (Xiang et al., 2009a), clinical intracranial electrocorticography (ECoG) data were retrospectively analyzed with the MEG results. Of the 10 patients, the 8 patients reported here had implantation of subdural electrodes and CCTV/EEG (VEEG) monitoring according to standard protocol at our hospital. Digital photos were taken before and during the operation to record the placements of the electrodes.

Implementation of the Algorithms

The aforementioned method for reconstruction of brain activity was implemented in MEG Processor with C/C++ on Windows platform (Xiang et al., 2010; Gummadavelli et al., 2013). MEG Processor was driven by its Windows interface. From the user perspective, its organization is contextual rather than linear: the multiple features from the software were not listed in long menus, they were accessible only when needed and were typically suggested within contextual popup menus or specific interface windows. This structure provided faster and easier access to requested functions.

Data Analyses

MEG data were visually inspected for artifacts. MEG waveforms with identifiable artifacts (amplitude >6 pT) were excluded from data analyses. Similar to previous reports (Xiang et al., 2009a), accumulated spectrograms, global spectrograms and spectral contour maps for all subjects were computed and analyzed. Before reconstructing brain activity for human MEG data, the head was modeled as a homogenous conducting sphere in order to account for volume-conducted return currents. The sphere model used in this study was a multiple local-sphere model, where each sphere (one per MEG sensor) was fit to a small patch of the head model (directly under the sensor) in order to better model the local return currents (Huang et al., 1999). The conducting boundary was defined with individual MRI, which was the inner skull. In other words, the best-fit sphere was fit to the scalp. From this head model, a whole-brain, subject-specific lead field was computed and used for magnetic source reconstruction. Accumulated source imaging and conventional beamforming (Vrba and Robinson, 2001) were implemented in MEG Processor for source estimation (Kotecha et al., 2009; Chen et al., 2010; Gummadavelli et al., 2013). CTF software package (VSM MedTech Systems Inc., Coquitlam, BC, Canada) was used to perform dipole fit analyses (Robinson et al., 2004; Kirsch et al., 2006). We used MNE (Gramfort et al., 2014) and Brainstorm (Tadel et al., 2011) to perform source estimation with Minimum-norm and multiple signal classification (MUSIC) algorithms, respectively.

To quantify the results, electrocorticography (ECoG) was used as the “gold standard” for defining epileptic zones. MEG sources were overlapped onto individual MRI data. Cerebral landmarks including the central sulcus, Sylvian fissure and the somatosensory cortex were used to define specific anatomical cortical brain regions (Agirre-Arrizubieta et al., 2009). The brain regions were the central, parietal, and occipital lobes. The frontal lobe was divided in the frontal superior, medial, inferior, and fronto-orbital regions; the temporal lobe into the lateral and mesial regions, the latter comprising the amygdala, the hippocampus, the parahippocampal gyrus, and the temporal-basal area. The inter-hemispheric region consisted of the mesial surface of the frontal, parietal, and occipital lobes (De Gooijer-Van De Groep, 2013). Similar to previous reports (Agirre-Arrizubieta et al., 2009; De Gooijer-Van De Groep, 2013), the concordance between MEG sources and ECoG was measured by determining if the interictal ECoG and MEG source locations were anatomically matched in the brain regions. We defined the sensitivity and specificity of the methods as followings.

$\begin{array}{cc}Sensitivity=\frac{TP}{TP+FN}& (18)\end{array}$

$\begin{array}{cc}Specificity=\frac{TN}{TN+FP}& (19)\end{array}$

Where TP represents the number of true positive (both MEG and ECoG showed epileptic foci); FN indicates the number of false negative (ECoG showed epileptic foci while MEG showed no epileptic focus); TN represents the number of true negative (both MEG method and ECoG showed no epileptic foci); FP represents number of false positive (MEG showed epileptic foci, but ECoG showed no epileptic focus).

Statistical Analysis

The comparisons of spectral and source data for epilepsy subjects and controls were performed with paired Student T-tests. The odds ratios of activity in brain areas in epilepsy subjects other than the areas identified in control groups for each frequency band were analyzed with Fisher's exact tests. Significance was accepted at the level of p < 0.05 for one comparison. Since multiple frequency bands and more than one source were analyzed, Bonferroni multiple comparison corrections were applied. Specifically, if multiple comparisons were to be taken into account then the significance level for any one of these comparisons was reduced from 0.05 to 0.05/parameter (e.g., for 9 frequency bands, p < 0.005).

Results

The size of 2 min MEG data digitized at a sampling rate of 4000 Hz (CTF MEG system, 275 sensors) was 0.597 GB. For time-frequency analyses, if the frequency bin of time-frequency transform was 600, the size of the time-frequency representation of 0.597 GB waveform data were 358 GB (600 × 0.597). When we computed the CxC data with time frequency data, the size of time-frequency based covariance matrices were 128164 GB (358 × 358 GB), which was approximately 125 TB. Of note, the source data computed from the time-frequency data would also be larger (>125 TB). Since the physical memory limit for windows 7 (64 bits, professional version) was 192 GB, the spectral data computed with the conventional time-frequency analysis method could not be stored in our Windows workstations because as it clearly exceeded the upper limits of the operating system. Alternatively, with accumulated spectrogram, we were able to limit the size of the spectrogram to 3 GB. Noticeably, the size and time required for computing an accumulated spectrogram mainly depended on the dimension of the accumulated spectrogram (number of frequency bins and time slices) and frequency ranges which could be adjusted by users. For an accumulated spectrogram with a dimension of 600 × 600 (600 frequency bins, 600 time slices) for a 2 min recording (sampling rate 4000 Hz), it took approximately 8.1 ± 0.03 h for data in 70–200 Hz, 1.3 ± 0.002 h for data in 1–70 Hz. Of note, the processing time would also depend on the speed of CPU and GPU, the number of programs running, the optimization of software compiling. In this study, we used two CPU (Intel Xeon, E4506, 2.13 Hz, each CPU has four cores). If GPU was used, the times were shortened to 42.5 ± 0.31 min for data in 70–200 Hz and 12.2 ± 0.009 min for 1–70 Hz, respectively. GPU could significantly shorten the computing time (p < 0.0001). Examples of accumulated spectrograms are shown in Figures 3–6. To identify high-frequency signals in multiple frequency bands with visual inspection, it took 2–3 days for a neurologist with 8 years of EEG/MEG experience.

FIGURE 3

Figure 3. Accumulated global spectrograms in three conditions. “Magnetic Noise” was computed with MEG data recorded without subjects. “Control Subject” was computed with MEG data recorded from a healthy child. “Epilepsy Subject” was computed with MEG data recorded from a child with epilepsy between seizures (interictal). The sampling rate of all MEG recordings was 6000 Hz. An accumulated global spectrogram represents the “spatial summation” of the entire MEG sensor array accumulated spectrograms. The three spectrograms show that the epilepsy subject has elevated spectral power as compared to the control subject.

FIGURE 4

Figure 4. Accumulated spectrograms show the well-known alpha activity in a healthy subject and an epilepsy subject. Noticeably, healthy subject (“Healthy Subject”) has a clear activity around 8–12 Hz (alpha activity). However, the epilepsy subject (“Epilepsy Subject”) has incrased activity in 2–4 Hz (low-frequency activity). The color bar shows the color coding of spectral power.

FIGURE 5

Figure 5. Global accumulated spectrograms from 10 epilepsy subjects and 1 healthy subject show the main frequency components of neuromagnetic signals in 70–200 Hz in epilepsy patients. Noticeably, the activity patterns vary across patients. The color bar shows the color coding of spectral power for all the global accumulated spectrograms.

FIGURE 6

Figure 6. Accumulated spectral contour maps from 10 epilepsy subjects and 1 healthy subject show the spatial distributions. Noticeably, the spectral distribution varies across patients. All the contour maps have the same orientation defined by the arrows: the “L” indicates the left side of the head and the “R” indicates the right side of the head. The “F” indicates that the upper part of the contour map represents the frontal region of the head; the “B” indicates that the lower part of the contour map represents the posterior region of the head. Each small circle represents one MEG sensor. The color bar shows the color coding of spectral power for all the contour maps.

The processing time for source scan with the conventional dynamic multi-dipole modeling (finding the 13 dipole for each time-slice) in multi-frequency ranges for recording at a sampling rate of 4000 Hz took 92.3 ± 0.4 h. However, our accumulated source imaging, which automatically scanned the entire brain for the same dataset took 12.7 ± 0.4 h. Of note, the approach was approximately 7.6 times faster than the conventional approach (p < 0.0001). If GPU was used, the time could be significantly shortened to approximately 6.3 ± 0.1 h. However, the use of GPU slowed down the user responses in our tests.

The global spectrograms of MEG datasets recorded from three conditions (no subject, healthy subjects and epilepsy subjects) showed that the epilepsy subjects had significantly increased spectral power. Figure 3 shows an example of global spectrograms in the three conditions. We noted that accumulated spectrograms revealed a clear alpha activity (approximately 8–12 Hz) in all healthy subjects (10/10, 100%) (Figure 4). Out of the 10 epilepsy patients, 9 patients showed increased spectral power in 70–200 Hz (9/10, 90%). Further analyses revealed that increased spectral power were around 106, 140, and 168 Hz in epilepsy patients (Figure 5). Figure 6 shows the spatial distributions of accumulated spectrograms in spectral contour maps.

Accumulated source imaging revealed focal increase of spectral power (Figure 7). Accumulated source imaging in low frequency ranges revealed that brain activities in 8–12 Hz (alpha) were localized to the occipital cortex in all the healthy subjects (10/10, 100%). However, alpha activity were localized to the occipital cortex in five epilepsy patients (5/10, 50%) and in nonoccipital cortices in other five patients (see Figure 8 for example). The five patients all had strong epileptic activity, which overshadowed and/or interrupted alpha activity. We noted that the increased spectral power at source levels varied among epilepsy patients.

FIGURE 7

Figure 7. An illustration of the basic principle of accumulated source imaging. The top waveforms show MEG data at sensor levels. The bottom images show individual structural magnetic resonance image and the region of interest (ROI, blue lines) for source scanning. MEG sensor data are firstly divided into small segments (e.g., “Sensor Data Segment 1,” “Sensor Data Segment 2,” “Sensor Data Segment N”). Volumetric sources are then produced by scanning the entire ROI with each segment of sensor data. The red, yellow and white small cubes indicate the sources (or voxels) identified. For illustration purposes, a very low resolution (12 millimeter) spatial resolution was used. An accumulated source image is generated by spatially adding all volumetric sources together. Of note, only sources reach certain thresholds (in this case, 75%) are added to accumulated source images, which differentiate this accumulating process from averaging. The color bar indicates the color coding of the source strength.

FIGURE 8

Figure 8. Accumulated source imaging shows low frequency brain activity in 8–12 Hz (alpha) in an epilepsy subject (“Epilepsy Subject”) and a healthy subject (“Control Subject”). Alpha activity is localized to the occipital cortex in the healthy subject. However, alpha activity is overshadowed by epileptic activity in the epilepsy subject. The epileptic activity is localized to the left and right parietal cortices in the epilepsy subject, which is concordant with clinical findings.

Accumulated source imaging showed that 9 out of the 10 epilepsy patients (9/10, 90%) had increased focal spectral power in high-frequency ranges at source levels. The epileptic areas localized by accumulated source imaging were concordant with clinical data. Figure 9 shows an example of epileptic foci volumetrically localized with high-frequency accumulated source imaging (70–200 Hz). The sensitivity and specificity of all subjects are shown in Table 2.

FIGURE 9

Figure 9. A digital photo of intracranial recording (“ECoG”) and accumulated source imaging (“ASI”) show the concordance of the two methods. The two images are placed in the similar orientation. “Frontal” indicate the frontal cortex; “Temporal” indicates temporal lobe. The left green arrow points to the epileptic area invasively defined by intracranial recording; the right blue arrow points to the epileptic area noninvasively localized with high-frequency neuromagnetic signals. Noticeably, the areas are matched in gyrus level. The color bar shows the color coding of accumulated source imaging. The value of the source voxel is normalized T value (no unit).

TABLE 2

Table 2. The sensitivity and specificity of five MEG source localization methods.

Discussion

The present study demonstrated an approach for detecting both low- and high-frequency neuromagnetic signals by integrating time-frequency transform, source localization, accumulation and MPPL algorithms into a comprehensive and systematic processing package. The strengths of our methodologies are reflected by the major features of our signal processing algorithms as well as their abilities to resolve the difficulties associated with the large data volume, multi-modality data and its clinical applicability.

Features of Our Signal Processing Algorithms

One of the unique features in our wavelet transform algorithm was that the sigma value (number of waves) could be dynamically changed so as to match the neurophysiological patterns. For example, neuromagnetic signals from a brain area may appear in multiple frequency ranges but in a similar time window. The conventional wavelet algorithm typically gives a wide time-window for low signals and a narrow time-window for high-frequency signals, which is not well-suited for analysis of brain activity. The improved wavelet transform algorithm in the present study could solve this problem by dynamically changing the sigma values so as to adjust the time-window for a better analysis of brain activity.

Accumulating algorithms in the computing of accumulated spectrograms provides a novel method for handling the large datasets obtained when analyzing both low and high-frequency MEG signals. Integration of time-frequency analysis and accumulation into a workflow system is a novel neuroimaging data processing algorithm technique that can summarize and visualize high-frequency signals with a few images.

CxC matrices and functions are critical to the study of neural HFOs. Since high-frequency signals are typically obscured by low-frequency signals (Xiang et al., 2009a), time-frequency representations were normalized according to the magnitude of each frequency bin across all MEG sensors to ensure that all frequency bins contributed equally to the source reconstruction. The time-frequency matrix allows for matrix operations such as subtraction of control state MEG signals from the activation state MEG signals, whose purpose is to increase signal-to-noise ratio (SNR) or to maximize the signal power at a peak frequency (or a frequency of interest) while simultaneously minimizing it at the neighboring surrounding frequency bins. CxC matrices based on the time-frequency data provide unique spatial patterns and gradients of magnetic fields for determining high-frequency sources.

The major differences between our technique and existing methods of volumetric imaging such as beamforming, minimum-norm are the features of accumulation and MPPL, which are more than a source localization algorithm. To our knowledge, none of the existing methods have the features of accumulation and MPPL. It is necessary to point it out that, some existing methods have internally fixed frequency ranges (e.g., 20–70 Hz) (Robinson et al., 2004), which could not be directly compared in our tests because our method was designed to analyze both low- and high-frequency signals (multi-frequencies). Of note, each method has its strengths and weaknesses (De Gooijer-Van De Groep, 2013). According to on our clinical experience, the unique features of the method are clinically important and necessary. For example, the conventional beamforming could also been used to detect multi-frequency signals (Vrba and Robinson, 2001). However, the conventional beamforming is based on covariance matrices, which are computed from data in long time-windows (if the time-windows are short, the sizes of the source data would be a problem). The process assumes that the brain activity is stationary in the time-window (Vrba and Robinson, 2001), which may be not true for real epileptic activity (Zijlmans et al., 2012b). By using accumulation algorithms, our approach does not making any assumption about the stationarity of the sources. Another example is SAM (g2). SAM (g2) is an outstanding method for detecting excess kurtosis (Kirsch et al., 2006). SAM (g2) is designed for detecting rare events (spikiness activity). It has been shown that combining SAM(g2) and other methods such as MUSIC gives the best clinical results (De Gooijer-Van De Groep, 2013). The development of MPPL in our method enables us to implement both kurtosis and other algorithms by using multiple parameters during source analyses. Consequently, both rare events (kurtosis) and common events (frequent spikes) could be detected by our methods.

Several Major Challenges Resolved with the Current Methodologies

Data volume challenge

Our results are consistent with previous reports (Blanco et al., 2011), the size of high sampling rate MEG/EEG data could be in the magnitude of TB (>12 TB). This was particularly true for multi-frequency spectral data (>125 TB). However, by using accumulating techniques, we were able to minimize the size of MEG data to less than 10 GB without losing the high-frequency information. Although accumulated spectrogram was utilized with MEG in the present study, the same technology can also be used in the analysis of EEG and intracranial EEG. In current clinical research, the detection and labeling of interictal and ictal epileptiform activity in intracranial EEG recordings is performed by expert review. This manual method has been known to be associated with a poor inter-reviewer reliability (Benbadis et al., 2009). In addition, manual review is not feasible for large data sets because it is very time consuming and labor-intensive (Restuccia et al., 2011; Andrade-Valenca et al., 2012; Dumpelmann et al., 2012; Jacobs et al., 2012; Zijlmans et al., 2012a; Haegelen et al., 2013; Stacey et al., 2013). Alternatively, accumulated source imaging can automatically analyze large data sets and provide images for experts to review. This new method can be used in combination with experts' review to verify its sensitivity and specificity and to advance our understanding of the relationship between HFO and epilepsy. According to our data, the new method can be further developed as a fully automatic detector with high specificity and sensitivity.

The development of methods for analysis of a substantial amount of MEG/EEG data has become an important research area. For example, data mining has been developed for the analysis of HFOs in epilepsy patients (Blanco et al., 2011; Worrell et al., 2012). Blanco et al. (2011) reported a quantitative analysis of HFOs and their rates of occurrence in 9 patients with neocortical epilepsy and two control patients with no history of seizures (sampling rate: 32,556 Hz). Using the data mining approach, they found that a cluster of ripple frequency oscillations with a median spectral centroid of 137 Hz is increased in the seizure-onset zone more frequently than a cluster of fast ripple frequency oscillations (median spectral centroid = 305 Hz). Our results are consistent with their findings. The relative rate of ripple frequency oscillations is an interesting potential biomarker for the epileptic neocortex, but larger prospective studies correlating HFOs rates with seizure-onset zones, resected tissue and surgical outcomes are required to determine the true predictive value of this line of research (Montazeri et al., 2009; Blanco et al., 2011; Worrell et al., 2012). However, to our understanding, algorithmic requirements differ substantially for data mining and for topological (feature) data analysis. In particular, little is known about the locations of high frequency brain signals and their relationship to neurological disorders. In this regard, one of the unique features of accumulated source imaging is its ability to localize and visualize epileptic activity in both low- and high-frequency ranges for correlating locations of brain signals to neurological disorders.

Multi-modality imaging data challenge

Our data have also shown that functional MEG data can be seamlessly integrated into structural MRI data. In comparison to conventional source imaging, one important feature of accumulated source imaging is MPPL. MPPL analysis results in multi-values per voxel in 3D images. One parameter is dedicated to the frequency signature, which is important for visualizing HFOs. For example, HFOs may be a band-limited event (Crepon et al., 2010) or they can be a broadband event (Staba and Bragin, 2011; Worrell et al., 2012; Zijlmans et al., 2012b). By visualizing the frequency in the imaging data, we can better address many current questions in the study of HFOs (Engel et al., 2009; Staba and Bragin, 2011; Worrell et al., 2012; Zijlmans et al., 2012b). For example, if HFOs are band-limited, should there be specific spectral boundaries? In other words, should a HFO be defined as an isolated event in the time–frequency map, or could it contain a variety of frequencies within a range? Since spontaneous activity can occur in multi-frequency ranges, we consider the development of accumulated source imaging with MPPL to be important for multi-modality integration because the unique parameters from MEG is encoded in each location (voxel) and can be easily integrated into other modalities without losing any information.

Integration of accumulating and source localization into a systematic approach is a powerful neuroimaging data processing technique that could simplify multi-modality analyses. For example, epileptic foci defined by HFOs are not time-locked and can spontaneously occur at any time point or time window. Without accumulation, thousands of MEG source images may need to be integrated into structural MRIs, which would be time-consuming. With accumulated source imaging, epileptic foci can be captured and summarized as a few images, which could be easily integrated into structural MRIs. As demonstrated in the Results section, accumulated source imaging is a potential powerful technique for multi-modality analysis of epileptic foci. Importantly, our software packages and libraries are based on C/C++. This methodology can be similarly implemented in more advanced computer systems such as cluster/GPU/FPGA or cloud/HPC. By using those advanced computer technologies, accumulated source imaging can be computed in a timely manner and routinely used in clinical practice in the future.

Clinical applicability challenge

Building on previous reports (Robinson et al., 2004; Kirsch et al., 2006) and our clinical observation, we developed the aforementioned method for detecting both low- and high-frequency brain signals. The proposed framework and architecture may also solve a few problems occurring in our clinical practice. First, this method provides an objective means of data analysis. The existing and currently practiced method for identifying epileptic spikes relies on visual inspection, which is subjective. Second, the proposed method provides meaningful quantitative source data, which are not available in conventional visual identification of epileptic spikes. Third, the new method can provide novel frequency descriptions about aberrant brain activity. In addition, the newly developed method semi-automatically quantifies MEG spectral power and source activity. Moreover, the new method has the capability of detecting and localizing high-frequency epileptic signals, a feat impossible to achieve with the conventional visual inspection of waveforms.

The results of spectral data showed that accumulated source imaging may play a key role in the differentiation of true HFOs from environmental noise in pre-operative workup for epilepsy surgery. It is well known that low-frequency signals may generate high-frequency harmonics. Since any harmonic of a high-frequency signal will localize to the corresponding low-frequency component, accumulated source imaging (which encodes both frequency and spatial information) can automatically reveal the main frequency by comparing the spectral power in the location of question. If HFOs were localized to a brain area which did not have low-frequency signals, the location would be an index for true HFOs. Since digital filtering may be used in the analysis of HFOs, the filter characteristics must be taken into account to avoid the detection of false oscillations (Benar et al., 2010). It has been noted that sharp transients with spectral content in HFO bands but without actual HFO in the raw data may be generated by filtering. According to our observation, such false oscillations are typically the result of the additive superposition of harmonics. They do not have a consistent spatial pattern in CxC and cannot be consistently localized to a location in the brain.

Accumulated source imaging may also play a key role in the differentiation of brain HFOs from artifacts in clinical practice. Raw MEG data contain a mixture of high-frequency brain signals and a variety of artifacts and noise. A major obstacle to HFO research is the unfortunate fact that various muscle activities typically result in prominent increases in gamma power (>25 Hz), and contaminate the recorded signal in the HFO spectrum. Myogenic activity interferes with the detection of HFO and represents a significant and often under-estimated challenge in clinical and basic research. For many years, intracranial EEG recordings were assumed to be largely, if not completely, immune to eye movements and muscle artifacts. This assumption has recently been proven to be erroneous (Ball et al., 2009; Jerbi et al., 2009; Kovach et al., 2011). To solve these problems, we tried a different approach. Since high-frequency signals are typically obscured by low-frequency signals, time-frequency representations could be normalized according to the magnitude of each frequency bin across all MEG sensors to ensure that all frequency bins contributed equally to the source reconstruction. The time-frequency matrix allows for matrix operations to maximize the signal power at a peak frequency while simultaneously minimizing it at the neighboring surrounding frequency bins. CxC matrices based on the time-frequency data provide spatial patterns and gradients of magnetic fields for accurately determining epileptic foci for clinical purposes. In addition, the method was able to show multiple metrics of source analyses. It is important to be able to visualize different metrics of the source data because the frequency, strength, reliability/probability and kurtosis are important for us to correctly interpret the results. According to our pilot data, very strong high-frequency sources (>100 Hz) typically pinpointed to the epileptogenic zones and the removal of these zones would likely result in good surgical outcomes and ultimately seizure freedom (Xiang et al., 2009a). Thus, we postulate that these multiple metrics of source data will allow discrimination among pathological, benign or artifactural source signals in the future.

Although our newly developed method showed promising results for detecting both low- and high-frequency brain signals, several weaknesses and problems have been identified and need to be addressed in the future. Specifically, we used multiple local spheres in the computation of forward solution, which did not address the effect of the inferior conductive boundary of the skull that is not proximal to any MEG sensors. This leads to questions as to the accuracy of the forward model for “deep” sources, as may be encountered in temporal lobe epilepsy. A model based on the superior and lateral curvature of the head may mitigate this problem. The Sarvas forward solution, as applied to each sensor's sphere origin could only compute the field due to the tangential components of the dipole moment. Given that there were multiple sphere origins that might be in the vicinity of one another, the dipole orientation and therefore the weights might be “confused” by the rapid change (with location) of the tangential orientation. The number of sensors in node-beam lead field was experimentally determined by using from 3 to 275 sensors. Since three sensors have only two degrees of freedom left to attenuate unwanted interference or brain signal, at least five sensors are necessary to discriminate between brain source and interference (with 3rd gradient compensation). Of note, we continue to perform research into improving our methodology to overcome some of these limitations.

We also noted that source activities in a few subjects were close to the brain-stem. Those activities might be an artifact or localization problem or real sources. According to our data, the activity in the center of the brain is more than likely real for several reasons: (1) MEG data recorded without subjects did not show similar sources. Thus, it is unlikely that the sources are from system artifacts (e.g., hardware, software or localization algorithms); (2) the shape of the volumetric sources appears to mimic the structure of individual structural MRIs in subjects. If the deep sources are system artifacts or localization problems, they should not mimic the structure of individual MRI. (3) There are reports showing that MEG can detect and localize source in the deep brain areas. (4) We were very careful to exclude artifact by excluding subjects with magnetic artifacts, recording MEG data with third-gradient noise cancellation and by visually inspecting the MEG data. However, MEG was not sensitive to sources in the brain areas and the source images of the deep source are more diffuse as compared to the surface sources. Therefore, further investigation and verification are necessary. We consider that closely following the guidance recommended by Gross and colleagues may further improve the quality of the data (Gross et al., 2013). In particular, by using multiple metrics of source analysis, we can incorporate new algorithms into data analysis by adding one source parameter. For example, Fatima and colleagues have developed a novel method to significantly improve the detection and localization of MEG sources by using independent component analysis (ICA) (Fatima et al., 2013), which can be incorporated into our source analysis pipeline to correct artifact and improve source localization. For resampling, one could also reduce the number of samples by high pass filtering the raw data and heterodyning it down to baseband, followed by decimation. The software and supplementary materials, which implemented the aforementioned algorithms, are freely available from the following website (http://sdrv.ms/PHenGK) for other researchers to test, reproduce, and improve the methods.

Summary

In summary, the present study has demonstrated that accumulated source imaging is a new powerful technique for quantitatively and objectively analyzing MEG signals at source levels. By volumetrically scanning sources and accumulating source information, accumulated source imaging could handle very large datasets and extract meaningful spatial information about brain activity. Accumulated source imaging based on HFO detection may play a key role in differentiating brain activity from environmental noise and muscle artifacts. Though further verification is necessary, we believe that the next study should focus on using more advanced computer systems such as cluster/GPU/FPGA/cloud/HPC to significantly improve the performance of the proposed methods for clinical applications in the future.

Author Contributions

Jing Xiang: study design and methodological development. Hisako Fujiwara, Nat Hemasilpin, Douglas F. Rose: data acquisition. Abraham Korman, Fawen Zhang, Jing Xiang: data analysis and manuscript preparation. Qian Luo, Rupesh Kotecha, Abraham Korman, Huan Luo: manuscript preparation.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This study was partially supported by a Trustee Grant from Cincinnati Children's Hospital Medical Center and Grant Number R21NS072817 from the National Institute of Neurological Disorders and Stroke (NINDS), the National Institutes of Health. We thank Drs. Lu Tang, Hongxing Liu, Xiaoshan Wang and Jingping Shi for helping with data analyses and sharing MEG and clinical data. We are very grateful for the two anonymous reviewers for very important comments and suggestions for improving the manuscript and improving our methods.

References

Agirre-Arrizubieta, Z., Huiskamp, G. J., Ferrier, C. H., Van Huffelen, A. C., and Leijten, F. S. (2009). Interictal magnetoencephalography and the irritative zone in the electrocorticogram. Brain 132, 3060–3071. doi: 10.1093/brain/awp137

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Andrade-Valenca, L., Mari, F., Jacobs, J., Zijlmans, M., Olivier, A., Gotman, J., et al. (2012). Interictal high frequency oscillations (HFOs) in patients with focal epilepsy and normal MRI. Clin. Neurophysiol. 123, 100–105. doi: 10.1016/j.clinph.2011.06.004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ball, T., Kern, M., Mutschler, I., Aertsen, A., and Schulze-Bonhage, A. (2009). Signal quality of simultaneously recorded invasive and non-invasive EEG. Neuroimage 46, 708–716. doi: 10.1016/j.neuroimage.2009.02.028

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Benar, C. G., Chauviere, L., Bartolomei, F., and Wendling, F. (2010). Pitfalls of high-pass filtering for detecting epileptic oscillations: a technical note on “false” ripples. Clin. Neurophysiol. 121, 301–310. doi: 10.1016/j.clinph.2009.10.019

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Benbadis, S. R., Lafrance, W. C. Jr., Papandonatos, G. D., Korabathina, K., Lin, K., Kraemer, H. C., et al. (2009). Interrater reliability of EEG-video monitoring. Neurology 73, 843–846. doi: 10.1212/WNL.0b013e3181b78425

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Blanco, J. A., Stead, M., Krieger, A., Stacey, W., Maus, D., Marsh, E., et al. (2011). Data mining neocortical high-frequency oscillations in epilepsy and controls. Brain 134, 2948–2959. doi: 10.1093/brain/awr212

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Brinkmann, B. H., Bower, M. R., Stengel, K. A., Worrell, G. A., and Stead, M. (2009). Large-scale electrophysiology: acquisition, compression, encryption, and storage of big data. J. Neurosci. Methods 180, 185–192. doi: 10.1016/j.jneumeth.2009.03.022

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Chen, Y., Xiang, J., Kirtman, E. G., Wang, Y., Kotecha, R., and Liu, Y. (2010). Neuromagnetic biomarkers of visuocortical development in healthy children. Clin. Neurophysiol. 121, 1555–1562. doi: 10.1016/j.clinph.2010.03.029

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Crepon, B., Navarro, V., Hasboun, D., Clemenceau, S., Martinerie, J., Baulac, M., et al. (2010). Mapping interictal oscillations greater than 200 Hz recorded with intracranial macroelectrodes in human epilepsy. Brain 133, 33–45. doi: 10.1093/brain/awp277

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Dalal, S. S., Guggisberg, A. G., Edwards, E., Sekihara, K., Findlay, A. M., Canolty, R. T., et al. (2008). Five-dimensional neuroimaging: localization of the time-frequency dynamics of cortical activity. Neuroimage 40, 1686–1700. doi: 10.1016/j.neuroimage.2008.01.023

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

De Gooijer-Van De Groep, K. L., Leijten, F. S., Ferrier, C. H., and Huiskamp, G. J. (2013). Inverse modeling in magnetic source imaging: comparison of MUSIC, SAM(g2), and sLORETA to interictal intracranial EEG. Hum. Brain Mapp. 34, 2032–2044. doi: 10.1002/hbm.22049

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

De Munck, J. C., and Bijma, F. (2009). Three-way matrix analysis, the MUSIC algorithm and the coupled dipole model. J. Neurosci. Methods 183, 63–71. doi: 10.1016/j.jneumeth.2009.06.040

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Dumpelmann, M., Jacobs, J., Kerber, K., and Schulze-Bonhage, A. (2012). Automatic 80–250 Hz “ripple” high frequency oscillation detection in invasive subdural grid and strip recordings in epilepsy by a radial basis function neural network. Clin. Neurophysiol. 123, 1721–1731. doi: 10.1016/j.clinph.2012.02.072

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Engel, J. Jr., Bragin, A., Staba, R., and Mody, I. (2009). High-frequency oscillations: what is normal and what is not? Epilepsia 50, 598–604. doi: 10.1111/j.1528-1167.2008.01917.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fatima, Z., Quraan, M. A., Kovacevic, N., and McIntosh, A. R. (2013). ICA-based artifact correction improves spatial localization of adaptive spatial filters in MEG. Neuroimage 78, 284–294. doi: 10.1016/j.neuroimage.2013.04.033

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ghuman, A. S., Mcdaniel, J. R., and Martin, A. (2011). A wavelet-based method for measuring the oscillatory dynamics of resting-state functional connectivity in MEG. Neuroimage 56, 69–77. doi: 10.1016/j.neuroimage.2011.01.046

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gotman, J. (2010). High frequency oscillations: the new EEG frontier? Epilepsia 51(Suppl. 1), 63–65. doi: 10.1111/j.1528-1167.2009.02449.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gramfort, A., Luessi, M., Larson, E., Engemann, D., Strohmeier, D.A., Brodbeck, C., et al. (2014). MNE software for processing MEG and EEG data. Neuroimage 86, 446–460. doi: 10.1016/j.neuroimage.2013.10.027

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gross, J., Baillet, S., Barnes, G. R., Henson, R. N., Hillebrand, A., Jensen, O., et al. (2013). Good practice for conducting and reporting MEG research. Neuroimage 65, 349–363. doi: 10.1016/j.neuroimage.2012.10.001

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Guggisberg, A. G., Dalal, S. S., Findlay, A. M., and Nagarajan, S. S. (2007). High-frequency oscillations in distributed neural networks reveal the dynamics of human decision making. Front. Hum. Neurosci 1:14. doi: 10.3389/neuro.09.014.2007

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gummadavelli, A., Wang, Y., Guo, X., Pardos, M., Chu, H., Liu, Y., et al. (2013). Spatiotemporal and frequency signatures of word recognition in the developing brain: a magnetoencephalographic study. Brain Res. 1498, 20–32. doi: 10.1016/j.brainres.2013.01.001

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Haegelen, C., Perucca, P., Chatillon, C. E., Andrade-Valenca, L., Zelmann, R., Jacobs, J., et al. (2013). High-frequency oscillations, extent of surgical resection, and surgical outcome in drug-resistant focal epilepsy. Epilepsia 54, 848–857. doi: 10.1111/epi.12075

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hamalainen, M. S., and Sarvas, J. (1987). Feasibility of the homogeneous head model in the interpretation of neuromagnetic fields. Phys. Med. Biol. 32, 91–97. doi: 10.1088/0031-9155/32/1/014

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Huang, M. X., Mosher, J. C., and Leahy, R. M. (1999). A sensor-weighted overlapping-sphere head model and exhaustive head model comparison for MEG. Phys. Med. Biol. 44, 423–440. doi: 10.1088/0031-9155/44/2/010

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Huang, M. X., Shih, J. J., Lee, R. R., Harrington, D. L., Thoma, R. J., Weisend, M. P., et al. (2004). Commonalities and differences among vectorized beamformers in electromagnetic source imaging. Brain Topogr. 16, 139–158. doi: 10.1023/B:BRAT.0000019183.92439.51

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jacobs, J., Staba, R., Asano, E., Otsubo, H., Wu, J. Y., Zijlmans, M., et al. (2012). High-frequency oscillations (HFOs) in clinical epilepsy. Prog. Neurobiol. 98, 302–315. doi: 10.1016/j.pneurobio.2012.03.001

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jacobs, J., Zelmann, R., Jirsch, J., Chander, R., Dubeau, C. E., and Gotman, J. (2009). High frequency oscillations (80–500 Hz) in the preictal period in patients with focal seizures. Epilepsia 50, 1780–1792. doi: 10.1111/j.1528-1167.2009.02067.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jerbi, K., Freyermuth, S., Dalal, S., Kahane, P., Bertrand, O., Berthoz, A., and Lachaux, J. P. (2009). Saccade related gamma-band activity in intracerebral EEG: dissociating neural from ocular muscle activity. Brain Topogr. 22, 18–23. doi: 10.1007/s10548-009-0078-5

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jirsch, J. D., Urrestarazu, E., Levan, P., Olivier, A., Dubeau, F., and Gotman, J. (2006). High-frequency oscillations during human focal seizures. Brain 129, 1593–1608. doi: 10.1093/brain/awl085

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kirsch, H. E., Robinson, S. E., Mantle, M., and Nagarajan, S. (2006). Automated localization of magnetoencephalographic interictal spikes by adaptive spatial filtering. Clin. Neurophysiol. 117, 2264–2271. doi: 10.1016/j.clinph.2006.06.708

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kotecha, R., Xiang, J., Wang, Y., Huo, X., Hemasilpin, N., Fujiwara, H., et al. (2009). Time, frequency and volumetric differences of high-frequency neuromagnetic oscillation between left and right somatosensory cortices. Int. J. Psychophysiol. 72, 102–110. doi: 10.1016/j.ijpsycho.2008.10.009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kovach, C. K., Tsuchiya, N., Kawasaki, H., Oya, H., Howard, M. A. 3rd., and Adolphs, R. (2011). Manifestation of ocular-muscle EMG contamination in human intracranial recordings. Neuroimage 54, 213–233. doi: 10.1016/j.neuroimage.2010.08.002

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Le Van Quyen, M., Staba, R., Bragin, A., Dickson, C., Valderrama, M., Fried, I., et al. (2010). Large-scale microelectrode recordings of high-frequency gamma oscillations in human cortex during sleep. J. Neurosci. 30, 7770–7782. doi: 10.1523/JNEUROSCI.5049-09.2010

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Levesque, M., Bortel, A., Gotman, J., and Avoli, M. (2011). High-frequency (80–500 Hz) oscillations and epileptogenesis in temporal lobe epilepsy. Neurobiol. Dis. 42, 231–241. doi: 10.1016/j.nbd.2011.01.007

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Matsumoto, A., Brinkmann, B. H., Matthew Stead, S., Matsumoto, J., Kucewicz, M., Marsh, W. R., et al. (2013). Pathological and physiological high-frequency oscillations in focal human epilepsy. J. Neurophysiol. 110, 1958–1964. doi: 10.1152/jn.00341.2013

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Montazeri, N., Shamsollahi, M. B., and Hajipour, S. (2009). MEG based classification of wrist movement. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2009, 986–989. doi: 10.1109/IEMBS.2009.5334472

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mosher, J. C., and Leahy, R. M. (1998). Recursive MUSIC: a framework for EEG and MEG source localization. IEEE Trans. Biomed. Eng. 45, 1342–1354. doi: 10.1109/10.725331

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mosher, J. C., Leahy, R. M., and Lewis, P. S. (1999). EEG and MEG: forward solutions for inverse methods. IEEE Trans. Biomed. Eng. 46, 245–259. doi: 10.1109/10.748978

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ou, W., Hamalainen, M. S., and Golland, P. (2009). A distributed spatio-temporal EEG/MEG inverse solver. Neuroimage 44, 932–946. doi: 10.1016/j.neuroimage.2008.05.063

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pail, M., Halamek, J., Daniel, P., Kuba, R., Tyrlikova, I., Chrastina, J., et al. (2013). Intracerebrally recorded high frequency oscillations: simple visual assessment versus automated detection. Clin. Neurophysiol. 124, 1935–1942. doi: 10.1016/j.clinph.2013.03.032

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Papadelis, C., Poghosyan, V., Fenwick, P. B., and Ioannides, A. A. (2009). MEG's ability to localise accurately weak transient neural sources. Clin. Neurophysiol. 120, 1958–1970. doi: 10.1016/j.clinph.2009.08.018

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pulvermuller, F., Birbaumer, N., Lutzenberger, W., and Mohr, B. (1997). High-frequency brain activity: its possible role in attention, perception and language processing. Prog. Neurobiol. 52, 427–445. doi: 10.1016/S0301-0082(97)00023-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rau, R., Raschka, C., and Koch, H. J. (2002). Uniform decrease of alpha-global field power induced by intermittent photic stimulation of healthy subjects. Braz. J. Med. Biol. Res. 35, 605–611. doi: 10.1590/S0100-879X2002000500014

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Restuccia, D., Del Piero, I., Martucci, L., and Zanini, S. (2011). High-frequency oscillations after median-nerve stimulation do not undergo habituation: a new insight on their functional meaning? Clin. Neurophysiol. 122, 148–152. doi: 10.1016/j.clinph.2010.06.008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Robinson, S. E. (2004). Localization of event-related activity by SAM(erf). Neurol. Clin. Neurophysiol. 2004, 109.

Pubmed Abstract | Pubmed Full Text

Robinson, S. E., Nagarajan, S. S., Mantle, M., Gibbons, V., and Kirsch, H. (2004). Localization of interictal spikes using SAM(g2) and dipole fit. Neurol. Clin. Neurophysiol. 2004, 74.

Pubmed Abstract | Pubmed Full Text

Sarvas, J. (1987). Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem. Phys. Med. Biol. 32, 11–22. doi: 10.1088/0031-9155/32/1/004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Srejic, L. R., Valiante, T. A., Aarts, M. M., and Hutchison, W. D. (2013). High-frequency cortical activity associated with postischemic epileptiform discharges in an in vivo rat focal stroke model. J. Neurosurg. 118, 1098–1106. doi: 10.3171/2013.1.JNS121059

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Staba, R. J., and Bragin, A. (2011). High-frequency oscillations and other electrophysiological biomarkers of epilepsy: underlying mechanisms. Biomark. Med. 5, 545–556. doi: 10.2217/bmm.11.72

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stacey, W. C., Kellis, S., Greger, B., Butson, C. R., Patel, P. R., Assaf, T., et al. (2013). Potential for unreliable interpretation of EEG recorded with microelectrodes. Epilepsia 54, 1391–1401. doi: 10.1111/epi.12202

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tadel, F., Baillet, S., Mosher, J. C., Pantazis, D., and Leahy, R. M. (2011). Brainstorm: a user-friendly application for MEG/EEG analysis. Comput. Intell. Neurosci. 2011:879716. doi: 10.1155/2011/879716

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tort, A. B., Scheffer-Teixeira, R., Souza, B. C., Draguhn, A., and Brankack, J. (2013). Theta-associated high-frequency oscillations (110–160 Hz) in the hippocampus and neocortex. Prog. Neurobiol. 100, 1–14. doi: 10.1016/j.pneurobio.2012.09.002

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Uhlhaas, P. J., Pipa, G., Neuenschwander, S., Wibral, M., and Singer, W. (2011). A new look at gamma? High- (>60 Hz) gamma-band activity in cortical networks: function, mechanisms and impairment. Prog. Biophys. Mol. Biol. 105, 14–28. doi: 10.1016/j.pbiomolbio.2010.10.004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Uhlhaas, P. J., and Singer, W. (2013). High-frequency oscillations and the neurobiology of schizophrenia. Dialogues Clin. Neurosci. 15, 301–313.

Pubmed Abstract | Pubmed Full Text

Van Essen, D. C., Ugurbil, K., Auerbach, E., Barch, D., Behrens, T. E., Bucholz, R., et al. (2012). The Human Connectome Project: a data acquisition perspective. Neuroimage 62, 2222–2231. doi: 10.1016/j.neuroimage.2012.02.018

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Vrba, J., and Robinson, S. E. (2001). Signal processing in magnetoencephalography. Methods 25, 249–271. doi: 10.1006/meth.2001.1238

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Weiss, S. A., Banks, G. P., McKhann, G. M. Jr., Goodman, R. R., Emerson, R. G., Trevelyan, A. J., et al. (2013). Ictal high frequency oscillations distinguish two types of seizure territories in humans. Brain 136, 3796–3808. doi: 10.1093/brain/awt276

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Worrell, G. A., Jerbi, K., Kobayashi, K., Lina, J. M., Zelmann, R., and Le Van Quyen, M. (2012). Recording and analysis techniques for high-frequency oscillations. Prog. Neurobiol. 98, 265–278. doi: 10.1016/j.pneurobio.2012.02.006

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Xiang, J., Degrauw, X., Korostenskaja, M., Korman, A. M., O'Brien, H. L., Kabbouche, M. A., et al. (2013). Altered cortical activation in adolescents with acute migraine: a magnetoencephalography study. J. Pain 14, 1553–1563. doi: 10.1016/j.jpain.2013.04.009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Xiang, J., Holowka, S., Qiao, H., Sun, B., Xiao, Z., Jiang, Y., et al. (2004). Automatic localization of epileptic zones using magnetoencephalography. Neurol. Clin. Neurophysiol. 2004, 98.

Pubmed Abstract | Pubmed Full Text

Xiang, J., Liu, Y., Wang, Y., Kirtman, E. G., Kotecha, R., Chen, Y., et al. (2009a). Frequency and spatial characteristics of high-frequency neuromagnetic signals in childhood epilepsy. Epileptic Disord. 11, 113–125. doi: 10.1684/epd.2009.0253

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Xiang, J., Liu, Y., Wang, Y., Kotecha, R., Kirtman, E. G., Chen, Y., et al. (2009b). Neuromagnetic correlates of developmental changes in endogenous high-frequency brain oscillations in children: a wavelet-based beamformer study. Brain Res. 1274, 28–39. doi: 10.1016/j.brainres.2009.03.068

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Xiang, J., Wang, Y., Chen, Y., Liu, Y., Kotecha, R., Huo, X., et al. (2010). Noninvasive localization of epileptogenic zones with ictal high-frequency neuromagnetic signals. J. Neurosurg. Pediatr. 5, 113–122. doi: 10.3171/2009.8.PEDS09345

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Xiang, J., and Xiao, Z. (2009). Spatiotemporal and frequency signatures of noun and verb processing: a wavelet-based beamformer study. J. Clin. Exp. Neuropsychol. 31, 648–657. doi: 10.1080/13803390802448651

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Zafeiriou, D. I., and Vargiami, E. (2012). Noninvasive ultra high-frequency (1kHz) oscillations' recording: high-fidelity over somatosensory cortex. Clin. Neurophysiol. 123, 2323–2324. doi: 10.1016/j.clinph.2012.05.010

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Zijlmans, M., Huiskamp, G. M., Cremer, O. L., Ferrier, C. H., Van Huffelen, A. C., and Leijten, F. S. (2012a). Epileptic high-frequency oscillations in intraoperative electrocorticography: the effect of propofol. Epilepsia 53, 1799–1809. doi: 10.1111/j.1528-1167.2012.03650.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Zijlmans, M., Jiruska, P., Zelmann, R., Leijten, F. S., Jefferys, J. G., and Gotman, J. (2012b). High-frequency oscillations as a new biomarker in epilepsy. Ann. Neurol. 71, 169–178. doi: 10.1002/ana.22548

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text