|Home | About | Journals | Submit | Contact Us | Français|
The resistive or non-resistive nature of the extracellular space in the brain is still debated, and is an important issue for correctly modeling extracellular potentials. Here, we first show theoretically that if the medium is resistive, the frequency scaling should be the same for electroencephalogram (EEG) and magnetoencephalogram (MEG) signals at low frequencies (<10 Hz). To test this prediction, we analyzed the spectrum of simultaneous EEG and MEG measurements in four human subjects. The frequency scaling of EEG displays coherent variations across the brain, in general between 1/f and 1/f2, and tends to be smaller in parietal/temporal regions. In a given region, although the variability of the frequency scaling exponent was higher for MEG compared to EEG, both signals consistently scale with a different exponent. In some cases, the scaling was similar, but only when the signal-to-noise ratio of the MEG was low. Several methods of noise correction for environmental and instrumental noise were tested, and they all increased the difference between EEG and MEG scaling. In conclusion, there is a significant difference in frequency scaling between EEG and MEG, which can be explained if the extracellular medium (including other layers such as dura matter and skull) is globally non-resistive.
Electronic supplementary material The online version of this article (doi:10.1007/s10827-010-0263-2) contains supplementary material, which is available to authorized users.
An issue central to modeling local field potentials is whether the extracellular space around neurons can be considered as a resistive medium. A resistive medium is equivalent to replacing the medium by a simple resistance, which considerably simplifies the computation of local field potentials, as the equations to calculate extracellular fields are very simple and based on Coulomb’s law (Rall and Shepherd 1968; Nunez and Srinivasan 2005). Forward models of the EEG and inverse solution/source localization methods also assume that the medium is resistive (Sarvas 1987; Wolters and de Munck 2007; Ramirez 2008). However, if the medium is non-resistive, the equations governing the extracellular potential can be considerably more complex because the quasi-static approximation of Maxwell equations cannot be made (Bédard et al. 2004).
Experimental characterizations of extracellular resistivity are contradictory. Some experiments reported that the conductivity is strongly frequency dependent, and thus that the medium is non-resistive (Ranck 1963; Gabriel et al. 1996a, b, c). Other experiments reported that the medium was essentially resistive (Logothetis et al. 2007). However, both types of measurements used current intensities far larger than physiological currents, which can mask the filtering properties of the tissue by preventing phenomena such as ionic diffusion (Bédard and Destexhe 2009). Unfortunately, the issue is still open because there exists no measurements to date using (weak) current intensities that would be more compatible with biological current sources.
In the present paper, we propose an indirect method to estimate if extracellular space can be considered as a purely resistive medium. We start from Maxwell equations and show that if the medium was resistive, the frequency-scaling of electroencephalogram (EEG) and magnetoencephalogram (MEG) recordings should be the same. We then test this scaling on simultaneous EEG and MEG measurements in humans.
We recorded the electromagnetic field of the brain during quiet wakefulness (with alpha rhythm occasionally present) from four healthy adults (4 males ages 20–35). Participants had no neurological problems including sleep disorders, epilepsy, or substance dependence, were taking no medications and did not consume caffeine or alcohol on the day of the recording. We used a whole-head MEG scanner (Neuromag Elekta) within a magnetically shielded room (IMEDCO, Hagendorf, Switzerland) and recorded simultaneously with 60 channels of EEG and 306 MEG channels (Nenonen et al. 2004). MEG SQUID (super conducting quantum interference device) sensors are arranged as triplets at 102 locations; each location contains one “magnetometer” and two orthogonal planar “gradiometers” (GRAD1, GRAD2). Unless otherwise noted, MEG will be used here to refer to the magnetometer recordings. Locations of the EEG electrodes on the scalp of individual subjects were recorded using a 3D digitizer (Polhemus FastTrack). HPI (head position index) coils were used to measure the spatial relationship between the head and scanner. Electrode arrangements were constructed from the projection of 3D position of electrodes to a 2D plane in order to map the frequency scaling exponent in a topographical manner. All EEG recordings were monopolar with a common reference. Sampling rate was 1,000 Hz.
For all subjects, four types of consecutive recordings were obtained, in the following order: (1) Empty-room recording; (2) Awake “idle” recording where subjects were asked to stay comfortable, without movements in the scanner, and not to focus on anything specific; (3) a visual task; (4) sleep recordings. All idle recordings used here were made in awake subjects with eyes open, where the EEG was desynchronized. A few minutes of such idle time was recorded in the scanner. For each subject, three awake segments with duration of 60 s were selected from the idle recordings (see example signals in Fig. 1).
As electrocardiogram (ECG) noise often contaminates MEG recordings, Independent component analysis (ICA) algorithm was used to remove such contamination; either Infomax (Bell and Sejnowski 1995) or the “Jade algorithm” from the EEGLAB toolbox (Delorme and Makeig 2004) was used to achieve proper decontamination. In all recordings, the ECG component stood out very robustly. In order not to impose any change in the frequency content of the signal, we did not use the ICA to filter the data on any prominent independent oscillatory component and it was solely used to decontaminate the ECG noise. We verified that the removal of ECG did not change the scaling exponent (not shown).
In each recording session, just prior to brain recordings, we recorded a few minutes of the electromagnetic field present within the dewar in the magnetic shielded room. Similar to wake epochs, three segments of 60 s. duration were selected for each of the four recordings. This will be referred to “empty room” recordings and will be used in noise correction of the awake recordings.
In each subject, the power spectral density (PSD) was calculated by first computing the Fast Fourier transform (FFT) of three awake epochs, then averaging their respective PSDs (square modulus of the FFT). This averaged PSD was computed for all EEG and MEG channels in order to reduce the effects of spurious peaks due to random fluctuations. The same procedure was also followed for empty-room signals.
Because the environmental and instrumental sources of noise are potentially high in MEG recordings, we took advantage of the availability of empty-room recordings to correct for the presence of noise in the signal. We used five different methods for noise correction, based on different assumptions about the nature of the noise. We describe below these different correction methods, while all the details are given in Supplementary Methods.
A first procedure for noise correction, exponent subtraction (ES), assumes that the noise is intrinsic to the SQUID sensors. This is justified by the fact that the frequency scaling of some of the channels is identical to that of the corresponding empty-room recording (see Sections 4.2 and 4.3). In such a case, the scaling is assumed to entirely result from the “filtering” of the sensor, and thus the correction amounts to subtract the scaling exponents.
A second class of noise subtraction methods assume that the noise is of ambient nature and is uncorrelated with the signal. This characteristics, warrants the use of spectral subtraction (where one subtracts the PSD of the empty-room from that of the MEG recordings), prior to the calculation of the scaling exponent. The simplest form of spectral subtraction, linear multiband spectral subtraction (LMSS), treats the sensors individually and does not use any spatial/frequency-based statistics in its methodology (Boll 1979). An improved version, nonlinear multiband spectral subtraction (NMSS), takes into account the signal-to-noise ratio (SNR) and its spatial and frequency characteristics (Kamath and Loizou 2002; Loizou 2007). A third type, Wiener filtering (WF), uses a similar approach as the latter, but obtain an estimate of the noiseless signal from that of the noisy measurement through minimizing the Mean Square Error (MSE) between the desired and the measured signal (Lim and Oppenheim 1979; Abd El-Fattah et al. 2008).
A third type of noise subtraction, partial least squares (PLS) regression, combines Principal component analysis (PCA) methods with multiple linear regression (Abdi 2010; Garthwaite 1994). This methods finds the spectral patterns that are common in the MEG and the empty-room noise, and removes these patterns from the PSD.
The method to estimate the frequency scaling exponent was composed of steps: First, applying a spline to obtain a smooth FFT without losing the resolution (as can happen by using other spectral estimation methods); Second, using a simple polynomial fit to obtain the scaling exponent. To improve the slope estimation, we approximated the PSD data points using a spline, which is a series of piecewise polynomials with smooth transitions and where the break points (“knots”) are specified. We used the so-called “B-spline” (see details in de Boor 2001).
The knots were first defined as linearly related to logarithm of the frequency, which naturally gives more resolution to low frequencies, to which our theory applies. Next, in each frequency window (between consecutive knots), we find the closest PSD value to the mean PSD of that window. Then we use the corresponding frequency as the optimized knot in that frequency range, leading the final values of the knots. The resulting knots stay close to the initial distribution of frequency knots but are modified based on each sensor’s PSD data to provide the optimal knot points for that given sensor (Fig. 2(a)). We also use additional knots at the outer edges of the signal to avoid boundary effects (Eilers and Marx 1996). The applied method provides a reliable and automated approach that uses our enforced initial frequency segments with a high emphasis in low frequency and it optimizes itself based on the data. After obtaining a smooth B-spline curve, a simple 1st degree polynomial fit was used to estimate the slope of the curve between 0.1–10 Hz (the fit was limited to this frequency band in order to avoid the possible effects of the visible peak at 10 Hz on the estimated exponent). Using this method provides a reliable and robust estimate of the slope of the PSD in logarithmic scale, as shown in Fig. 2(b). For more details on the issue of automatic non-parametric fitting, and the rationale behind combining the polynomial with spline basis functions, we refer the reader to Magee (1998) as well as Royston and Altman (1994) and Katkovnik et al. (2006).
This procedure was realized on all channels automatically (102 channels for MEG, 60 channels for EEG, for each patient). Every single fit was further visually confirmed. In the case of MEG, noise correction is essential to validate the results. For doing so, we used different methods (as described above) to reduce the noise. Next, all the mentioned steps of frequency scaling exponents were carried out on the corrected PSD. Results are shown in Fig. 4.
Three ROIs were selected for statistical comparisons of the topographic plots. As shown in Fig. 4(f), FR (Frontal) ROI refers to the frontal ellipsoid, VX (Vertex) ROI refers to the central disk located on vertex and PT (Parietotemporal) refers to the horseshoe ROI.
We start from first principles (Maxwell equations) and derive equations to describe EEG and MEG signals. Note that the formalism we present here is different than the one usually given (as in Plonsey 1969; Gulrajani 1998), because the linking equations are here considered in their most general expression (convolution integrals), in the case of a linear medium (see Eq. 77.4 in Landau and Lifchitz 1984). This generality is essential for the problem we treat here, because our aim is to compare EEG and MEG signals with the predictions from the theory, and thus the theory must be as general as possible.
Maxwell equations can be written as
If we suppose that the brain is linear in the electromagnetic sense (which is most likely), then we have the two following linking equations. The first equation links the electric displacement with the electric field:
where ε is a symmetric second-order tensor.
A second equation links magnetic induction and the magnetic field:
where μ is a symmetric second-order tensor.
where σ is a symmetric second-order tensor.2 Because the effect of electric induction (Faraday’s law) is negligible, we can write:
This system is much simpler compared to above, because electric field and magnetic induction are decoupled.
where ω=2πf and
where the relation in Eq. (7) is the current density produced by the (primary) current sources in the extracellular medium. Note that in this formulation, the electromagnetic parameters εf, μf and σf depend on frequency.3 This generalization is essential if we want the formalism to be valid for media that are linear but non-resistive, which can expressed with frequency-dependent electric parameters. It is also consistent with the Kramers–Kronig relations (see (Landau and Lifchitz 1984; Foster and Schwan 1989).
is the current density of these sources in Fourier frequency space. This current density is composed of the axial current in dendrites and axons, as well as the transmembrane current. Of course, this expression is such that at any given point, there is only one of these two terms which is non-zero. This is a way of preserving the linearity of Maxwell equations. Such a procedure is legitimate because the sources are not affected by the field they produce.4
From Eq. (6) (Faraday’s law in Fourier space), we can write:
From Eq. (6) (Ampère–Maxwell’s law in Fourier space), we can write:
Setting γf=σf+iωεf, one obtains:
where is a source term and γf is a symmetric second-order tensor (3×3). Note that this tensor depends on position and frequency in general, and cannot be factorized. We will call this expression (Eq. (10)) the “first fundamental equation” of the problem.
From the mathematical identity
it is clear that this is sufficient to know the divergence and the curl of a field , because the solution of is unique with adequate boundary conditions.
As in the case of magnetic induction, the divergence is necessarily zero, it is sufficient to give an explicit expression of the curl as a function of the sources.
We consider the following boundary conditions:
The latter expression can be proven by calculating the total current and apply the divergence theorem (not shown).
The “second fundamental equation” above implies inverting γf, which is not possible in general, because it would require prior knowledge of both conductivity and permittivity in each point outside of the sources. If the medium is purely resistive (γf=γ where γ is independent of space and frequency), one can evaluate the electric field first, and next integrate using the quasi-static approximation (Ampère–Maxwell’s law). Because for low frequencies, we have necessarily , we obtain
which is also known as Ampère’s law in Fourier space.
Thus, for low frequencies, one can skip the second fundamental equation. Note that in case this quasi- static approximation cannot be made (such as for high frequencies), then one needs to solve the full system using both fundamental equations. Such high frequencies are, however, well beyond the physiological range, so for EEG and MEG signals, the quasi-static approximation holds if the extracellular medium is resistive, or more generally if the medium satisfies (see Eqs. (5) and (6)).
According to the quasi-static approximation, and using the linking equation between current density and (Eq. (7)), we can write:
Because the divergence of magnetic induction is zero, we have from Eq. (11):
This equation can be easily integrated using Poisson integral (“Poisson equation” for each component in Cartesian coordinates) In Fourier space, this integral is given by the following expression
If the medium is purely resistive (“ohmic”), then γ does not depend on the spatial position (see Bédard et al. 2004; Bédard and Destexhe 2009) nor on frequency, so that the solution for the magnetic induction is given by:
and does not depend on the nature of the medium.
For the electric potential, from Eq. (10), we obtain the solution:
Thus, when the two source terms and are white noise, the magnetic induction and electric field must have the same frequency dependence. Moreover, because the spatial dimensions of the sources are very small (see Appendices), we can suppose that the current density is given by a function of the form:
such that and have the same frequency dependence for low frequencies. Equation (19) constitutes the main assumption of this formalism.
In Appendix A, we provide a more detailed justification of this assumption, based on the differential expressions of the electric field and magnetic induction in a dendritic cable. Note that this assumption is most likely valid for states with low correlation such as desynchronized-EEG states or high-conductance states, and for low-frequencies, as we analyze here (see details in the Appendices).
Thus, the main prediction of this formalism is that if the extracellular medium is resistive, then the PSD of the magnetic induction and of the electric potential must have the same frequency dependence. In the next section, we will examine if this is the case for simultaneously recorded MEG and EEG signals.
A total of four subjects were used for the analysis. Figure 1 shows sample MEG and EEG channels from one of the subjects, during quiet wakefulness. Although the subjects had eyes open, a low-amplitude alpha rhythm was occasionally present (as visible in Fig. 1). There were also oscillations present in the empty-room signal, but these oscillations are evidently different from the alpha rhythm because of their low amplitude and the fact that they do not appear in gradiometers (see Supplementary Fig. S1).
In the next sections, we start by briefly presenting the method that was used to estimate the frequency scaling of the PSDs. Then we report the scaling exponents for 0.1–10 Hz frequency bands and their differences in EEG and MEG recordings.
Because of the large number of signals in the EEG and MEG recordings, we used an automatic non-parametric procedure to estimate the frequency scaling (see Section 2). We used a B-spline approximation by interpolation with boundary conditions to find a curve which best represents the data (see Section 2). A high density of knots was given to the low-frequency band (0.1–10 Hz), to have an accurate representation of the PSD in this band, and calculate the frequency scaling. An example of optimized knots to an individual sensor is shown in Fig. 2(a); note that this distribution of knots is specific to this particular sensor. The resulting B-spline curves were used to estimate the frequency scaling exponent using a 1st degree polynomial fit. Figure 2(b) shows the result of the B-spline analysis with optimized knots (in green) capturing the essence of the data better than the usual approximation of the slope using polynomials (in red). The goodness of fit showed a robust estimation of the slope using B-spline method. Residuals were −0.01 ± 0.6 for empty-room, 0.2 ± 0.65 for MEG awake, 0.05 ± 0.6 for LMSS, 0.005 ± 0.64 for NMSS, 0.08 ± 0.5 for WF,0.001 ± 0.02 for PLS, and −0.02 ± 0.28 for EEG B-spline (all numbers to be multiplied by 10−14).
Figure 3 shows the results of the B-spline curve fits to the log-log PSD vs frequency for all sensors of all subjects. In this figure, and only for the ease of visual comparison, these curves were normalized to the value of the log(PSD) of the highest frequency. As can be appreciated, all MEG sensors (in red) show a different slope than that of the EEG sensors (in blue). The frequency scaling exponent of the EEG is close to 1 (1/f scaling), while MEG seems to scale differently. Thus, this representation already shows clear differences of scaling between EEG and MEG signals.
However, MEG signals may be affected by ambient or instrumental noise. To check for this, we have analyzed the empty-room signals using the same representation and techniques as for MEG, amd the results are represented in Fig. 3 (insets). Empty-room recordings always scale very closely to the MEG signal, and thus the scaling observed in MEG may be due in part to environmental noise or noise intrinsic to the detectors. This emphasizes that it is essential to use empty-room recordings made during the same experiment to correct the frequency scaling exponent of MEG recordings.
To correct for this bias, we have used five different procedures (see Section 2). The first class of procedure (ES) considers that the scaling of the MEG is entirely due to filtering by the sensors, which would explain the similar scaling between MEG and empty-room recordings. In this case, however, nearly all the scaling would be abolished, and the corrected MEG signal would be similar to white noise (scaling exponent close to zero). Because the similar scaling may be coincidental, we have used two other classes of noise correction procedures to comply with different assumptions about the nature of the noise. The second class, is composed of spectral subtraction (LMSS and NMSS) or Wiener filtering (see Section 2). These methods are well-established in other fields such as acoustics. The third class, uses statistical patterns of noise to enhance PSD (PLS method, for details see Section 2).
We applied the above methods to all channels and represented the scaling exponents in topographic plots in Fig. 4. This figure portrays that both MEG and EEG do not show a homogenous pattern of the scaling exponent, confirming the differences of scaling seen in Fig. 3. The EEG (Fig. 4(a)) shows that areas in the midline have values closer to one, while those at the margin can deviate from 1/f scaling. MEG on the other hand shows higher values of the exponent in the frontal area and a horseshoe pattern of low value exponents in parietotemporal regions (Fig. 4(b)). As anticipated above, empty-room recordings scale more or less uniformly with values close to 1/f (Fig. 4(c)), thus necessitating the correction for this phenomena to estimate the correct MEG frequency scaling exponent. Different methods for noise reduction are shown in Fig. 4: spectral subtraction methods, such as LMSS (Fig. 4(d)), NMSS (Fig. 4(e)), WF enhancement (Fig. 4(f)). These corrections preserve the pattern seen in Fig. 4(b), but tend to increase the difference with EEG scaling: one method (LMSS) yields minimal correction while the other two (NMSS and WF) use band-specific SNR information in order to cancel the effects of background colored-noise (see Supplementary Fig. S2), and achieve higher degree of correction (see Supplementary Methods for details). Figure 4(g) portrays the use of PLS to obtain a noiseless signal based on the noise measurements. The degree of correction achieved by this method is higher than what is achieved by spectral subtraction and WF methods. Exponent subtraction is shown in Fig. 4(h). This correction supposes that the scaling is due to the frequency response of the sensors, and nearly abolishes all the frequency scaling (see also Supplementary Fig. S3 for a comparison of different methods of noise subtraction).
Figure 5(a) represents the overall pattern providing evidence on the general difference and the wider variability in MEG recordings. The next three panels relate to the individual ROIs. Of the spectral subtraction methods, NMSS achieves a higher degree of correction in comparison with LMSS (see Fig 4(c), (d) as well as Supplementary Fig. S3). Because NMSS takes into account the effects of the background colored-noise (Supplementary Fig. S2), it is certainly more relevant to the type of signals analyzed here. The results of NMSS and WF are almost identical and confirm one another (see Fig. 4(e), as well as Supplementary Fig. S3). Therefore, of this family of noise correction, only NMSS is portrayed here. Of the methods dealing with different assumptions about the nature of the noise, the “Exponent subtraction” almost abolishes the frequency scaling (Also see in Fig. 4(h), as well as Supplementary Fig. S3). Applying PLS yields values in between “Exponent subtraction” and that of NMSS and is portrayed in Fig. 5.
In the Frontal region (Fig. 5(b)), the EEG scaling exponents show higher variance by comparison to MEG. Also, EEG shows some overlaps with the distribution curve of non-corrected MEG; this overlap becomes limited to the tail end of the NMSS correction and is abolished in the case of PLS correction. As can be appreciated, VX (Fig. 5(c)) shows both similar values and similar distribution for EEG and non-corrected MEG. These similarities, in terms of regional overall values and distribution curve, are further enhanced after NMSS correction. It is to be noted that, in contrast to these similarities, the one-to-one correlation of NMSS and EEG at VX ROI are very low (see below, Table 1(B), (C)). The values of PLS noise correction are very different from that of EEG and have a similar, but narrower, distribution curve shape. Two other ROIs show distinctively different values and distribution in comparing EEG and MEG. Both NMSS and PLS agree on this with PLS showing more extreme cases. Figure 5(d) reveals a bimodal distribution of MEG exponents in the parietotemporal region (PT ROI). This region has also the highest variance (in MEG scaling exponents) compared to other ROIS. The distinction between EEG and MEG is enhanced in PLS estimates; however, the variance of PT is reduced in comparison to NMSS while the bimodality is still preserved but weakened. The values of mean and standard deviation for these ROIs’ exponents are provided in Table 1(A) (mean ± standard deviation).
The box-plots of Fig. 5-plots further show the difference between the medians, lower/upper quartile and interquartile range. The overall difference is that the uncorrected MEG has much wider variance compared to EEG and corrected MEG (in case of PLS correction); the absolute value of the median of MEG (uncorrected, or corrected with either NMSS or PLS) is always smaller than that of EEG. The VX region is an exception to the above rules; interestingly, the one-to-one correlation of VX happens to be the lowest of all (see below). In the case of NMSS-corrected MEG, the shape of the pdf is preserved. However, PLS narrows the distribution curve of MEG but further enhances the differences between MEG and EEG. Therefore, median and lower/upper quartiles will have different value than that of EEG.
Correlation values (Table 1(B) (C)) show that, although VX ROI has the closest similarity in terms of its central tendency and probability distribution, it provides the lowest correlation in a pairwise fashion. P-values (for testing the hypothesis of no correlation against the alternative that there is a nonzero correlation) for Pearson’s correlation were calculated using a Student’s t-distribution for a transformation of the correlation and they were all significant (less than 10−15 for α = 0.05). Similarly, a non-parametric statistic Kendall tau rank correlation was used to measure the degree of correspondence between two rankings and assessing the significance of this correspondence between MEG and EEG in the selected ROIs (Table 1(C)). P-values for Kendall’s tau and Spearman’s rho calculate using the exact permutation distributions were all significant (less than 10−15 for α = 0.05). Kendall tau shows that the rank correlation for all areas considered together as well as for PT, show a lesser correlation than that is shown by Pearson linear correlation. Furthermore, we carried out a Kruskal-Wallis nonparametric version of one-way analysis of variance. We used this test to avoid bias in ANOVA (KruskalWallis assumes that the measurements come from a continuous distribution, but not necessarily a normal distribution as is assumed in ANOVA). KruskalWallis uses analysis of variance on the ranks of the data values, not the data values themselves and therefore is an appropriate test for comparison of the homogeneity of pattern between ROIs of two image as well as their statistical median. As shown in Table 1(D), all p-values were significant emphasizing the difference between the spatial aspect of the spectral nature of MEG and EEG. Note that the difference of scaling exponent of EEG and MEG was also confirmed by nonlinear spatial kendall correlation analysis, independently of the ROIs classification (not shown).
Noise correction does not affect all the sensors in a same fashion. As presented in Supplementary Fig. S3, the simple linear spectral subtraction (LMSS) may lead to an increment or decrement of the scaling exponent. In any case, the correction achieved by this method is minimal. This is due to the fact that LMSS ignores the complex non-linear patterns of the SNR in different channels (Supplementary Fig. S2). We show that for all subjects, as the frequency goes up, the SNR goes down. It is also noticable that in each defined frequency band, i.e. 0–10 Hz (Slow, Delta and Theta), 11–30 Hz (Beta), 30–80 Hz (Gamma), 80–200 Hz (Fast oscillation), 200–500 Hz (Ultra-fast oscillation), there is an observable sensor-to-sensor SNR variability. This variability is at its maximum in the band with the highest SNR (i.e. 1–10 Hz). All together, the non-linear nature of MEG SNR shows that a linear spectral subtraction could behave non-optimally, leading to minimal correction. This also conveys that the optimal spectral correction can be achieved only by non-linear methods that explicitly take into account the SNR information of the data. Therefore the correction achieved by NMSS and WF have higher validity, in agreement with the fact that both methods yield similar results in terms of values and spatial distribution (Fig. 4(e), (f)).
In this paper, we have used a combination of theoretical and experimental analyses to investigate the spectral structure of EEG and MEG signals. In the first part of the paper, we presented a theoretical investigation showing that if the extracellular medium is purely resistive, the equations of the frequency dependence of electric field and magnetic induction take a simple form, because the admittance tensor does not depend on spatial coordinates. Thus, the macroscopic magnetic induction does not depend on the electric field outside the neuronal sources, but only depends on currents inside neurons. In this case, the frequency scaling of the PSD should be the same for EEG and MEG signals. This conclusion is only valid in the linear regime, and for low frequencies.
An assumption behind this formalism is that the spatial and frequency dependence of the current density factorize (Eq. (19)). We have shown in the Appendices that this is equivalent to consider the different current sources as independent. Thus, the formalism will best apply to states where the activity of synapses is intense and of very low correlation. This is the case for desynchronized-EEG states or more generally “high-conductance states”, in which the activity of neurons is intense, of low correlation, and the neuronal membrane is dominated by synaptic conductances (Destexhe et al. 2003). In such conditions, the dendrites are bombarded by intense synaptic inputs which are essentially uncorrelated, and one can consider the current sources as independent (Bédard et al. 2010). In the present paper, we analyzed EEG and MEG recordings in such desynchronized states, where this formalism best applies.
Note that the above reasoning neglects the possible effect of abrupt variations of impedances between different media (e.g., between dura matter and cerebrospinal fluid). However, there is evidence that this may not be influential. First, our previous modeling work (Bédard et al. 2004) showed that abrupt variations of impedance have a negligible effect on low frequencies, suggesting that even in the presence of such abrupt variations should not play a role at low frequencies. Second, in the frequency range considered here, the skull and the skin are very close to be resistive at low frequencies (Gabriel et al. 1996b), so it is very unlikely that they play a role in the frequency scaling in EEG and MEG power spectra even at high frequencies.
In the second part of the paper, we have analyzed simultaneous EEG and MEG signals recorded in four healthy human subjects while awake and eyes open (with desynchronized EEG). Because of the large number of channels involved, we used an automatic procedure (B-splines analysis) to calculate the frequency scaling. As found in previous studies (Pritchard 1992; Freeman et al. 2000; Bédard et al. 2006a), we confirm here that the EEG displays frequency scaling close to 1/f at low frequencies.5 However, this 1/f scaling was most typical of the midline channels, while temporal and frontal leads tended to scale with slightly larger exponents, up to 1/f2 (see Fig. 4(a)). The same pattern was observed in all four patients.
This approach differs from previous studies in two aspects. First, in contrast to prior studies (such as Novikov et al. 1997; Linkenkaer-Hansen et al. 2001), we calculated the frequency scaling of all the sensors and did not confine our analysis to a specific region. Second, unlike other investigators (such as Hwa and Ferree 2002a, b), we did not limit our evaluations to either EEG or MEG alone, but rather analyzed the scaling of both type of signals simultaneously. Such a strategy enables us to provide an extended spatial analysis of the frequency scaling. It also provides a chance to compare the scaling properties of these signals in relation to their physical differences.
For the MEG recordings, the frequency scaling at low frequencies was significantly lower compared to the EEG (see Fig. 3). This difference in frequency scaling was also accompanied by spatial variability patterns (see Fig. 4) showing three distinct regions:(1) a frontal area where the exponents had their highest values in the case of MEG; (2) a central area where the values of exponents of EEG and MEG get closer to each other and (3) a parietotemporal horseshoe region showing the lowest exponents for MEG with bimodal characteristics (Fig. 5). In some cases, the scaling of the uncorrected and corrected MEG signal was also close to 1/f, as reported previously (Novikov et al. 1997). In the frontal area (FR mask), the scaling exponent of the EEG was generally larger. At Vertex (VX mask), EEG and MEG had similar values and at the Parietotemporal region (PT mask), MEG showed a bimodal property with a much broader range of scaling exponent in comparison to EEG (see Fig. 4). Note that to avoid the effect of spurious peaks, Novikov et al. used the spectrum of signal differences and argued for the existence of a local similarity regime in brain activity. This approach fundamentally changes the spectral characteristics of Magnetometers (which measure the absolute magnitude of the magnetic induction) into a measure that only for the neighboring sensors approximates the behavior of the gradiometers (which measures the gradient of the magnetic induction). So it is not clear how to relate their values to the ones obtained here.
To make sure that the differences of frequency scaling between EEG and MEG were not due to environmental or instrumental noise, we have used five different methods to remove the effect of noise. These methods are based on different assumptions about the nature and effect of the noise. A first possibility is to correct for the noise induced by the MEG sensors. It is known that the SQUID detectors used in MEG recordings are very sensitive to environmental noise and they can produce 1/f noise (Hämäläinen et al. 1993). Under this assumption, part of the scaling of the MEG could be due to “filtering” by the sensors themselves, which justifies a simple subtraction of scaling exponents to remove the effects of this filtering. Note that such empty-room recordings were not possible for the EEG, although the noise from the recording setup could be estimated (see Miller et al. 2009 for example). Because in some cases both MEG and emptyroom signals have similar scaling, a simple correction by subtracting the exponents would almost entirely abolish the frequency scaling while in other cases it may even revert the sign of the scaling exponent (see Fig. 4(h), Supplementary Fig. S3).
However, if noise is not due to the sensors but is of additive uncorrelated nature, then another method for noise correction must be used. For this reason, we have used a second class of well-established methods consisting of spectral subtraction (Boll 1979; Sim et al. 1998). Using three of such methods (LMSS, NMSS and WF) changed the scaling exponent, without fundamentally changing its spatial pattern (Fig. 4(d)–(f)). The largest correction was obtained by non-linear methods which take into account the SNR information in the MEG signal. We also applied another class of method which uses the collective characteristics of all frequencies in noise correction (PLS). Similar to exponent subtraction, this method nearly abolished all the scaling of the MEG (Fig. 4(g)). In conclusion, although different methods for noise subtraction give rise to different predictions about frequency scaling, all of the used methods enhanced the difference between EEG and MEG scaling. Thus, we conclude that the difference of EEG and MEG scaling cannot be attributed to noise, but is significant, therefore reinforcing the conclusion that the medium must be non-resistive.
An alternative method to investigate this is the “Detrended Fluctuation Analysis” (DFA; see Kantelhardt et al. 2001; Linkenkaer-Hansen et al. 2001; Hwa and Ferree 2002a, b). Like many nonlinear approaches, DFA results are very vulnerable to the selection of certain parameters. Different filters severely affect the scaling properties of the electromagnetic signals to different extents, and therefore the parameters estimated through the DFA analysis could be false or lead to distorted interpretations of real phenomena (Valencia et al. 2008), and these effects are especially prominent for lower frequencies, which are precisely our focus of investigation here. One of the fields for which DFA can provide robust results is to analyze surrogate data with known characteristics. Although the use of DFA to evaluate the scaling exponents of EEG was vigorously criticized (Valencia et al. 2008), a previous analysis (Hwa and Ferree 2002a, b) reported two different regions, a central and a more frontal, which somehow correlate with the FR and VX regions identified in our analysis. Similarly, a study by Buiatti et al. (Buiatti et al. (2007)) using DFA provided evidence for topographical differences in scaling exponents of EEG recordings. They report that scaling exponents were homogeneous over the posterior half of the scalp and became more pronounced toward the frontal areas. In contrast to Linkenkaer-Hansen et al. (2001) (where envelope of alpha oscillations was used for DFA estimation), this study uses the raw signal in its DFA analysis and yields values closer to those reported here.
Both uncorrected signals and empty-room correction show that there is a fundamentally different frequency scaling between EEG and MEG signals, with near-1/f scaling in EEG, while MEG shows a wider range at low frequencies. Although it is possible that non-neuronal sources affect the lower end (<1Hz) of the evaluated frequency domain (Voipio et al. 1989), the solution to avoid these possible effects remain limited to invasive methods such as inserting the electrode into the scalp (Ferree et al. 2001) or using intracranial EEG recordings (similar to Miller et al. 2009). This approach would render wide range spatial recording as well as simultaneous invasive EEG and MEG recordings technically demanding or impractical. However, if technically feasible, such methods could provide a way to bypass non-neuronal effects at very low frequency. It could also provide a chance to evaluate the effects of spatial correlation on spectral structure at a multiscale level.
The power spectral structure we observe here is consistent with a scenario proposed previously (Bédard et al. 2006a): the 1/f structure of the EEG and LFP signals is essentially due to a frequency-filtering effect of the signal through extracellular space; this type of scaling can be explained by ionic diffusion and its associated Warburg impedance6 (see Bédard and Destexhe 2009). It is also consistent with the matching of LFPs with multi-unit extracellular recordings, which can be reconciled only assuming a 1/f filter (Bédard et al. 2006a). Finally, it is also consistent with the recent evidence from the transfer function calculated from intracellular and LFP recordings, which also showed that the extracellular medium is well described by a Warburg impedance (Bédard et al. 2010, , submitted to this issue). If this non-resistive aspect of extracellular media is confirmed, it may influence the results of models of source localization, which may need to be reformulated by including more realistic extracellular impedances.
In conclusion, the present theoretical and experimental analysis suggests that the scaling of EEG and MEG signals cannot be reconciled using a resistive extracellular medium. The 1/f structure of EEG with smaller scaling exponents for MEG is consistent with non-resistive extracellular impedances, such as capacitive media or diffusion (Warburg) impedances. Including such impedances in the formalism is non trivial because these impedances are strongly frequency dependent. The Poisson integral (the solution of Poisson’s law ) would not apply anymore (see Bédard et al. 2004; Bédard and Destexhe 2009). Work is under way to generalize the formalism and include frequency-dependent impedances.
Finally, it is arguable that the scaling could also be influenced by the cancellation and the extent of spatial averaging of microscopic signals, which are different in EEG and MEG (for more details on cancellation see Ahlfors et al. 2010; for details on spatial sensitivity profile see Cuffin and Cohen 1979). Such a possible role of the complex geometrical arrangement of underlying current sources should be investigated by 3D models which could test specific assumptions about the geometry of the current sources and dipoles, and their possible effect on frequency scaling. Such a scenario constitutes another possible extension of the present study.
We thank Philip Louizo for comments on spectral subtraction methods and Hervé Abdi for comments on Partial least square methods. Research supported by the Centre National de la Recherche Scientifique (CNRS, France), Agence Nationale de la Recherche (ANR, France), the Future and Emerging Technologies program (FET, European Union; FACETS project) and the National Institutes of Health (NIH grants NS18741, EB009282 and NS44623). N.D. is supported by a fellowship from Ecole de Neurosciences de Paris (ENP). Additional information is available at http://cns.iaf.cnrs-gif.fr.
Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.
To compare the frequency dependence of magnetic induction and electric field, we evaluate them in a dendritic cable, expressed differentially. For a differential element of dendrite, in Fourier space, the current produced by a magnetic field (Ampère–Laplace law) is given by the following expression (see Appendix B):
when the extracellular medium is resistive. Note that the source of magnetic induction is essentially given by the component of along the axial direction () within each differential element of dendrite because the perpendicular (membrane) current does not participate to producing the magnetic induction if we assume a cylindrical symmetry.
For the electric potential, we have the following differential expression for a resistive medium (see Appendix C):
where is the transmembrane current per unit of surface.
If we consider the differential expressions for the magnetic induction (Eq. (20)) and electric potential (Eq. (21)), one can see that the frequency dependence of the ratio of their modulus is completely determined by the frequency dependence of the ratio of current density and . In Appendix D, we show that this ratio is quasi-independent of frequency for a resistive medium, for low frequencies (smaller that ~10 Hz), and if the current sources are of very low correlation.
Thus, magnetic induction and electric potential can be very well approximated by:
for sufficiently small differential dendritic elements (N/l large).
where , and . Thus, the scaling of the PSDs of the electric potential and magnetic induction must be the same for low frequencies (smaller than ~10 Hz) if the medium is resistive and when the current sources have very low correlation.
According to Maxwell equations, the magnetic induction is given by:
where dv′=dx′1dx′2dx′3 and
for a perfectly resistive medium.
We now show that this expression is equivalent to Ampere-Laplace law.
From the identity , where , we can write:
Moreover, we also have the following identity
where is a unitary vector perpendicular to the integration surface and going outwards from that surface. Extending the volume integral outside the head, the surface integral is certainly zero because the current is zero outside of the head. It follows that:
where dv′=dx′1dx′2dx′3 because
Equation (27) is called the Ampère–Laplace law (see Eq. 13 in Hämäläinen et al. 1993). It is important to note that this expression for the magnetic induction is not valid when the medium is not resistive.
Finally, from the last expression, the magnetic induction for a differential element of dendrite can be written as:
In this appendix, we derive the differential expression for the electric field. Starting from Eq. (10), we obtain the solution for the electric potential:
It follows that the electric field produced by the ensemble of sources can be expressed as:
such that every differential element of dendrite produces the following electric field:
The transmembrane current obeys because we are in a quasi-stationary regime in a differential dendritic element. Taking into account the differential law of charge conservation , we have:
where is the density of transmembrane current per unit surface and δS′ is the surface area of a differential dendritic element. This approximation is certainly valid for frequencies lower than 1,000 Hz because the Maxwell-Wagner time (see Bédard et al. 2006b) of the cytoplasm () is much smaller than the typical membrane time constant of a neuron (τm ~ 5–20 ms).
Finally the contribution of a differential element of dendrite to the electric potential at position is given by
We note that the expressions for the electric field and potential produced by each differential element of dendrite have the same frequency dependence because it is directly proportional to for the two expressions. Also note that if the medium is resistive, then γf=γ and the frequency dependence of the electric field and potential are solely determined by that of the transmembrane current .
For each differential element of dendrite, we consider the standard cable model, in which the impedance of the medium is usually neglected (it is usually considered negligible compared to the membrane impedance). In this case, we have:
where , , , cm, rm et ri are respectively the membrane potential, the current density in the axial direction, the transmembrane current density, the specific capacitance (F/m2), the specific membrane resistance (Ω.m2) and the cytoplasm resistivity (Ω.m).
It follows that
Under in vivo–like conditions, the activity of neurons is intense and of very low correlation. This is the case for desynchronized-EEG states, such as awake eyes-open conditions, where the activity of neurons is characterized by very low levels of correlations. There is also evidence that in such conditions, neurons are in “high-conductance states” (Destexhe et al. 2003), in which the synaptic activity dominates the conductance of the membrane and primes over intrinsic currents. In such conditions, we can assume that the synaptic current sources are essentially uncorrelated and dominant, such that the deterministic link between current sources will be small and can be neglected (see Bédard et al. 2010). Further assuming that the electric properties of extracellular medium are homogeneous, then each differential element of dendrite can be considered as independent and the voltages Vm have similar power spectra.
In such conditions, we have:
Note that this expression implies that we have in general for each differential element of dendrite:
according to Eq. (34).
It follows that
Thus, for frequencies smaller than 1/(ωτm) (about 10 to 30 Hz for τm of 5-20 ms), the ratio will be frequency independent, and for each differential element of dendrite, we have:
for frequencies smaller than ~10 Hz.
1Examples of nonlinear effects are variations of the macroscopic conductivity σf with the magnitude of electric field E. Such variations could appear due to ephaptic (electric-field) interactions for example. In addition, any type of linear reactivity of the medium to the electric field or magnetic induction can lead to frequency-dependent electric parameters σ,ε,μ (for a detailed discussion of such effects, see Bédard and Destexhe 2009).
2Note that in textbooks, these linking equations (Eqs. (2)–(4)) are often algebraic and independent of time (for example, see Eqs. 5.2-6, 5.2-7 and 5.2-8 in Gulrajani 1998). The present formulation is more general, more in the line of Landau and Lifchitz (1984).
3In textbooks, the electric parameters are sometimes considered as complex numbers, for example with the notion of phasor (see Section 5.3 in Gulrajani 1998), but they are usually considered frequency independent.
4If it was not the case, then the source terms would be a function of the produced field, which would result in more complicated equations.
5Note that to compare scaling exponents between studies one must take into account that the electrode montage may influence the scaling. For example, in bipolar (differential) EEG recordings, if two leads are scaling as 1/(A+f) and 1/(B+f), the difference will have regions scaling as 1/f2.
N. Dehghani and C. Bédard are co-first authors.