EEG Bands of Wakeful Rest, Slow-Wave and Rapid-Eye-Movement Sleep at Different Brain Areas in Rats

Accumulating evidence reveals that neuronal oscillations with various frequency bands in the brain have different physiological functions. However, the frequency band divisions in rats were typically based on empirical spectral distribution from limited channels information. In the present study, functionally relevant frequency bands across vigilance states and brain regions were identified using factor analysis based on 9 channels EEG signals recorded from multiple brain areas in rats. We found that frequency band divisions varied both across vigilance states and brain regions. In particular, theta oscillations during REM sleep were subdivided into two bands, 5–7 and 8–11 Hz corresponding to the tonic and phasic stages, respectively. The spindle activities of SWS were different along the anterior-posterior axis, lower oscillations (~16 Hz) in frontal regions and higher in parietal (~21 Hz). The delta and theta activities co-varied in the visual and auditory cortex during wakeful rest. In addition, power spectra of beta oscillations were significantly decreased in association cortex during REM sleep compared with wakeful rest. These results provide us some new insights into understand the brain oscillations across vigilance states, and also indicate that the spatial factor should not be ignored when considering the frequency band divisions in rats.


INTRODUCTION
Neuronal oscillations in cortical networks are supposed to be essential for information communication, thereby underlying the fundamental brain functions (Basar et al., 2001;Buzsaki and Draguhn, 2004). What's more, it has been recognized that neural oscillations in the brain are often separated into several oscillatory bands, and each frequency band may originate in different cortical and subcortical structures, such as delta and spindle oscillations originate from the thalamo-cortical network (Dossi et al., 1992;Steriade et al., 1993), while theta oscillations are the most prominent in septo-hippocampal system (Nunez et al., 1987;Leung and Yim, 1991). Accumulating evidence suggests that specific brain states related to kinds of cognitive functions, such as perception, attention, and memory (Klimesch, 1999;Basar et al., 2001;Thut and Miniussi, 2009), are tightly associated with these oscillatory activities (Klimesch, 1999;Engel et al., 2001;Nicolelis et al., 2008).
In humans, the electroencephalogram (EEG) is widely used in the physiological and pathological studies (Babiloni et al., 2007;Basar and Güntekin, 2008;Lomas et al., 2015), in which EEG is often manually subdivided into broad frequency bands through visual inspection (Nunez and Cutillo, 1995). Similarly, in animals, frequency bands are always divided arbitrarily and determined empirically, based on the presumed homologous frequency bands in humans (Lancel and Kerkhof, 1989;Trampus et al., 1990;Grasing and Szeto, 1992). To our knowledge, only few studies determine the relevant functional frequency bands with multivariate statistical methods in humans (Klimesch, 1999;Corsi-Cabrera et al., 2000), rats (Corsi-Cabrera et al., 2001) and frogs (Fang et al., 2012). Actually, previous studies show that the frequency bands acquired from principal component analysis (PCA) are more consistent with the underlying physiological mechanisms in rats (Corsi-Cabrera et al., 2001). Note that electrophysiological oscillations are sensitive to recording sites (Siapas and Wilson, 1998;Klimesch, 1999;Sirota et al., 2003), indicating the possible involvement of different brain oscillatory generators (Mitchell et al., 2008;Young and McNaughton, 2009;Fang et al., 2010;Timofeev and Chauvette, 2013). Therefore, considering the spatial factor in the division of frequency bands would be more fulfilling (Klimesch, 1999).
Slow waves and spindle activities are the main landmarks of slow-wave sleep (SWS) (Steriade et al., 1993). In rats, studies of slow oscillations have been frequently restricted to 1-4 Hz, while frequency ranges of sleep spindles are shown in a wider distribution, from 6-12 to 14-20 Hz (Siapas and Wilson, 1998;Eschenko et al., 2006;Mölle et al., 2006;Johnson et al., 2010;Peyrache et al., 2011). During wakefulness and rapid-eyemovement sleep (REMS), a series of neurotransmitter systems abolish sleep spindles and slow waves (Steriade et al., 1993), which are replaced by "desynchronized" low-amplitude activities of beta and gamma from 10-15 to 80-120 Hz (Corsi-Cabrera et al., 2001;Brown et al., 2012). Unlike the division of theta and alpha bands in humans, alpha band is replaced by a wider range of the theta band (usually 4-12 Hz) in rats (Corsi-Cabrera et al., 2001;Buzsáki, 2006).
Factor analysis, a multivariate statistical method, facilitates identification of electrophysiological patterns neglected by visual inspection (Fang et al., 2012). It was hypothesized that the bands extracted by factor analysis would fit better the functional oscillations within the specific brain regions during specific behavior states. With this purpose, we concentrated on several brain regions of rats, using the factor analysis method, to extract the covariance of different frequencies across vigilance states, including wakeful rest (WR), SWS, and REMS. The results showed that the band divisions differ across vigilance states and brain regions. Unlike the traditional frequency band divisions, which here are more in accordance with the potential functional oscillatory activities.

Animals
Thirteen male Wistar rats (weighting 260-290 g) were used in this experiment. Before the surgery, rats were housed in small groups with food and water ad libitum and were maintained on a 12 h light/dark cycle (white lights on at 8:00). All experiments were approved by the Ethical Committee on Animal Experimentation of the University of Electronic Science and Technology of China (UESTC).

Surgery
Chronic electrode implantation was performed under general anesthesia (sodium pentobarbital 60 mg/kg body weight, i.p.), complemented with 0.6 ml atropine sulfate (0.5 mg/ml, s.c.) to prevent excessive secretion of the respiratory tract. Additional pentobarbital (15 mg/kg) was given intraperitoneally when required. Before and after the resection of the temporal muscle, local analgesia was administered with lignocaine (2%). All stereotactic coordinates were relative to bregma with the skull surface flat, according to Paxinos and Watson (2005). Each rat received five epidural cortical electrodes (stainlesssteel screw; diameter, 500 µm) and four depth electrodes (insulated nichrome wires; diameter, 200 µm). Reference was set at the cerebellum. The coordinates of the electrodes and the 9-electrode montage with their typical EEG during three states (WR, SWS, and REMS) were shown in Table 1 and Figure 1, respectively. The temporal electrode implantation followed the procedure introduced by Meeren et al. (2001). Two electromyogram (EMG) electrodes were implanted bilaterally in the dorsal neck muscles. After the surgical procedure, penicillin G was used for anti-infection. All rats were given at least 2 weeks to recover before the recording sessions started.

Recordings
Animals were habituated to the experimental environment 2 days prior to the recording. The connector on the rat skull was connected to a long recording lead, which allowed the subject to move freely in a homemade glass box (40 × 50 × 60 cm). Video recording was synchronized with the signal acquisition. Local field potential (LFP), electrocorticogram (ECoG) and EMG were recorded with a signal acquisition system (Chengyi, RM62160, China). Data were continuously recorded for 72 h (recording onset at 12:00). Bandpass filters were set between 0.16 and 100 Hz for EEG and between 8.3 and 500 Hz for EMG (50 Hz notch filter). The sample frequency was set at 1000 Hz. All recordings were stored on a hard disk (Lenovo Company, USA) for further analysis. All experiments were performed in a noise-attenuated room in which the environmental background noise was 32.2 ± 3.0 dB (mean ± SD). Other environmental variables (light: 12-h light/dark cycle [white lights on at 8:00]; food and water: ad libitum access; temperature: 25 ± 0.5 degree centigrade) were maintained. The experimenter entered the noise-attenuated room daily to clean cages and replace food and water at 12:00.

Histological Tests
To precisely determine the anatomical location of the depth electrodes, the animals were deeply anesthetized with chloral hydrate (300 mg/kg) after experiments. Then, the animals were perfused intracardially with saline followed by a 4% paraformaldehyde phosphate buffered solution. Brains were removed and stored in the paraformaldehyde phosphate buffer for 1 day and gradient dehydration with sucrose (15 and 30%) 2 days before sectioning. Serial coronal sections of 30 µm were cut on a freezing microtome, attached on Poly-L-Lysine-coated slides, and stained with ferricchloride solution (Riboni et al., 1991). After the histological inspection, the rats with any electrode out of the designed anatomical location were excluded (3 rats), and data from the other 10 rats were included for further analyses ( Figure S1).

Data Selection
To ensure that the rats were sufficiently adapted to the recording environment before obtaining data, only the data acquired on the third day (the last 24 h) were included in our analysis. Both SWS and REMS waves were selected from 8:00 to 12:00 (the last 4 of the 24 h, this period was rich in sleep activity). For each SWS and REMS segment, only state period >60 s, then the epoch was selected excluding the two sides of the data about 20 s (the epoch was selected as an integral multiple of 5 s). The 12:00 to 15:00 period was excluded because of the disturbance of the experimenter. If the selected data segments were less than the assumed 10 segments, the data from 15:00 to 20:00 were then included (from the 4th to the 8th h of the last 24 h). The wakeful rest data were chosen from the rest of the 17 h daily cycle, except from 8:00 to 15:00 (from the 4th to the 20th h of the last 24 h). For the wakeful rest data, only epochs during this state and more than 7 s were used to select the median 5 s segments for further analysis. The state selection was based on characteristics of EEG, MEG, and behaviors, which was systematically summarized in Table 2 (also see Figure 1).

Data Processing
Matlab (release 2013a) and SPSS software (release 19.0) were combined for the data processing. First, each segment of SWS and REMS was divided into 5-s epochs. Each epoch for the three states was band-pass filtered (0.5-34 Hz, to exclude potential highfrequency noise including muscle activities), down sampled (at 256 Hz) and de-trended (remove the linear trend). Power spectra with 1 Hz resolution were calculated by the Welch's method with the Hamming window. All the power spectra were first logtransformed and then averaged over brain areas, states and rats (one value for each brain area, state and rat) for further statistical analysis.
In this study, for each rat the log-transformed power spectra of each segment acquired from each electrode in SWS and REMS was averaged and the 60 epochs of WR were averaged in blocks of 6. Then 10 averaged values were computed for each state, electrode and rat. To compute frequency covariant with each other and to extract the frequency bands of each brain area across states, these 10 values were pooled for rats and brain areas for factor analysis (in each state, 200 values for PrL, CG, V2, and TE; 100 values for RSC). Each brain area under each state was submitted to factor analysis with EEG frequencies as variables and the PCA as the extraction method. The following criteria were used: the Kaiser-Guttman criterion (eigenvalues higher than 1 for eigenvectors) combined with the scree slope method and frequency bands with loading factor above 0.5 were included. The varimax method was used for factor rotation. To further check the independence of factors, the promax method was also used. Although most of the correlations between factors reached a meaningful significance level, which reflects that the factors were correlated with one another (Figure S2), the results from the two methods were almost identical. Therefore, only the promax-based results are presented.

Statistical Analysis
In order to evaluate differences in power spectra, 3-way repeated measures ANOVA was used for within-subject variables (state, brain region, and frequency). Both main effects and interactions were examined. Simple or simple-simple effects were further applied when the interactions were significant. Partial η 2 were employed for estimating of effect size of ANOVAs, and the values of 0.2, 0.5, and 0.8 were successively corresponding to small, medium, and large effect size (Cohen, 1992). For post-hoc multiple comparisons, the paired t-tests were used. The threshold of significance was set at a Bonferroni-adjusted p-value of 0.05.

Frequency Bands for Different States and Brain Regions
To guarantee the suitability of our data for factor analysis and to ensure the reliability of the results, both Kaiser-Meyer-Olkin (KMO) test and Bartlett's test were performed.  Table 4).
Combining with the Kaiser-Guttman criterion and the scree slope method (Figure 2), during WR, a fast frequency band from 21-22 to 32 Hz and an intermediate band from 8-9 to 19-21 Hz were detected for all regions; two slow bands from 1 to 4 Hz and from 4-5 to 6-8 Hz were extracted for PrL, CG, and RSC, whereas, only one slow band from 1 to 6 Hz was identified for V2 and TE (Table 4 and Figure 3).
For SWS, a consistent slow band from 1 to 6 Hz was detected for PrL, CG, RSC, and V2, but a slow band from 1 to 7 Hz  Figure 3). During REMS, two slow bands from 1 to 4-5 Hz and 5-6 to 7 Hz were detected for PrL, CG, RSC, and V2; a narrow band from 8 to 10-11 Hz that co-varied with higher frequencies from 13-15 to 16-18 Hz was also identified for PrL, CG, RSC, and V2; a fast band from 12-18 to 32 Hz was extracted for PrL, CG, and V2, whereas, two fast bands of 18-27 and 28-32 were identified for RSC, respectively; the division of frequency bands for TE was unique, including 1-15, 15-19, and 20-32 Hz ( Table 4 and Figure 3).
To further reveal the distribution of power spectra across brain areas, the positive peaks in each region were calculated. Coincident slow oscillation frequency peaks occurred at 2 Hz during WR and SWS, and coincident frequency peaks occurred at 7 and 14 Hz during REMS. However, the 6 Hz peaks were found among regions (except PrL) during WR, the 2 Hz peaks were identified among regions (except CG) during REMS, and the 23 Hz peaks were detected among regions (except TE) during REMS (Figure 5).
To further confirm the subdivision of spindle and theta oscillations, epochs selected from SWS and REMS were examined with time-frequency analysis. Original EEG activities and their spectrograms of RSC reveal two types of spindles with different waveforms and spectrograms (Figures 6A,B), and two types of theta activities with discrete spectrograms ( Figure 6C). Additionally, these diverse oscillation activities happened at separate moments (Figure 6), which further support the frequency band subdivisions of spindle and theta.

DISCUSSIONS
The present study used factor analysis to extract the frequency bands across brain regions during WR, SWS, and REMS. We found that frequency band divisions varied both across vigilance states and brain regions. These results suggested that the frequency bands divisions were determined by vigilance states and specific neural mass in distinct brain regions.

The Frequency Band Divisions of EEG
In both human and rat, two types of REMS can be differentiated, i.e., tonic and phasic REMS, which the latter is characterized by ponto-geniculo-occipital (PGO) waves (pontine waves in FIGURE 2 | Scree plots derived from the scree slope method for the five brain regions during WR, SWS, and REMS, respectively. rats), increased in the frequency of theta waves, and muscle twitches (Horne, 2000;Karashima et al., 2005;Montgomery et al., 2008). In the present work, obvious theta oscillations with 7 Hz peaks were found across all the brain regions during REMS. However, unlike the customary frequency band division in rats, the theta was further divided into two bands (5-7 and 8-11 Hz) in most regions. These divisions are in line with behavioral evidence of two kinds of REMS (Karashima et al., 2005;Montgomery et al., 2008). In humans, thalamo-cortical network specifically activates during phasic REMS, accompanying with lower alertness compared with tonic REMS (Wehrle et al., 2007). In rats, although both types of REMS are speculated to be involved in mnemonic process (Karashima et al., 2005), recent studies suggested the diverse roles of them in the mnemonic process (Montgomery et al., 2008;Brankack et al., 2012). Thus, the subdivision of theta bands may be more coincident with the brain functional states.
FIGURE 3 | Illustrations of the frequency bands and power spectra. Log transformed EEG power spectra obtained from the mean of the entire group (n = 10). Colored areas represent the frequency bands clustered by factor analysis, corresponding to the band divisions in Table 4. The saw-tooth patterns indicate the overlaps of frequency bands.
Previous studies have reported that there are two different sleep spindles, i.e., low-frequency and high-voltage spindles (HVSs: 6-10 Hz, primarily 7-8 Hz) and high-frequency and lowvoltage spindles (LVSs: 6-20 Hz, primarily 10-20 Hz) (Kandel and Buzsaki, 1997;Johnson et al., 2010), and the former has a larger amplitude (about 3 to 5 times) than the latter (Buzsaki, 1991). Consistent with these findings, we have identified similar two bands (6-10 and 9-21 Hz) with different power spectra during SWS using both factor analysis ( Table 4) and time-frequency analysis (Figure 6). Despite that HVSs and LVSs occupy similar or sometimes overlapping frequency bands and share many of the cellular mechanisms (Kandel and Buzsaki, 1997), their functional significances are quite different. LVSs are believed to be correlated with memory consolidation, whereas HVSs are not or even associated with memory deficits (Johnson et al., 2010). Compared with the classical division of the spindle, the exact subdivision of spindle bands may be more meaningful. FIGURE 4 | Mean (n = 10) and standard deviation (shaded areas) of EEG power spectra during WR, SWS, and REMS for each brain region. All the solid symbols indicate significant differences between states (p < 0.05, Bonferroni corrected).

Variation of Functional Band Divisions Across Brain Regions
During SWS, the frequency band divisions of LVSs were different along the anterior-posterior regions, i.e., the "slow" spindles (∼16 Hz) that were restricted to frontal regions (PrL and CG) and the "fast" spindles (∼21 Hz) that were found in parietal areas (RSC and V2). To the best of our knowledge, no studies related to the spindle activities have distinguished the two spindle components in rats. What's more, from the perspective of the spatial distribution, the current findings are consistent with the results in humans (9-12 and 12-15 Hz for slow and fast spindles, respectively) and cats (12-14 and 14-16 Hz for slow and fast spindles, respectively) (Mölle et al., 2011;Timofeev and Chauvette, 2013). Moreover, these two types of spindles may originate from different sources (Mölle et al., 2011;Timofeev and Chauvette, 2013), and take part in the diverse phases of memory consolidation (Mölle et al., 2011).
During WR, unlike the other regions, functional binding of delta and theta was found in both V2 and TE (Figure 3). In this state, rat hippocampus shows large irregular activity (Vanderwolf and Robinson, 1981), moreover, theta activity can originate in RSC independent from hippocampus (Talk et al., 2004;Young and McNaughton, 2009). It is reasonable to infer that the theta oscillation in WR might originate from RSC (note that this speculation doesn't exclude other potential generators). Furthermore, RSC tightly connected to the anterior thalamic nuclei (Vann et al., 2009;Dalrymple-Alford et al., 2015), and theta-like unit activity can be found in RSC and anterior thalamic nuclei simultaneously (Talk et al., 2004). Ultimately, the co-variation of delta and theta in V2 and TE might reflect a coordinated activities of the thalamus and RSC. Functionally, this co-variation might facilitate the information integration across sensory regions, supporting the past remembering and/or future imaging in this behavior state (Vann et al., 2009;Mednick et al., 2011).
For REMS, the phasic RMES theta co-varied with frequencies from 13-15 to 16-18 Hz across regions (except TE). Because phase synchronizations may occur for harmonic oscillation (Pletzer et al., 2010), and the 14 Hz peaks (the second harmonic of 7 Hz) were consistently found across regions. These co-variations of the phasic REMS theta bands with additional frequency bands might only indicate the effect of harmonic frequency relationships rather than the functional oscillation binding. There is no clear explanation for the co-variations in TE both FIGURE 5 | Mean (n = 10) of EEG power spectra in the five brain areas during WR, SWS, and REMS. Shaded areas denote coincident power peaks across brain areas, and arrows indicate the majority of peaks among regions. during SWS and REMS, but it cannot be attributed to probable artifacts, since both SWS and REMS have lower muscle activities compared with WR. As for PrL, the HVSs were not identified in this region. However, the same frequency division (7-14 Hz for spindles) was also found in the previous study (Siapas and Wilson, 1998), reflecting probably PrL has not been involved in HVS activity.
Compared with other brain regions, the theta division of RSC perfectly fit the power spectra during WR (Figure 3). Interestingly, the RSC may be the hypothetical generator (discussed above) of this theta activity in WR. Similarly, hippocampal theta activity in rodents can be recorded in other neocortex (Petsche and Stumpf, 1960;Winson, 1974;Bland and Whishaw, 1976). In the current study, band divisions of structures (RSC and V2) above hippocampus have better fit with power spectra than others during REMS. These findings indicate that the results of factor analysis may be affected by the constructions and locations related to the underlying generators to some extent, for which the difference of the division may be caused by the different signal-to-noise ratios. However, it should be noted that this speculation cannot exclude the influences of other potential generators. Actually, many independent theta generators have been observed in rats (Bilkey and Heinemann, 1999;Kahana et al., 2001;Seidenbecher et al., 2003;Young and McNaughton, 2009). Nevertheless, RSC and hippocampus may dominate the theta activities during WR and REMS, respectively. Therefore, the results of factor analysis might be influenced by the potential generators, and spatial factor should be considered in frequency band divisions when using factor analysis in rats.
The specific frequency bands identified by factor analysis across brain regions and states reflect the potential generators and functional significances. Therefore, it is worth noting that frequency band divisions based on distributed brain regions are meaningful for brain functional studies.

Comparison with a Previous Study
Based on our knowledge, there is one closely related study in rats (Corsi-Cabrera et al., 2001). The present study has three obvious In the spectrogram of the EEG (bottom), LVSs correspond to a peak range from 10 to 20 Hz. (C) EEG of tonic and phasic REM occur during periods of REMS (top). In the spectrogram of the EEG (bottom), tonic REM is marked by a strong peak around 6 Hz, however, phasic REM is marked around 10 Hz. Red lines indicate the onset of epochs. Spectrograms were calculated with EEGLAB (http://sccn.ucsd.edu/eeglab/; default parameter). Bars are normalized and hot colors reflect high power.
differences from the previous one. First of all, during WR, the delta activities were identified only in the present work. This is quite possibly caused by the difference between wakefulness and wakeful rest, and the delta activity is consistent with the findings in Young and McNaughton (2009). Secondly, during REMS, the theta oscillations were divided into two bands in this study. This difference may be due to the different data samples. In the work of Corsi-Cabrera et al., the data size is relatively small, and the phasic REMS occupy only approximately 5% of the total REMS (Montgomery et al., 2008). Finally, despite the similarity of statistical changes of the power spectra (Corsi-Cabrera et al., 2001), however, in the current study, we found a significant decrease of power (range from 10 to 28 Hz) during REMS compared with WR in PrL, CG, and RSC. These decreased power are supported by the brain imaging study in humans, in which a lack of increased blood flow in the homologous brain regions during REMS (Maquet et al., 1996). Moreover, these homologous brain regions in humans densely connected with the medial pulvinar nucleus (PuM) of thalamus, and decreased power (beta and gamma) was also found in PuM during REMS (Magnin et al., 2004). Interestingly, unlike in PrL, CG, and RSC, this decrease was not found in V2 and TE (also in results of Corsi-Cabrera et al.). It is possible that PrL, CG, and RSC are all association cortex with extensive connections (Vogt and Miller, 1983;Conde et al., 1990Conde et al., , 1995Vann et al., 2009;Vann, 2013). This situation could reflect an incomplete interaction between brain regions, a phenomenon resulting in a distracting state of REMS (Horne, 2000).
Although our method used in the current study is similar to the former one, we further considered the spatial factors on the frequency band divisions and used a bigger data size, with which the results observed here are also interesting. Especially, in the present work with a bigger data size, we found the subdivision of spindle and theta activities, which were not found in the former one. Furthermore, we also found functional relevant frequencies (delta and theta) co-vary in some brain regions and changed power of some frequency bands limited in other brain regions. We believe these findings could provide us some new insights into brain oscillations across brain regions and vigilance states.

Limitations
Although there are numerous advantages in our current study using factor analysis with considerations of spatial factors, we cannot exclude several limitations. Firstly, the results of factor analysis are based on power spectra, it should be noted that the power spectra are greatly affected by reference electrode. However, this limitation should not affect the conclusions since (1) the absolute mean power of cerebellar activity is several folds lower than that at the cerebral level, cerebellum as referenced is more like a rest reference at infinity than other cortexes (Culic et al., 2003) and is widely used in the rat electrophysiological study; (2) the power spectra in current work are similar to that in the former study, in which the reference electrode was not placed above the cerebellum (Corsi-Cabrera et al., 2001); (3) in this work, factor analysis was also conducted with average reference, and the results were similar (data not shown).
Secondly, present results were obtained with two recording technologies (ECoG and LFP), which allows for different spatial resolutions. However, this limitation should not affect the conclusions since (1) each brain region was used only one recording technology and independently conducted with factor analysis; (2) differences of band divisions widely distributed across regions, rather than limited between two different recording technologies.
Finally, it should be pointed out that EEG properties differ along the left-right axes of the rat brain (Fang et al., 2010). Whereas, lateralization of frequency band divisions was not taken into consideration in the current work, this question will be addressed in future work.

Summary
In this study, we used factor analysis on averaged EEG power spectra derived from multiple brain regions across vigilance states in rats. We found that the frequency band divisions changed both across states and brain regions. These findings not only provide more precise frequency band divisions, but also promote us to consider the impacts of spatial factors on the frequency band divisions.

AUTHOR CONTRIBUTIONS
DY, YX, and WJ: conceived and designed the study. WJ, YW, and MX: conducted the experiments. WJ, MC, DG, and GF: performed data analysis and prepared the draft manuscript. WJ, DY, and GF: reviewed data interpretation. DY and YX: edited and approved final manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fncom. 2016.00079