The influences of tissue anisotropy and source activity on power and phase stability of low-frequency EEG rhythms: a mathematical observation of the forward problem model

Spontaneous low-frequency brain rhythms could be associated with many neural functions, which were usually evaluated by physical parameters, such as power, frequency, and phase. However, how the parameters might be affected by the electrical features of brain tissue and the electrical activity of the rhythm source is not yet clearly known. To address this issue, the electrical significance of power and narrow-band phase stability (NBPS) was investigated by simulated rhythms. The rhythms were derived from the oscillatory field potentials (FPs) on a homogeneous sphere model of the solution to the electroencephalogram (EEG) forward problem. The sphere’s electrical feature was set as isotropic or anisotropic conductivity. The source was set as a quasi-static dipole current, whose activity was representative of a low-frequency sine oscillation with a nonlinear phase course, and the source location was changeable. After the instantaneous power and phase of simulated rhythms were estimated by the Hilbert transformation, the NBPS was calculated and the statistical properties of power and NBPS were analyzed. It was found that only nonlinear phase dynamics could lower the NBPS. However, power depended on many factors, such as conducting anisotropy, amplitude and positon of dipole current, and meshes on the sphere model. We hypothesized that NBPS would map the influence of nonlinear phase dynamics on the brain rhythms, but be independent of power. This research might highlight the nonlinear phase dynamics of intrinsic low-frequency oscillations in the brain. The results might be beneficial for the measurement and analysis of spontaneous EEG rhythms.


Introduction
Brain rhythms in distinct frequency bands, such as delta (0.5-4 Hz), theta(4-8 Hz), alpha(8-12 Hz), beta (12-30 Hz) and gamma(30-80 Hz) could play an important role in a variety of brain functions (Buzsaki 2006). In neural networks, the rhythms reflected the summated excitatory post-synaptic potentials of tens of thousands of coherently active neurons (Lopes da Silva 1991). In fact, the oscillations can be derived from EEG at many levels, e.g., scalp EEG, deep brain EEG in vivo, or intracranial EEG. Some parameters (e.g., power, mean phase, mean frequency, phase synchronization, and phase resetting) were widely utilized to investigate the scale of synchronized neurons, neural communication across areas, stimulus processing, memory formation and cognitive control of information input, as well as human and animal behavior (Winson 1978, Lachaux et al 1999, Mormann et al 2005, Schnitzler and Gross 2005, Ali 2006, Uhlhaas and Singer 2006, Kraskov et al 2007, Chauviere et al 2009, Wang 2010, Ding and Simon 2013, Ge et al 2013. Recently, low-frequency oscillations, particularly oscillations in the delta and theta frequency bands, have attracted more attention , Ng et al 2013. Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Traditionally, broad-band EEG signals have been decomposed into narrow-band oscillations using spectral analysis with complex wavelets, the Hilbert transformation, or filtering (Le Van Quyen et al 2001, Yeung et al 2004, Ding and Simon 2013. Neuronal oscillations could be characterized by three distinct parameters: amplitude, frequency, and phase. The amplitude, or power fluctuations, could reflect the extent of synchrony of the neurons in a local assembly . Frequency was determined by the physiological time constants of the relevant neuronal assemblies (Buzsaki 2006, Wang 2010. The phase pattern might reflect neuronal synchronization and reveal the coordination or communication of neuronal activity across brain regions, providing more information than the amplitude could be (Schyns et al 2011, Ng et al 2013. While there has been extensive research on broadband EEG signal processing, the significance of the phase dynamics of narrow-band oscillations has not been clearly understood. Early in the 1950s, narrow-band phase coherence was applied to brain rhythms (Wiener 1956(Wiener , 1957. This research indicated that nonlinearity might affect the interaction between sources of brain rhythms, so that these sources tended to oscillate at the same frequency. Recently, the information flow observed by the narrow-band method was suggested to be related to synaptic coordination (Zheng et al 2011, Zheng et al 2012, and a phase-based entropy in the narrow-band was claimed to be suitable to estimate directed connectivity in large-scale investigations of the human functional connectome (Lobier et al 2014). Using deep brain EEG, we found that both theta power and theta NBPS were related to synaptic reorganization during the development of epilepsy (epileptogenesis) (Ge et al 2013).
We hypothesized that NBPS could be used to map the influence of nonlinear phase dynamics on the brain rhythms, whose significance might be different from that of power. Therefore, we explored this issue in key areas: (i) Simulations of low-frequency rhythms; (ii) Evaluation of simulated rhythms by power and NBPS; and (iii) Global analysis of the properties of power and NBPS in some comparative studies, such as linear and nonlinear phase dynamics, and isotropic and anisotropic conductivity values.
The theory of equivalent dipole current has been widely used in the generation of EEG (forward problem) (Hjorth 1975, Nunez 1981, Lopes da Silva and Van Rotterdam 1982, de Munck et al 1988, Kearfott et al 1991, Yan et al 1991, Kim et al 2003. Given the position and moments of a static equivalent dipole current generator, and the geometry and electrical conductivity profile of the volume conductor (i.e., model of the head), the electrical FPs were usually regarded as EEG, which could be theoretically stated by the Poisson's equation and the Neumann boundary condition on the head surface, and traditionally computed by the finite element method (FEM) (Nunez et al 1997, Gulrajani 1998, Haueisen et al 2002, Zhang et al 2006. When referring to the solution to the EEG forward problem, the head was regarded as a homogeneous sphere model, with the electrical features of brain tissue described by (an)isotropic electrical conductivity, and the brain rhythm source depicted as a quasi-static dipole current, whose activity was representative of a sine oscillation at a low frequency. Thus, by changing the conductivity tensor and including anisotropy conductivity, amplitude, phase dynamics, and the placement of the dipole current, respectively, a distribution of oscillatory FPs produced by the dipole current at a time-point could be estimated by the FEM on the sphere model. During a time period, a simulated rhythm could be formed by the time series of oscillatory FPs.
In this study, the instantaneous power and phase of simulated rhythms were estimated by the Hilbert transformation. Then, all the values of power and relative phase on the sphere model were averaged over time. Finally, statistical tests were employed to evaluate the electrical significance of power and phase stability, respectively. To highlight the electrical significance of NBPS, we theoretically investigated the influence of electrical features of brain tissue and the electrical activity of the brain rhythm source, and compared it to the electrical significance of power, in order to know whether NBPS was independent of power. It was noted that three statistical methods were used to depict NBPS, including a histogram of the mean relative phase on a rose plane, the value of the mean relative phase, called the phase preserved index (PPI), and relative phase sorting strips with the help of EEGLAB.

Methods
The power and phase stability of the simulated lowfrequency rhythms were analyzed by changing some factors: the conductivity tensor and inclusion of anisotropy conductivity, the activity of dipole current, i.e., amplitude, frequency and phase, duration, and placement of the dipole current. When a factor was analyzed as a variable, the other factors were defaults in a computation.

Data derivations of low-frequency rhythms
The rhythms were simulated by oscillating electrical FPs in a 3D sphere conductor (radius=10 cm, homogeneous medium), which were estimated by the FEM in the software COMSOL Multiphysics. Given a sine dipole current fixed at the center (the default), and 2, 4, 6 and 8 cm from the center, along the X-axis (the negative Xaxis was the dipole moment). The dipole current activity could be described by three factors: amplitude (1 ( nA the default), 2 ) nA , frequency (6 Hz), and initial phase-0 radian (the default), 4t 2 radian, t was the time-during a time window (1050 ms as a default, 350 ms). The electrical field produced by the dipole current was quasistatic, and the distribution of FPs driven by the dipole current at each time-point were produced on the generated meshes (a total of 3611 meshes were automatically generated under the mode of user-controlled mesh/standard in the software). A 1000 ms time course of dipole current was split into 128 time points. The oscillating electrical potentials could be regarded as the simulated rhythm on each generated mesh during a time course. It should be noted that the initial phase of 0 radian and 4t 2 radian of dipole current could, respectively, produce the linear and nonlinear phase dynamics of simulated rhythms, assuming the conductivity was constant.
The meshes (i.e., placement of a recording electrode, location of the dipole current, zero potential surface, and iso-potential lines) were illustrated in figure 1.

Conductivity
In the brain, the conducting medium was ideally isotropic ( ) s s s = = = -0.14 S m y z x 1 (Wolters et al 2006), and the conductivity was presumed to increase equably from -0.14 10 S m 1 / to -0.14 S m 1 (Here, called an isotropic medium.). However, the conducting medium was really anisotropic was still fixed at -0.04427 S m 1 (e.g., increase in the ratio s x :s z from 1:2 to 1:9) (Here, called anisotropic medium 2). The influence of anisotropic conductivity was compared to that of isotropic conductivity.

Power and phase stability
The instantaneous power and phase of simulated rhythms were estimated by the Hilbert transformation. On each mesh, after instantaneous power was averaged over time, the power change with the meshes could be observed as a spherical distribution. The subsequent statistical significance was performed between two groups of relative phase values corresponding to nonlinear phase dynamics and linear phase dynamics, anisotropic conductivity and isotropic conductivity, and two positions of dipole current, respectively. Thus, the power change for the conductivity could be observed globally. In addition, in order to make this apparent, the power values were further averaged over the meshes, and were displayed in a histogram. A student's unpaired t-test was used to test significant differences. NBPS was estimated by the statistics on the relative phase values defined in a circular space at the frequency of interest (i.e., 6 Hz): First, the relative phase values of simulated rhythms on all the meshes and time points were displayed in a circular space by the rose statistics method (rose, Matlab). The level of significance was indicated by the scale of the rose petals in the concentrated part (called the preferred relative phase). Obviously, it could represent the phase stability, with a greater scale implying more stable (or less split) phase dynamics, and vice versa. Here, the phase stability (denoted as PS) was defined in a normalized way as the ratio of the scale of preferred relative phase to the total.
Second, a method of phase sorting in EEGLAB was employed to display phase stability (Delorme and Makeig 2004). Here, the relative phase (as shown in equation (1)) was first embedded in a normalized rhythm (a simulated rhythm was divided by its amplitude) by fast Fourier transform (FFT) and inverse fast Fourier transform (iFFT), and then an image toolbox of phase synchronization (Phase sorting, EEGLAB) was performed on the normalized rhythm. The strips showed the sorting of the relative phase, with the width of the strips or the sorting concentration, as indicated by the color bar. The wider the strips or the greater the color-bar value, the more stable the phase stability was, and vice versa. The time course of phase stability could be observed intuitively by this technology.
A simulation of the phase stability depicted by two kinds of methods was shown in figure 2. Assuming a 6 Hz sine function as a reference signal with a time course of phase: 2π×6×t, two simulated signals of a 6 Hz sine function were produced with a time course of phase: 2π×6×t, 2π×6×t+4t 2 ; that is to say, the time course of the relative phase of two simulated signals was, respectively, 0 radian and 4t 2 radian, as observed in figure 2(A) (Here, the 1000 ms time course was split into 128 time points). The phase stability represented by PS in the rose plane is shown in figure 2(B). By adding a tiny phase bias (i.e., kk×0.01 radian, kk is the number of oscillator, kk=1,2,3, KK100) to the above time courses of 0 radian and 4t 2 radian, respectively, 100 oscillators with similar phase dynamics were formed, which were then sorted by EEGLAB as shown in figures 2(C) and (D) respectively. In figure 2(C), the width of the sorting strips was uniformly wide during the whole time course, implying greater phase stability. Inversely, the width of the sorting strips became gradually narrow with the time course, implying less phase stability as shown in figure 2(D).
Third, NBPS was described by the mean relative phase value Φ over a time window of simulated rhythms (denoted as phase preserved index, PPI), as shown in equation (2): In the three methods, NBPS was in the range of 0-1, where 1 indicated perfectly stable phases, and 0 indicated independent phases. The bootstrap statistics (SPSS) on all relative phase values showed the significance. (2) The power values became monotonic in the isotropic medium, however, they dispersed in the anisotropic medium; and (3) Overall, power decreased with increasing conductivity (p<0.001). Figure 3(B) showed the spherical distribution of power values in relation to altering the anisotropy extent. It was found that the power values generally rose with the anisotropic extent increasing.

Influence of amplitude of dipole current
By doubling the amplitude of the dipole current, the spherical distribution of the power values was shown in figure 3(C). It was found that the power values were proportional to the square of the amplitude of the dipole current, not related to the phase dynamics of the dipole current (not shown here).

Influence of location of dipole current
By altering the placement of the dipole current, the spherical distribution of the power values and FPs were shown in figures 4 and 5, and the histograms of the power values in relation to the conductivity and the location of the dipole current were shown in figure 6. The findings were that: (1) Both greater power values and greater FPs were concentrated near the dipole current as shown in figures4 and 5; (2) The more off-central the placement of the dipole current, the greater the power values around the dipole current were, as shown in figures 4 and 6.
These results implied that the electrical significance of the power of brain rhythms measured by EEG might be complicated and that they might vary greatly with many factors, such as placement of the recording electrodes, anisotropy including the anisotropic extent of the brain tissue, and the strength and location of the rhythm source.
3.2. Phase stability of simulated rhythms 3.2.1. Influence of conductivity By altering the conductivity value, the phase stability as indicated by the PPI changed with the meshes (in figure 7(A)). It was found that the PPI values were related neither to the (an)isotropic medium conductivity nor to the meshes. The finding converged with the phase stability, as indicated by PS in the rose figure (on the above panels in figures 7(B), (C)), and it also converged with the phase stability, as indicated by the sorting strips (on the below panels in figures 7(B), (C)).

Influence of initial phase of dipole current
The reduction in PPI values (from 1.0 to 0.372) (shown in figure 7(A)) was found to converge with the decline in phase stability values (from 1.0 to 0.372), presented by the rose petals (on the above panels in figures 7(B), (C)), and with the decrease in relative phase sorting concentration (from 0.45 to 0.39) (on the below panels in figures 7(B), (C)), presented by the color-bar values when the phase dynamics of the simulated rhythms changed from linearity (corresponding to the initial phase with the 0 radian of the dipole current, as described in the Methods) to nonlinearity (corresponding to the initial phase with the 4t 2 radian of the dipole current, as described in the Methods). The results indicated that phase stability could be sensitive to the nonlinear phase dynamics of simulated rhythm.

Influence of location of dipole current
By altering the location of the dipole current (here, the initial phase and duration of the dipole current were defined as the 4t 2 radian and 350 ms, respectively), the PPI values changed with the location of the dipole current, as shown in figure 8(A). It was found that phase stability was not sensitive to the location of the dipole current.

Influence of duration of dipole current
There was little reduction in the phase stability of simulated rhythms (PS) (from 1 to 0.994), as shown in figure 8(B), and there was little change in the phase stability of simulated rhythms, according to relative phase sorting concentration (around 0.88 in the colorbar values) as shown in figure 8(C), when the phase dynamics of simulated rhythms changed from linearity (corresponding to initial phase 0 radian of the dipole current, as described in the Methods) to nonlinearity (corresponding to initial phase 4t 2 radian of the dipole current as described in the Methods).
Given the 4t 2 radian initial phase of the dipole current, the reduction in PPI changed from 0.978 (shown in figure 8(A)) to 0.372 (shown in figure 7(A)) when the duration of the simulated rhythms increased from 350 to 1050 ms. A convergent finding was the  (C)). The results indicated that phase stability decreased monotonically with the duration of the dipole current. In addition, it was necessary that the time window was long enough to observe the nonlinear phase dynamics.
In summary, NBPS seemed to be sensitive only to phase dynamics; the nonlinearity of the phase dynamics might reduce it. In the computation, the longer the time window, the lower the value of NBPS was; that is, an optimized long enough time window was necessary for NBPS to reflect the nonlinear phase dynamics.

Conclusions and discussions
The electrical significance of NBPS and power were analyzed based on simulated low-frequency oscillatory FPs. It was concluded that NBPS was related only to nonlinear phase dynamics of simulated rhythms.
However, power of an electrode might be substantially affected by many factors, such as strength of the rhythm source, distance from the rhythm source, and anisotropic conductivity.

Phase stability
The EEG phase patterns might reflect the selectivity of neural firing (Ng et al 2013). Some researches indicated that the mean phase was related to nonlinearity in the neural system. Many factors affected nonlinearity, such as electrophysiology, pathology in the brain, and external stimuli (Lachaux et al 1999, Mormann et al 2000, Ding and Simon 2013. The conclusions derived from the mean-phase-based model (in equations (1) and (2)) were supported by these investigations. In addition, the conclusion on phase stability agreed with the earlier hypothesis: nonlinearity will affect the interaction between the rhythms sources, leading to an alteration in narrow-band phase stability (Wiener 1956(Wiener , 1957. Finally, the theory of a self-sustained system argued that the weak coupling could affect the phase synchrony of the oscillators, regardless of the amplitude of the oscillators (Rosenblum et al 1996, Pikovsky 2001, Rosenblum et al 2002). The conclusion that the phase stability was independent of power was supported by this theory.

Anisotropy conductivity
Brain tissue was composed of gray matter and white matter, which could define the feature of conductivity. Anisotropy could be measured by some technologies, such as electromagnetically induced transparency (EIT) and diffusion tensor imaging (DTI) (Tuch et al 2001, Abascal et al 2008, and it could be associated with physiology and pathology in the brain (McAdams and Jossinet 1995, Haueisen et al 2002, Wolters et al 2006. Given these findings, we produced many values of anisotropy conductivity in the simulation. In addition, we only considered the influence of nonlinear phase dynamics on NBPS, not considering the influence of nonlinearity on either conductivity or dielectric constant. However, we predicted that NBPS would be lowered if the nonlinearity of either conductivity or dielectric constant could make the phase split. Finally, in spite of the fact that NBPS was investigated at a given frequency in the current study, nonlinear frequencies could produce results on phase dynamics similar to those described in this paper.

Power
In line with Ohm's law, the FPs should be inversely proportional to electrical conductivity. As power should reflect the strength of simulated rhythms, therefore, it was inversely proportional to conductivity. The influence of (an)isotropic conductivity on power was convergent with the analytical solution of electrical field conduction. In the solution, the FPs produced by the dipole current were proportional to the amplitude of the dipole current, but inversely proportional to the square of the distance to the dipole current (Brody et al 1973). Given the isotropic conductivity, despite the fact that the electrical conduction was the same in a 3D space, the FP was only determined by the distance between a mesh when the location of the dipole current was fixed. So, the monotonic curves of power variation with the meshes (shown in figure 3) were in accord with the analytical solution. Given the anisotropic conductivity, the electrical conduction was different in 3D directions, therefore, the electrical field vector sum was very different in a 3D space, implying various values of FPs, so the shapes of the curves of power change with the meshes were very divergent. In order to explain the influence of the location of the dipole current by the solution, we estimated that the greater the eccentricity of the dipole current, that is the closer the distance between a mesh and the dipole current, the stronger the rhythms were, therefore the greater the power values that could be observed, as shown in figures 4 and 5.
In addition, as the FPs were proportional to the amplitude of the dipole current, and the power value should be proportional to the square of the amplitude of the dipole current, so the power value was proportional to the amplitude of the dipole current (shown in figure 3(C)). Finally, in the neural measurement, spontaneous power should reflect the strength of the brain rhythms, which were related to many factors, such as brain region, scale of synchronous neurons, brain function, pathology, and physiology (Mormann et al 2005, Buzsaki 2006, Chauviere et al 2009. The finding that the power value depended on various factors in the simulation offered a theoretical support to these experimental measurements.

The dipole current
The concept of a 'dipole current' was abstract, in that it was not as tangible as an actual current source. In addition, the dipole moment was very complex in a real neural system. In a typical regular arrangement, the dipole moment was aligned with the main axis of the cell, and the magnitude of the dipole moment was thought to reflect the longitudinal current flowing along the somato-dendritic axis as a function of the dendritic and somatic membrane. In a population of neurons, the summation of IPSCs and EPSCs could be equivalent to a dipole current, leading to the generation of extracellularly recorded FPs, which could only be measured by macro-electrodes at a distance from the dipole current (Lopes da Silva and Van Rotterdam 1982). In addition, different types of neuron pools could generally correspond to different fields, such as closed, open, and open-closed fields in the central nervous system (Lopes da Silva and Van Rotterdam 1982). Here, we employed an open field to examine the electrical significance of power and phase stability. In addition, for the sake of intuitive observation, we not only designed the amplitude, phase dynamics, and location of the dipole current, but also directly linked the dipole moment to the X-axis. However, the work was not suitable for a dipole current oscillating at a high frequency. (B) phase stability indicated by preferred rose petal (above panel) (bootstrapping technique; * p<0.05), and relative phase sorting indicated by strips (below panel) in the isotropic medium; (C) phase stability indicated by preferred rose petal (above panel) (bootstrapping technique; * p<0.05), and relative phase sorting indicated by strips (below panel) in the anisotropic medium.

Simulations and real EEG rhythms
An inhomogeneous head model (called the four shells head model) can be divided into multiple layers, such as scalp, skull, cerebrospinal fluid and cerebral cortex. In a homogeneous model, only the cerebral cortex was considered. However, the shape of the real cerebral cortex was complicated. Furthermore, the mechanism of EEG rhythms in an actual neural system has not yet been clearly understood. The simulations on EEG rhythms were approximated to a great extent. Finally, the real EEG signals could contain noise. However, considering that the oscillatory FPs were computed at a special frequency, we could not mimic the effect of noise in a narrow frequency band.
Further work will be undertaken to compare other similar models of phase stability to this model.