Estimation of nonlinear neural source interactions via sliced bicoherence

Biomedical Signal Processing and Control 30 (2016) 43–52

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

Estimation of nonlinear neural source interactions via sliced bicoherence Tolga Esat Özkurt Department of Health Informatics, Graduate School of Informatics, Middle East Technical University, Üniversiteler Mahallesi, 06800 C¸ankaya, Ankara, Turkey

a r t i c l e

i n f o

Article history: Received 6 September 2015 Received in revised form 25 December 2015 Accepted 5 May 2016 Keywords: Bicoherence Connectivity Cross-frequency coupling Higher-order spectrum Neural oscillations Synchronization

a b s t r a c t Neural oscillations and their spatiotemporal interactions are of interest for the description of brain mechanisms. This study offers a novel third order spectral coupling measure named “sliced bicoherence”. It is the diagonal slice of cross-bicoherence allowing an efficient quantification of the nonlinear interactions between neural sources. Our methodology comprises an indirect estimation method, a parametric confidence level formula and a subtracted version for robustness to volume conduction. The methodology provides an efficient estimation of third-order nonlinear cross relations reducing the complexity to the same order of second-order coherence computation. Unlike other bispectral measures, the suggested measure solely holds terms related to cross relations between channel sources and omits the possible strong autobispectral relations. Feasibility and robustness of the methodology are demonstrated both on simulated and publicly available MEG data. The latter were collected for a motor task and an eyes-open resting state. Analytical confidence level marked the non-significant couplings. Simulations confirmed that the subtracted bicoherence enabled robustness to volume conduction by avoiding the spurious nearby channel couplings. Central regions were shown to be coupled with muscular activity by sliced bicoherence. Couplings for spontaneous data occurred particularly at theta and alpha bands. Volume-conduction related bicoherence values originated especially from the low frequencies below 5 Hz. The suggested nonlinear measure is promising to be a part of the rich collection of the multichannel electrophysiological brain connectivity metrics. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction Neural oscillations have been assumed to play a fundamental role during the functioning of normal and pathological brain. Specifically, neural synchronization is known to be an essential means for the communication of spatially distant oscillatory networks in the brain [1,2]. Synchronization is achieved through interaction of massive neural populations either at one frequency or multiple frequencies. Coherence is a well-known measure to obtain linear phase and/or amplitude relations between two sites at one specific frequency. Recent years have also witnessed the ubiquitous use of cross-frequency measures [3] such as phase-amplitude coupling which quantifies low frequency phase and high frequency amplitude synchronization [4] and amplitude-amplitude correlation [5,6]. Various studies have shown that synchronization measures may signify genuine neural interactions characterizing

E-mail addresses: [email protected], [email protected] http://dx.doi.org/10.1016/j.bspc.2016.05.001 1746-8094/© 2016 Elsevier Ltd. All rights reserved.

the fundamental mechanisms of the brain for medical [7], taskrelated [8] and resting states [9]. Rapid estimation of neural synchronization measures is of utmost importance while revealing the brain networks from macro scale electrophysiological multichannel data. Thus, employed methods need not only yield computationally affordable estimations of synchronization measures but should also allow efficient determination of statistical significance. Examples for the latter utilize the simple parametric formulae assuming standard distributions (Gaussian, uniform etc.) for the assessment of coherence [10] and phase-amplitude coupling [11] confidence limits. There have been various studies using bispectrum and its normalized version named as bicoherence on EEG and MEG data [12–16]. We would like to emphasize the expensive computational cost as one of the main obstacles over the bispectral measures. A pairwise bispectral analysis of multichannel data requires estimates in the order of ∼N2 × M2 , where N and M denote the number of channels and the number of sampled frequencies, respectively. In order to alleviate the huge computational cost of bispectral analysis, Chella et al. [16] applied principal component analysis to reduce

44

T.E. Özkurt / Biomedical Signal Processing and Control 30 (2016) 43–52

the data dimension prior to the parameter estimation. However, this operation alone still cannot decrease the degree of the complexity but rather reduces the number of channels to the number of principal components. Current study presents an efficient indirect estimation of diagonal slice of cross-bicoherence (normalized cross-bispectrum), hence reducing the complexity to the same order of coherence computation. We also present simple formulae to obtain the confidence limits of the estimate. The suggested procedure allows feasible extractions of higher-order spectra related brain connectivity from data with high number of channels and low computational cost, making it computationally comparable to that of the classical coherence and the power spectrum. 2. Methods The 1,1,2 cross-bicoherence for two signals x1 (n) and x2 (n) is given by b12 (f1 , f2 ) =

|E{X1 (f1 ) X1 (f2 ) X2 ∗ (f1 + f2 )}| S1 (f1 ) S1 (f2 ) S2 (f1 + f2 )

2

(1)

where superscript * stands for the complex conjugate and E{} is the statistical expectation operator [17]. Here X and S denote Fourier coefficients and spectra respectively. Please note that the rest of the text shall consider a particular slice of bispectrum, i.e., the diagonal slice where f1 = f2 . Despite not being explicitly stated, higherorder couplings in this slice were commonly observed in empirical neuroelectrophysiological studies (see Section 4, for some examples). Nevertheless, in principle, it is straightforward to select any other slice of the cross-bispectrum to proceed with the suggested methodological approach. 2.1. Diagonally sliced 1,1,2 cross-bicoherence: definition and relation to coherence For the sake of brevity, the measure in the subtitle shall be called “sliced bicoherence” throughout this paper. We would like to consider sliced bicoherence in analogy with the well-known coupling measure of coherence, which is defined for two signals x1 (n) and x2 (n) as [18,19]:





2

2 |E X1 (f )X2 ∗ (f ) | |S12 (f )| = c12 (f ) = . S1 (f )S2 (f ) S1 (f )S2 (f )

(2)

Similar to the coherence, sliced bicoherence for the signals x1 (n) and x2 (n) can be formulated as: b12 (f ) =

|B112 (f, f )|



2

S1 2 (f )S2 (2f )

=



|E X1 2 (f )X2 ∗ (2f ) | S1 2 (f )S2 (2f )

2

(3)

where B denotes (cross) bispectrum. Please notice the structural similarity of coherence (Eq. (2)) and sliced bicoherence (Eq. (3)) formulations. Some essential distinctive properties may be listed as follows: 1. Coherence is a 2nd order cross-spectral measure. While sliced bicoherence is a 3rd order one. 2. Coherence is a linear measure. Input and output signals of a linear filter are also perfectly coherent. While sliced bicoherence is a nonlinear measure. 3. Unlike coherence, sliced bicoherence is not symmetric, i.e., b12 (f) = / b21 (f). As it is clear from Eq. (3) that b12 (f) implies the relation of the frequency component at f in the first signal to the component at 2f in the second one, while the other way around is the case for b21 (f). 4. Two signals are perfectly coherent if they are 1:1 linear phase locked at the same frequency f. While two signals are perfectly

sliced bicoherent if they are 1:2 quadratic phase locked at a frequency f for one signal and its double 2f for the other signal. 5. A signal is perfectly coherent with itself, i.e., c11 (f) = c22 (f) = 1 for ∀f. This is not the case for sliced bicoherence. 2.2. Estimation The IDFT of the numerator of Eq. (3) for “autobispectrum”, i.e., B111 (f, f ) was called as “sum-of-cumulants” (SOC) in signal processing [20]. It has been used for various purposes such as system identification [21], signal reconstruction [22], detection of nonlinearity [23] and robust speaker recognition [24]. It should be noted that these studies utilized SOC to capture only “auto” bispectral features in the aforementioned literature. The current study instead aims at seeking cross-frequency couplings via the cross-bispectrum related measures. In order to estimate b12 , one may compute the average over a number of segments, just like in coherence. Computing the whole total bispectrum to obtain the numerator of Eq. (3) would take considerable time, as its dimension is higher than the spectrum. Instead, we suggest using an indirect formula computing the crossbispectral slice:

where stands for Discrete Fourier Transform, * is the convolution operator and



y2 (n) =

x2 (N − 1 − n/2), n is even 0,

.

(5)

n is odd

A similar formula to Eq. (4) has been used for the computation of “auto” bispectrum in some studies [22,24]. The validity of Eq. (4) can be shown using the convolution theorem (see Appendix A). 2.3. Confidence limit Halliday et al. [10] suggested a plain parametric formula depending on the number of segments in order to compute a statistical limit for coherence. Analogously, Özkurt [11] derived a formula that gives a statistical limit for direct phase-amplitude coupling measure. The latter formula solely depends on data length. Both studies assumed standard and reasonable probability distributions (such as Gaussian and uniform) that are fairly easy to manipulate and hence obtain simple confidence limits. They enable efficient statistical thresholds for the coupling estimates. In a similar fashion, Haubrich [25] showed that the 0.95 confidence level for bicoherence is 3/N for signals with Gaussian probability distributions, where N denotes the number of segments. One should note that Gaussian signals ideally would have zero bicoherence. Hence the derived formula provides a lowest limit while deciding the reliability of the resultant bicoherence value. Elgar and Guza [26] compared this formula to the results of numerical simulations. We here show that Haubrich’s formula is only valid for the bicoherence b12 (f1 , f2 ) of different frequencies, that is for f1 = / f2 . Our numerical derivation indicated that specifically for the diagonal sliced bicoherence b12 (f ), the confidence limit becomes 6/ (N + 1) instead, that is for f1 = f2 . Please see Appendix B for the derivation and note that we additionally validated the confidence level derivation with a Monte Carlo simulation. 2.4. Use in electrophysiology One can reliably estimate the oscillatory coupling between a frequency of a signal in a channel and the double frequency in another

T.E. Özkurt / Biomedical Signal Processing and Control 30 (2016) 43–52

“independent” channel using Eq. (4). For example, muscular – cortical couplings or basal ganglia – cortical couplings [27] could be obtained from the given bispectral slice formula. However, using this cross-frequency measure over raw EEG or MEG channels may be problematic because of volume conduction. Volume conduction and magnetic field spread lead to common existent oscillatory components in the EEG/MEG channel signals. This in turn causes spurious coupling estimations between those channels. Although the degree of volume conduction may be reduced with source identification methods such as Dynamic Imaging of Coherent Sources (DICS) beamformer [28]; its hazardous effect still cannot be totally eliminated in the projected source space [29,30]. Regarding the volume conduction problem with respect to coherence, Nolte et al. [31] suggested the use of imaginary part of complex coherency. Imaginary coherency achieves the alleviation of the volume conduction effect by disregarding the zero-phase sources in the coherence estimation. The only price paid is henceforth the elimination of “real” zero-phase coherence as well, which may genuinely exist between brain sources. The effect of volume conduction on the cross-channel sliced bicoherence should be approached cautiously. Common source signals in two channels would lead to very high spurious “coherence” values between them. Since a signal has to be perfectly coherent with itself, i.e., c11 (f ) = 1 for ∀f . However, this is not the case for the bicoherence, i.e., b11 (f ) = / 1 holds in general. Hence, volume conduction leads to high values in bicoherence only if there exists an intrinsic nonlinear coupling in the common source signals shared by the channel pair. 2.5. Cross-bispectrum slice robust to volume conduction: subtracted bicoherence A recent study [16] asserted a bispectral metric that takes the antisymmetric part of cross-bispectrum in order to get rid of the volume conduction effect. This corresponds to |B112 − B211 | for the bispectral slice b12 (f). Although, the antisymmetric part may erase “autobispectra” of the sources, it still keeps the autobispectra of the channels themselves leaked into the cross-bispectra of the two distinct channels. Here is a simple example case of two channels x1 and x2 sharing a common source s2 : x1 (n) = s1 (n) + s2 (n)

(6)

x2 (n) = s2 (n) + s3 (n) Let us insert Eq.(6) into    B112 − B211 = E X1 2 (f )X2 ∗ (2f ) − E X2 (f )X1 (f )X1 ∗ (2f ) , to obtain a resultant comprising source terms Bs1s1s2 , Bs1s1s3 ,Bs2s2s3 , Bs1s2s2 , Bs1s2s3 (see Appendix C). Hence, the antisymmetric part eliminates Bs2s2s2 but still maintains terms Bs1s1s2 and Bs1s2s2 which originate from the autobispectra of x1 . This would mix autobispectral bias terms into the cross-channel coupling estimate. As a result, a cross channel bispectral relation can be overshadowed. Moreover, other mixed cross bispectral terms such as Bs1s2s1 and Bs3s2s2 have nothing to do with the cross relation x1 → x2 (1,1,2 relation), which we aim to estimate by using b12 (Eq. (3)) in this study. In order to alleviate this problem, we suggest an alternative metric of sliced bicoherence being robust to volume conduction. In addition, the suggested metric solely holds terms related to “cross” relations between channels and omits the possible strong autobispectral relations that are greater in the channel x1 itself:



b˜ 12 (f ) =

b12 (f ) − b11 (f ), if b12 (f ) − b11 (f ) > 0 0

, if b12 (f ) − b11 (f ) ≤ 0

(7)

45

If a spurious cross-bispectral relation occurs because of a common source, its effect is either totally eradicated (b12 (f ) − b11 (f ) ≤ 0) by Eq. (7) or reduced down to a level depending on the strength of that source in the first channel x1 (for the case b12 (f ) − b11 (f ) > 0). Hence it provides a “safe” cross-bispectral estimate by avoiding autobispectral contributions from the channels. This metric is called “subtracted bicoherence” throughout the following text. Please note that, except the robustness to volume conduction, the suggested subtracted bicoherence measure has a different aim than the antisymmetric bispectral measure suggested by [16]. The former (i) ignores any bispectral relation between channels stemming from the common sources and more importantly, (ii) accounts only for pairwise 1,1,2 bispectral relations and omits the others unlike the latter. Also note that if there are distinct sources belonging to the pairwise channels and the auto bispectral relation in Channel 1 occurs at the very same frequency with a cross-relation for Channel1 → Channel 2, then there is a chance that the latter is lessened or even removed through subtracted bicoherence measure. This potential “price” may be thought analogously to other volume-conduction immune measures such as imaginary coherency [31], which is blind to source coherences occurring at zero phase. 3. Experiments 3.1. Simulation for coherence vs. sliced bicoherence We simulated two signals having similar spectra such that they both comprise harmonic components at 5 Hz, 10 Hz and 20 Hz (Fig. 1A). The sampling frequency was chosen as 256 Hz. We generated them such that the couplings are between 5–5 Hz and 10–20 Hz. While coherence is better suited to capture the former, the sliced bicoherence is supposedly suitable for the latter. In order to obtain the simulated harmonics, white Gaussian noise was passed through two-way least-squares finite impulse response filters centered at around the harmonic frequencies. The filter bandwidths were made 1 Hz for the harmonic frequencies 5 Hz and 10 Hz, while the bandwidth was 2 Hz for the harmonic at 20 Hz. The supposed linear phase relation at 5 Hz was produced by time delaying one of the corresponding harmonics. The coupling relation at 10–20 Hz was produced simply by squaring the lower frequency harmonic as this is known to cause a quadratic phase coupling [32]. We added extra white non-Gaussian noise (exponential noise) to both signals at the level of SNR = 0 dB, i.e., the simulated signal and noise possess the same amount of energy. We removed the mean value from the final signal. The simulation procedure was repeated 100 times. Coherence and sliced bicoherence estimations were realized with non-overlapped windows of length 128, which corresponds to the estimator frequency resolution limit of 2 Hz. We used the Matlab® (Mathworks Inc., Natick, MA) built-in function mscohere for the coherence and applied our own routine implemented according to the “indirect” estimator described in the previous section for the sliced bicoherence. The Matlab routines for the sliced and subtracted bicoherence estimation were made available in the Supplementary material. The coupling relations between the simulated signals is given in Fig. 1A. Coherence and sliced bicoherence measures identify the assigned relations correctly as shown by Fig. 1B and C. 3.2. Simulation for cross-channel bispectrum Another simulation was realized in order to demonstrate the cross-bispectral (sliced, subtracted and antisymmetric) measures between channels containing mixtures of sources. In this direc-

46

T.E. Özkurt / Biomedical Signal Processing and Control 30 (2016) 43–52

Fig. 1. Simulated signals s1 and s2 (A) both have spectral peaks at 5, 10 and 20 Hz and they are coupled at 5 Hz (linear time shifting) and 10 Hz (quadratic phase coupling). (B) Coherence captures the linear phase relation at 5 Hz, while (C) sliced bicoherence identifies the nonlinear phase relation at 10 Hz. Simulation was repeated 100 times. Shaded areas indicate standard error of the mean.

Fig. 2. (A) Simulated sources (s1 ,s2 ,s3 )—channels (y1 ,y2 ) and (B) the produced bispectral coupling (1,1,2) relations between sources. For instance, an arrow directed from s1 to s2 denotes the simulated relation b12 (f) at the indicated frequency f. It was defined in Eq. (2). Channels y1 and y2 contain mixtures of sources weighted according to Sarvas formula. Source and channel properties are provided in Table 1.

Table 1 Simulated dipole source and channel parameters. Sources and channels

Locations (cm)

Orientations

s1 s2 s3 y1 y2

(0,−2,8.5) (0,0,9) (0,2,9.5) (0,−2,11) (0,2,11)

(2,1,1) (1,2,1) (1,1,2) (1,1,1) (1,1,1)

tion, we produced 3 sources and they were mapped to 2 channels (Table 1 and Fig. 2A) via Sarvas formula [33]. The channel y1 was positioned such that it mainly picks s1 but very little of s3 , as the latter is located far away from y1 . On the contrary, the channel y2 mainly picks s3 and ignorable contribution from s1 . The common source s2 was collected by both channels with comparable amounts of weights.

Quadratic phase coupling relations between sources were produced as depicted in Fig. 2B. The cross-channel relation is supposed to be only for y1 → y2 at 15 Hz. All other relations are supposedly irrelevant. This could only be revealed by the subtracted bicoherence measure (Fig. 3A; top). As expected, the sliced bicoherence yielded additional spurious peaks at 25 Hz and 35 Hz exceeding the confidence level (Fig. 3A; middle). These stem from spurious common source s2 auto bispectral relation at 35 Hz and the bispectral relation from s1 to s2 at 25 Hz. Another measure suggested by Chella et al. [16] namely “antisymmetric bispectrum” could eradicate the common source auto bispectral relation at 35 Hz. However, it still kept the peak at 25 Hz indicating the relation from s1 to s2 (Fig. 3A; bottom). We also investigated the cross-channel inversely directed relation y2 → y1 . In fact, according to our designed simulation, there is no relation from s3 (reflected by y2 ) to s1 (reflected by y1 ). Accord-

T.E. Özkurt / Biomedical Signal Processing and Control 30 (2016) 43–52

47

Fig. 3. The robustness of subtracted bicoherence (top row) measure when compared with “sliced” bicoherence (middle row) and “antisymmetric” (bottom row) bispectrum measures. Simulated source properties were given by Fig. 2. (A) shows the couplings directed from Channel 1 → Channel 2, while (B) shows the estimated couplings in the inverted direction Channel 2 → Channel 1. Only subtracted coherence estimate keeps the cross-channel relation at 15 Hz from Channel 1 to Channel 2 (above the confidence level threshold); while the other measures include spurious or irrelevant coupling estimates.

ingly, the proposed measure, subtracted bicoherence did not reveal any peaks over the confidence level (Fig. 3B; top). However, the other measures showed spurious peaks (Fig. 3B; middle and bottom). These results confirm that the subtracted bicoherence can successfully reveal cross-channel bispectral relation b12 reflecting cross-source relations distinctly belonging to those channels and being robust to volume conduction as theoretically explained in Section 2.5.

3.3. Sliced bicoherence with an external reference channel We conducted coupling analysis on simultaneously recorded MEG–EMG data provided by Fieldtrip toolbox [34]. The subject was instructed to raise her left hand while data were acquired. EMG data were reported to be obtained from the extensor carpi radialis longus muscle in the arms. MEG data were recorded via 151 channel CTF Omega System (Port Coquitlam, Canada) with the sampling

48

T.E. Özkurt / Biomedical Signal Processing and Control 30 (2016) 43–52

Fig. 4. Coupling relations between MEG channels and an EMG channel while a subject raises his left hand. (A) Corticomuscular coherence being prominent at beta band (15–30 Hz) is contralateral to the movement and strongest in channels located above the primary motor cortex. (B) Corticomuscular sliced bicoherence at alpha–beta bands (6–30 Hz) was more prominent on central channels and (C) it is significant at 6 Hz with its two first harmonics. The dashed line indicates the analytically found 95% confidence level.

rate of 1200 Hz. All other details regarding data acquisition may be found at the Fieldtrip site tutorial (http://www.fieldtriptoolbox. org). It is common practice to use coherence as a synchronization measure for the quantification of corticomuscular activity, especially for neuromuscular disorders such as Parkinson’s disease [27]. The coherent site is expected to be the primary motor cortex at the contralateral site. The coherence between the left EMG and MEG channels at the beta band (15–30 Hz) can be observed in Fig. 4A. The topographic plot agrees well with the established coherent region of primary motor cortex. Similarly, we applied “sliced bicoherence” measure to reveal the possible cortical sites coupled with EMG activity on the left arm. The segment length was taken as 0.25 s. The couplings were found to be significant at around 6 Hz and its first two harmonics at 12 Hz and 24 Hz (Fig. 4C). The topographic plot in Fig. 4B shows the coupled site (6–30 Hz), which seems more central when compared to the right hemisphere lateralized coherent region depicted in Fig. 4A. The cerebral network responsible for motor activity comprises various cortical and subcortical structures [35]. Hence, it is reasonable that except the peripheral regions, more central regions such as supplementary motor areas (SMA) are also involved in corticomuscular synchronization as evidenced by sliced bicoherent couplings.

3.4. Channel level sliced and subtracted bicoherence We used MEG data released by Human Connectome Project (HCP; MEG2 release), which included 67 subjects [36]. The MEG system was a whole-head Magnes 3000 (4D Neuroimaging, San Diego, CA) with 248 magnetometer channels. Subjects were instructed to keep their eyes open in a darkened room in supine position during acquisition. The data length amounted approximately to 18 min in 3 epochs, each of which took about 6 min. A subset of resting data comprising arbitrarily chosen 9 subjects was under sliced bispectral analysis in our study. We used the readily channel level Fieldtrip preprocessed data also included in the same dataset. The raw data were passed through Independent Component Analysis in order to eliminate noisy channels and data segments of 2 s length. The sampling frequency for the preprocessed data was 508.63 Hz. The details of acquisition and preprocessing pipeline can be found in “WU-Minn HCP 500 Subjects + MEG2 Data Release: Reference Manual” (http://www. humanconnectome.org/documentation/S500/HCP S500+MEG2 Release Reference Manual.pdf). Both sliced and subtracted bicoherence measures were applied on each epoch for each subject. We used the first 200 s of each epoch for the estimation. The segment length was taken as 512 which corresponds to ∼1 s. Henceforth the confidence level becomes 0.03.

T.E. Özkurt / Biomedical Signal Processing and Control 30 (2016) 43–52

49

Fig. 5. Coupling topography for a representative subject. Siginificant bispectral couplings with the medial frontal channel at theta (epoch 1, epoch 3) and alpha (epoch 2) band are ubiqitious at resting state. Subtracted bicoherence measure alleviates volume conduction related spurious estimated couplings by sliced bicoherence.

A reference seed channel corresponding to central prefrontal cortex (channel label: A37) was selected. Couplings between the reference channel and remaining channels were assessed. Please note that we did not aim for a thorough bispectral coupling analysis for this study but wanted to demonstrate the feasibility of the suggested methodology. We observed sliced bispectral couplings particularly at theta and alpha bands for all subjects at various channels corresponding to temporal, parietal and occipital cortices. Importantly, the coupling locations and strengths seem to change temporally reflecting spatial dynamic nonlinear relationships. Assessed bispectral channel based topography for an exemplar subject is exhibited in Fig. 5 for all three epochs. The topographies were depicted for either theta (Fig. 5 left and right; 4–7 Hz) and alpha (Fig. 5 middle; 8–12 Hz) bands. For illustration, we chose the band that had the greater sliced bicoherence peak. Aforementioned temporal evolution of the coupled channel locations can also be observed from the figures. The coupling masked by common source autobispectra nearby channels to the seed could be revealed thanks to the subtracted bispectrum measure. This can especially be realized from the figures for epoch 1, where the spurious autobispectral source contribution is eliminated and medial parieto-temporal cross-relations are accentuated (Fig. 5, top left). The same spot is less obvious for sliced bicoherence measure without subtraction (Fig. 5, bottom left). We also applied the antisymmetric bispectral measure to compare with the suggested measures. It could eliminate volume conduction effect and gave comparable results to the subtracted bicoherence, however it omitted also 1,1,2 couplings quite away from the reference frontal channel (see Supplementary Fig. S1). Fig. 6 shows the subtracted and sliced bicoherence values with respect to the frequency for the same subject and epochs. Selected channel locations were indicated by circles. We observe that the spurious volume-conduction related bicoherence values stem especially from the low frequencies below 5 Hz. The determined confidence level indicated prominence of alpha and theta couplings as the dashed line successfully identifies the non-significant bicoherences for frequencies outside the range of alpha and theta. Similar results were also obtained for the other subjects (see Supplementary Figs. S2 and S3).

4. Discussion This study suggests diagonally “sliced bispectrum” to be added to the rich collection of the electrophysiological brain connectivity

metrics. In this direction, we proposed an efficient time-domain estimation method, an analytical confidence level formula and a “subtracted” version to improve robustness to volume conduction. Feasibility of the suggested methodology was demonstrated via both simulations and real MEG data. It is straightforward to apply our methodology to identify the nonlinear interactions in source level with the addition of an inverse source reconstruction method. Bispectral phase has been exploited for the estimation of timedelay [37], phase difference angle [16] or statistical properties of the time-domain signal [38]. Our study is restricted specifically to bispectral “magnitude”, which is to capture the coupling strength between brain source sites. Moreover, the suggested measure focuses on spatially distanced cross-source relations and hence excludes the auto couplings intrinsic to one local site. Interestingly, many human and animal brain studies using bispectral measures could detect conspicuous neural activity particularly on the diagonal slice, even though they investigated the whole bispectra. For instance, Wang et al. [15] could find out the only significant bispectral power at the (16–16 Hz) slice. The data were of a monkey LFP obtained while a visual Go/No-go task was realized. Ning and Bronzino [39] showed bispectral peaks at around (7 Hz, 7 Hz) at the hippocampus and at around (3 Hz, 3 Hz) at the frontal cortex of rats for the REM and quiet waking states, respectively. Saikia and Hazarika [40] computed EEG bispectra for various tasks regarding the imagination of hand movement. All apparent bispectral peaks were reported to occur exclusively on the diagonal slice. Bispectral couplings between spontaneous EEG data channels of nine subjects were inspected by [16]. The only significant common peaks they could identify was at (10 Hz, 10 Hz) and (12 Hz, 12 Hz). This was interpreted as a consequence of an interaction between alpha and beta bands in the resting state. A similar result as the bispectral coupling between these bands was also reported for eyes-closed data [41]. Furthermore, various phase synchronization studies detected phase related couplings between a frequency and its double, which could be reflected on bispectral diagonal slice. Tass et al. [42] identified a cross-frequency 1:2 phase synchronization between EMG and MEG data for a tremor dominant Parkinson disease patient. They showed that sensorimotor and premotor cortical activities at around 10 Hz were coupled to the muscular activity at 5 Hz. Schack et al. [43] found out 1:2 phase synchronization between posterior and anterior EEG channels during a working memory task. Frequencies of the coupled oscillations were 6 Hz and 12 Hz. The coupling strength of these oscillations was reported being varied according to the memory load.

50

T.E. Özkurt / Biomedical Signal Processing and Control 30 (2016) 43–52

Fig. 6. Bicoherence spectra for the same data given in Fig. 3 with selected coupled channels marked as red circles. The spurious couplings are eliminated with subtracted bicoherence measure especially at low frequencies. Remaining couplings happen to be significant in theta–alpha range. The dashed lines indicate 95% confidence level. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

These studies witness the ubiquity of couplings that could be revealed by the suggested measure in this study. It should be emphasized that a vast majority of higher-order spectra based studies aimed at estimating bispectrum in its all three dimensions for all pairwise frequency combinations, i.e., for f1 being not necessarily equivalent to f2 . However, in many cases, the prominent peaks could mostly be detected at the diagonal slice f1 = f2 for the real data. Surely, this is by no means to imply ignoring any potential benefits that could be obtained from the remaining regions of bispectrum, which is an open hall for future studies. In any case, results of the aforementioned real data studies underline the potential usefulness of the diagonally sliced bicoherence measure in neuroelectrophysiology. Accordingly, our study allows an efficient estimation of this measure supplemented with an analytic confidence level formula and volume conduction avoidance, which make it a promising methodology while exploring nonlinear connectivity in brain.

(#112E562). He would like to thank Natalia Melnik for her help to download MEG data. He also thanks Dr. Hüseyin Hacıhabibo˘glu for proofreading the manuscript. Data were provided in part by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. Appendix A. According to the convolution theorem:

where and denote discrete Fourier and inverse Fourier transforms, respectively. Let xu (n) be upsampled x2 (n) by 2:



Acknowledgements The author was supported in part by the Scientific and Technological Research Council of Turkey under TÜBI˙ TAK 3501

xu (n) =

x2 (n/2), n even 0,

n odd

.

(A1.2)

T.E. Özkurt / Biomedical Signal Processing and Control 30 (2016) 43–52

51

Then it’s easy to show that x2 (n) = s2 (n) + s3 (n)

Finally for a discrete causal signal with length N, Eq. (A13) turns into

For the sake of brevity, we assume the source strength to be equal in the channels. Then one can easily express the cross channel bispectrum in terms of source bispectra:





B112 = E [s1 (f ) + s2 (f )][s1 (f ) + s2 (f )][s2 ∗ (2f ) + s3 ∗ (2f )]

= Bs1s1s2 + Bs1s1s3 + Bs2s2s2 + Bs2s2s3 + 2Bs1s2s2 + 2Bs1s2s3 Appendix B.



Haubrich [25] derived a confidence level for zero bicoherence / f2 . Here, we demonstrate a similar derivation for the sliced for f1 = bicoherence where f1 = f2 . Assume that x1 (n) and x2 (n) are two zero mean white Gaussian random processes and that they are divided into N number of statistically independent segments for the estimation of the sliced bicoherence b12 . Thus, the Fourier coefficients akj + ibkj for a frequency k are also zero mean Gaussian and let their variance be  akj 2 =  bkj 2 = r2 . N 

1 | N2

bˆ 12 =

 1 N3

(akj 2 + i2akj bkj − bkj ) × (amj − ibmj )|

2

N

2

(akp 2 + bkp )

 ×

p=1

N  

amq 2 + bmq

2





(A2.1)

=



2

(akj 2 amj − bkj amj + 2akj bkj bmj )

 2  4 × E akj amj 2 + E N

×E



2

akj 2 bkj bmj

2





4

bkj amj 2

− 2×E

=

2 × (3r 6 + 3r 6 + 4r 6 − 2r 6 ) N

=

16r 6 . N

⎧ N ⎨ 

2



2

2 ⎫ ⎬

(akp + bkp )

p=1



2



+4 2

akj 2 bkj amj 2

×E

⎭

2

=

1 (4N + 4N 2 )r 4 × 2Nr N3

=

8r 6 (N + 1) . N



(A2.2)

 N 

2

amq + bmq

2





q=1

(A2.3)

By taking the ratio using Eqs. (A2.2) and (A2.3), we finally show 2 that bˆ 12 is ␹2 distributed with mean E(bˆ 12 ) = N+1 . This would correspond to

6 N+1

−(Bs1s2s1 + Bs2s1s2 + Bs2s2s1 + Bs3s1s1 + Bs1s3s2 + Bs2s3s1 +Bs3s2s2 ) Thus the subtraction eliminates the common source autobispectrum Bs2s2s2 but the channel autobispectra Bs1s1s2 , Bs1s2s2 , Bs1s2s1 and Bs2s2s1 are left embedded in the antisymmetric part. Appendix A. Supplementary data

Similarly one can work out the terms in the denominator and obtain 1 E N3 ⎩

B112 − B211 = Bs1s1s2 + Bs1s1s3 + Bs2s2s3 + 2Bs1s2s2 + 2Bs1s2s3

q=1

The expectation E(bˆ 12 ) can be computed by performing it separately for numerator and denominator of Eq. (A2.1). When multiplications inside the numerator are realized and all terms with zero expectations are omitted, one obtains: 2 ×E N

= Bs2s1s1 + Bs2s1s2 + Bs2s2s1 + Bs2s2s2 + Bs3s1s1 + Bs1s3s2 +Bs2s3s1 + Bs3s2s2

2 2

j=1





B211 = E [s2 (f ) + s3 (f )][s1 (f ) + s2 (f )][s1 ∗ (2f ) + s2 ∗ (2f )]

as 95% confidence level for zero sliced bicoherence.

Appendix C. Here we zoom into the antisymmetric part of “sliced” bicoherence B112 –B211 as suggested by Chella et al. [16] for two channels sharing a common source s2 : x1 (n) = s1 (n) + s2 (n)

Supplementary data associated with this article can be found, in the online version, at 10.1016/j.bspc.2016.05.001. References [1] A. Schnitzler, J. Gross, Normal and pathological oscillatory communication in the brain, Nat. Rev. Neurosci. 6 (4) (2005) 285–296, http://dx.doi.org/10.1038/ nrn1650. [2] M. Siegel, T.H. Donner, A.K. Enge, Spectral fingerprints of large-scale neuronal interactions, Nat. Rev. Neurosci. 13 (2012) 121–134, http://dx.doi.org/10. 1038/nrn3137. [3] V. Jirsa, V. Müller, Cross-frequency coupling in real and virtual brain networks, Front. Comput. Neurosci. 7 (78) (2013), http://dx.doi.org/10.3389/ fncom.2013.00078. [4] T.E. Özkurt, A. Schnitzler, A critical note on the definition of phase-amplitude cross-frequency coupling, J. Neurosci. Methods 201 (2) (2011) 6–11, http://dx. doi.org/10.1016/j.jneumeth.2011.08.014. [5] A. Bruns, R. Eckhorn, H. Jokeit, A. Ebner, Amplitude envelope correlation detects coupling among incoherent brain signals, NeuroReport 11 (7) (2000) 1509–1514. [6] J.F. Hipp, D.J. Hawellek, M. Corbetta, M. Siegel, A.K. Engel, Large-scale cortical correlation structure of spontaneous oscillatory activity, Nat. Neurosci. 15 (6) (2012) 884–890, http://dx.doi.org/10.1038/nn.3101. [7] C. de Hemptinne, E.S. Ryapolova-Webb, E.L. Air, et al., Exaggerated phase-amplitude coupling in the primary motor cortex in Parkinson disease, Proc. Natl. Acad. Sci. U. S. A. 110 (12) (2013) 4780–4785, http://dx.doi.org/10. 1073/pnas.1214546110. [8] J.F. Hipp, A.K. Engel, M. Siegel, Oscillatory synchronization in large-scale cortical networks predicts perception, Neuron 69 (2011) 387–396. [9] E. Florin, S. Baillet, The brain’s resting-state activity is shaped by synchronized cross-frequency coupling of oscillatory neural activity, NeuroImage 111 (2015) 26–35, http://dx.doi.org/10.1016/j.neuroimage.2015.01.054. [10] D.M. Halliday, J.M. Rosenburg, A.M. Amjad, P. Breeze, B.A. Conway, S.F. Farmer, A framework for the analysis of mixed time series/point process data—theory and application to the study of physiological tremor single motor unit discharges and electromyograms, Prog. Biophys. Mol. Biol. 64 (1995) 237–278. [11] T.E. Özkurt, Statistically reliable and fast direct estimation of phase-amplitude cross-frequency coupling, IEEE Trans. Biomed. Eng. 59 (7) (2012) 1943–1950. [12] T.P. Barnett, L.C. Johnson, P. Naitoh, N. Hicks, C. Nute, Bispectrum analysis of electroencephalogram signals during waking and sleeping, Science 172 (3981) (1971) 401–402. [13] G. Dumermuth, P.J. Hubera, B. Kleiner, T. Gassera, Analysis of the interrelations between frequency bands of the EEG by means of the bispectrum: a preliminary study, Electroencephalogr. Clin. Neurophysiol. 31 (2) (1971) 137–148. [14] T.H. Bullock, J.Z. Achimowicz, R.B. Duckrow, S.S. Spencer, V.J. Iragui-Madoz, Bicoherence of intracranial EEG in sleep, wakefulness and seizure, Electroencephalogr. Clin. Neurophysiol. 103 (6) (1997) 661–678.

52

T.E. Özkurt / Biomedical Signal Processing and Control 30 (2016) 43–52

[15] X. Wang, Y. Chen, M. Ding, Testing for statistical significance in bispectra: a surrogate data approach and application to neuroscience, IEEE Trans. Biomed. Eng. 54 (11) (2007) 1974–1982. [16] F. Chella, L. Marzetti, V. Pizzella, F. Zappasodi, G. Nolte, Third order spectral analysis robust to mixing artifacts for mapping cross-frequency interactions in EEG/MEG, NeuroImage 91 (2014) 146–161, http://dx.doi.org/10.1016/j. neuroimage.2013.12.064. [17] S. Elgar, J.W. van Atta, M. Gharib, Cross-bispectral analysis of a vibrating cylinder and its wake in low Reynolds number flow, J. Fluids Struct. 4 (1) (1990) 59–71. [18] J.S. Bendat, A.G. Piersol, Random Data: Analysis and Measurement Procedures, Wiley, New York, 1971. [19] A. Bruns, Fourier-, Hilbert- and wavelet-based signal analysis: are they really different approaches? J. Neurosci. Methods 137 (2) (2004) 321–332, http://dx. doi.org/10.1016/j.jneumeth.2004.03.002. [20] L. Luo, L.F. Chaparro, Parametric identification of systems using a frequency slice of the bispectrum, Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing 5 (1991) 3481–3484. [21] S.A. Alshebeili, A.E. C¸etin, A.N. Venetsanopoulos, An adaptive system identification method based on bispectrum, IEEE Trans. Circuits Syst. 38 (8) (1991) 967–969. [22] T. Akgül, A. El-Jaroudi, Reconstruction of mixed-phase signals from sum-of-autotriplecorrelations using least squares, IEEE Trans. Signal Process. 46 (1) (1998) 250–254. [23] T. Akgül, M. Sun, R.J. Sclabassi, A.E. C¸etin, Characterization of sleep spindles using higher order statistics and spectra, IEEE Trans. Biomed. Eng. 47 (8) (2000) 997–1009. [24] T.E. Özkurt, T. Akgül, Robust text-independent speaker identification using bispectrum slice, Proceedings of the IEEE Signal Processing and Communications Applications Conference (2004) 418–421, http://dx.doi.org/ 10.1109/78.651230. [25] R. Haubrich, Earth noises 5 to 500 millicycles per second, J. Geophys. Res. 70 (1965) 1415–1427. [26] S. Elgar, R.T. Guza, Statistics of bicoherence, IEEE Trans. Acoust. Speech Signal Process. 36 (10) (1988) 1667–1668. [27] J. Hirschmann, C.J. Hartmann, M. Butz, N. Hoogenboom, T.E. Özkurt, S. Elben, J. Vesper, L. Wojtecki, A. Schnitzler, A direct relationship between oscillatory subthalamic nucleus-cortex coupling and rest tremor in Parkinson’s disease, Brain 136 (Pt. 12) (2013) 3659–3670, http://dx.doi.org/10.1093/brain/awt271. [28] J. Gross, J. Kujala, M. Hamalainen, L. Timmermann, A. Schnitzler, R. Salmelin, Dynamic imaging of coherent sources: studying neural interactions in the human brain, Proc. Natl. Acad. Sci. U. S. A. 98 (2) (2001) 694–699, http://dx. doi.org/10.1073/pnas.98.2.694. [29] K. Sekihara, J.P. Owen, S. Trisno, S.S. Nagarajan, Removal of spurious coherence in MEG source-space coherence analysis, IEEE Trans. Biomed. Eng. 58 (11) (2011) 3121–3129, http://dx.doi.org/10.1109/TBME.2011.2162514. [30] M. Drakesmith, W. El-Deredy, S. Welbourne, Reconstructing coherent networks from electroencephalography and magnetoencephalography with

[31]

[32] [33] [34]

[35]

[36]

[37] [38]

[39] [40]

[41]

[42]

[43]

reduced contamination from volume conduction or magnetic field spread, PLoS One 8 (12) (2013) e81553, http://dx.doi.org/10.1371/journal.pone. 0081553. G. Nolte, O. Bai, L. Wheaton, Z. Mari, S. Vorbach, M. Hallett, Identifying true brain interaction from EEG data using the imaginary part of coherency, Clin. Neurophysiol. 115 (10) (2004) 2292–2307, http://dx.doi.org/10.1016/j.clinph. 2004.04.029. C.L. Nikias, J.M. Mendel, Signal processing with higher-order spectra, IEEE Signal. Process. Mag. 10 (3) (1993) 10–37. J. Sarvas, Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem, Phys. Med. Biol. 32 (1) (1987) 11–22. R. Oostenveld, P. Fries, E. Maris, J.M. Schoffelen, FieldTrip: open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data, Comput. Intell. Neurosci. 2011 (2011), http://dx.doi.org/10.1155/2011/ 156869. L. Timmermann, J. Gross, M. Dirks, J. Volkmann, H.J. Freund, A. Schnitzler, The cerebral oscillatory network of parkinsonian resting tremor, Brain 126 (1) (2002) 199–212, http://dx.doi.org/10.1093/brain/awg022. L.J. Larson-Prior, R. Oostenveld, S. Della Penna, et al., Adding dynamics to the Human Connectome Project with MEG, NeuroImage 80 (2013) 190–201, http://dx.doi.org/10.1016/j.neuroimage.2013.05.056. M.J. Hinich, G.R. Wilson, Time delay estimation using the cross bispectrum, IEEE Trans. Signal Process. 40 (1) (1992) 106–113. T.J. Maccarone, The biphase explained: understanding the asymmetries in coupled Fourier components of astronomical time series, Mon. Not. R. Astron. Soc. 435 (4) (2013) 3547–3558. T.K. Ning, J.D. Bronzino, Bispectral analysis of the rat EEG during various vigilance states, IEEE Trans. Biomed. Eng. 36 (4) (1989) 497–499. A. Saikia, S. Hazarika, Bispectrum analysis of EEG in estimation of hand movement, advances in computing and communications SE-12 Communications in Computer and Information Science, vol. 191, Springer, Berlin, Heidelberg, 2011, pp. 109–118, http://dx.doi.org/10.1007/978-3-64222714-1 12. F. Shahbazi, A. Ewald, G. Nolte, Univariate normalization of bispectrum using Hölder’s inequality, J. Neurosci. Methods 233 (2014) 177–186, http://dx.doi. org/10.1016/j.jneumeth.2014.05.030. P. Tass, M.G. Rosenblum, J. Weule, J. Kurths, A. Pikovsky, J. Volkmann, A. Schnitzler, H.J. Freund, Detection of n:m phase locking from noisy data: application to magnetoencephalography, Phys. Rev. Lett. 81 (15) (1998) 3291–3294. B. Schack, W. Klimesch, P. Sauseng, Phase synchronization between theta and upper alpha oscillations in a working memory task, Int. J. Psychophysiol. 57 (2) (2005) 105–114, http://dx.doi.org/10.1016/j.ijpsycho.2005.03.016.