Cross-frequency transfer in a stochastically driven mesoscopic neuronal model

1. Introduction

Brain activity, as registered in macroscopic recordings like EEG, does not yield power spectra composed of multiple isolated and narrow frequency peaks, but of broad frequency bands that are merged with each other in a seamless manner, and which are embedded in a 1/f^b background with large power at low frequencies and a fat tail at high frequencies (Freeman et al., 2000; Buzsáki and Draguhn, 2004). Nonlinearities in the interactions between neuronal populations (Friston, 2000) can be expected to lead to mixing within this continuum of frequencies. Indeed, transcranial stimulation of the brain at low frequencies (smaller than 1 Hz) has been seen to cause for instance an increase in oscillatory power at larger frequencies (5–10 Hz) (Marshall et al., 2006; Massimini et al., 2007). Much emphasis has been placed recently on quantifying and characterizing the transfer of spectral power across frequencies (known in what follows as cross-frequency coupling, Jirsa and Müller, 2013), and on identifying its functional roles in the brain (Jensen and Colgin, 2007; Canolty and Knight, 2010). In particular, power spectrum correlations have been observed for instance between theta and gamma rhythms in the rat hippocampus during memory retrieval (Shirvalkar et al., 2010), between posterior gamma and frontal alpha/beta oscillations in the human brain during motor imagery tasks (de Lange et al., 2008), and between the gamma and delta bands in different regions of the human visual cortex during a visual task (Bruns and Eckhorn, 2004). Phase coupling between delta and alpha bands has also been reported in human brains performing an orientation task (Isler et al., 2008). More common is the situation in which the oscillation power in a given frequency band is modulated by a second rhythm at lower frequency. Such phase-to-amplitude cross-frequency coupling has been observed for instance between alpha and gamma activities in humans during rest (Osipova et al., 2008) and between theta and gamma oscillations in rats during learning (Tort et al., 2009). Other behavioral correlates of cross-frequency coupling have been found, associated for instance with reward coding (Cohen et al., 2009a) and decision making (Cohen et al., 2009b) in humans. Also, recent work has shown that such cross-frequency coupling is modulated by behavioral tasks (Voytek et al., 2010). Cox et al. (2014) reported cross-frequency coupling between the phase of sleep spindles and the amplitude of higher frequency rhythms, in particular beta, recorded in EEG during sleep. This effect, in turn, was modulated in the frontal cortex by the phase of slow sleep oscillations.

Despite the large number of experimental studies pointing toward cross-frequency correlations, several difficulties arise when it comes to the interpretation of this phenomenon. As pointed out by Aru et al. (2015), the methodologies applied in a number of recent studies on cross-frequency coupling are not flawless and the results might have been overinterpreted. Therefore, further and stricter studies on the functional role of cross-frequency coupling are needed to confirm previous results. In particular, not all cross-frequency correlations are signatures of direct interaction between rhythms. When they are, such correlations may be explained by different mechanisms, which may be grouped into two broad scenarios. In one scenario, two neuronal oscillators operating at two different rhythms might be coupled bidirectionally to each other. This coupling could mediate an interaction that would result in each of the oscillators being affected in one way or another by the natural frequency of the other oscillator (Jirsa and Müller, 2013). When such bidirectional interaction occurs locally, it has been proposed to be mediated by the firing activity of the underlying neurons (Mazzoni et al., 2010). In such a way delta oscillations, for instance, control the level of local cortical excitability, which in turn modulates the excitatory-inhibitory balance that gives rise to the gamma rhythm (Mazzoni et al., 2011). In a second, somewhat simpler scenario, cross-frequency correlations might arise due to unidirectional coupling, through which the spectral features of the driving neuronal population would be directly transferred to the driven population. When the same external stimulus is encoded by two different rhythms, cross-frequency correlations can appear as a result of that common unidirectional driving (Mazzoni et al., 2008). In some cases, however, the stimulus does not necessarily affect directly the neurons underlying one of the rhythms. This might be the case of recent experimental work by Bayer et al. (2011), who examined the effect of rocking on sleep in human subjects. In that study, healthy volunteers were asked to lie down on a rocking bed that oscillated slowly, at a frequency of 0.25 Hz. This periodic stimulation was seen to ease the transition from waking to sleep, and to increase the power of cortical oscillations (measured via EEG) in the alpha range. Here we ask whether a cross-frequency transfer such as that reported by Bayer et al. (2011) can be the result of the low frequency input driving a mesoscopic broadband oscillator operating in the alpha range. To that end, we need a mesoscopic model of brain activity.

Brain dynamics at the mesoscale is frequently described by population models such as the neural mass model (NMM), originating from the works of Freeman (1972), Wilson and Cowan (1972), Amari (1974), Lopes da Silva et al. (1974), and Nunez (1974). These models aim to reproduce the average behavior of relatively large populations of cells. The dynamical unit in this model can be interpreted as a cortical hypercolumn (Jansen et al., 1993; Jansen and Rit, 1995), the model variables being the average postsynaptic potentials of the different neuronal populations (Faugeras et al., 2009). NMMs have been extensively used in recent years to describe a wide variety of brain behaviors including rhythm generation (Ursino et al., 2010) and propagation (Cona et al., 2011), spontaneous dynamics (Nguyen Trong et al., 2012), photic stimulation (Spiegler et al., 2011), criticality (Aburn et al., 2012), and even plasticity (Wang and Knösche, 2013). Aberrant dynamics in epilepsy has been described with NMMs both at the level of the generation (Wendling et al., 2000) and termination (Goodfellow et al., 2012; Freestone et al., 2013) of epileptic seizures. Coupled NMMs have also been used to examine the relationship between structural and functional connectivity in healthy and neurodegenerative conditions (Pons et al., 2010; Ponten et al., 2010). Generalizing the NMM description to continuous space gives rise to neural field models (Jirsa and Haken, 1996; Coombes, 2010), which have been used to study the spatiotemporal dynamics of cortical waves (Hutt and Atay, 2006; Bojak and Liley, 2010) and even to implement control strategies in robotics (Erlhagen and Bicho, 2006).

NMMs can be made to exhibit a variety of rhythms depending on the values of the dynamic parameters (David and Friston, 2003), but in most situations the oscillations obtained are relatively narrowband. This contrasts with the 1/f^b spectra usually observed experimentally (Freeman et al., 2000). In neural fields, such type of spectral behavior has been linked to the multiple spatial scales characteristic of spatially extended neuronal tissue, where it is observed near a stability threshold (Hutt and Frank, 2005). In this paper, we reproduce the spectral properties measured by Bayer et al. (2011) by applying a temporally correlated noise to a coupled NMM. This is a reasonable assumption, since the brain has multiple sources of noise (Faisal et al., 2008) that have a variety of functional roles (McDonnell and Ward, 2011). In microscopic models, a temporally correlated Ornstein-Uhlenbeck noise is known to reproduce the observed 1/f spectral profile of LFP activity (Sancristóbal et al., 2013). Our results show that a network of coupled neural masses subject to temporally correlated noise exhibits a well-defined rhythm (in the alpha range) embedded in a broadband spectral background similar to what is observed experimentally. We also show that this broadband oscillator reacts to periodic driving at a frequency much lower than its natural frequency, by increasing its activity at the latter in agreement with experimental observations.

2. Materials and Methods

2.1. Extending Jansen and Rit Model

The basic building block of our model is a cortical column (Hubel and Wiesel, 1977; Helmstaedter et al., 2007; Ts'o et al., 2009) that we describe in the way proposed by Jansen et al. (1993). In this model the neurons of a cortical column are classified into three neuronal populations: pyramidal neurons, excitatory interneurons, and inhibitory interneurons. The dynamics of each population is described using two simple transformations. The first one stands for synaptic processing: it describes how the presynaptic signal coming from interconnected populations translates into a postsynaptic membrane potential. This transformation is linear and is given by the convolution:

$\begin{array}{cc}\begin{array}{l}y(t)={\displaystyle {\int }_{-\infty }^{t}h}(t^{\prime })p_{\text{tot}}(t-t^{\prime })\text{d}{t}^{\prime },\hfill \end{array}& \left(1\right)\end{array}$

where p_tot(t) is a total input acting upon the population, expressed in terms of a firing rate, y(t) is the net postsynaptic membrane potential, and h(t) is an impulse response. The kernel of the transformation is valid for t > 0 and is defined for excitatory and inhibitory connections as follows:

$\begin{array}{cc}{h}_e(t)=Aat{e}^{-at}& \left(2\right)\end{array}$

$\begin{array}{cc}{h}_i(t)=Bbt{e}^{-bt},& \left(3\right)\end{array}$

where A and B are the maximum excitatory and inhibitory postsynaptic potential amplitudes, respectively, and a and b are inverse time constants that lump together all signal propagation delays. This transformation can be expressed in differential form using the Laplace transform, which leads to:

$\begin{array}{cc}\frac{d^{2}y(t)}{d{t}^2}+2a\frac{dy(t)}{dt}+a^{2}y(t)=Aa·p_{\text{tot},\text{exc}}(t),& \left(4\right)\end{array}$

where y(t) is a postsynaptic membrane potential averaged over all neurons in a given population, and p_{tot, exc} is the sum of firing rates of all excitatory signals coming into that population. An equivalent expression containing the sum of all incoming inhibitory signals p_{tot, inh}(t) with constants B and b describes inhibitory processing.

The second transformation describes how the postsynaptic membrane potential within a population translates into a firing rate in that population's output. This transformation is nonlinear and is given by:

$\begin{array}{cc}\text{Sigm}(y)=\frac{2e_0}{1+e^{r({\nu }_0-y)}},& \left(5\right)\end{array}$

where 2e₀ is the maximum firing rate, ν₀ is the potential for which the firing rate is equal to half of the maximum, and r determines the steepness (and nonlinearity) of the response.

Figure 1 shows a schematic representation of the model. It is built as a network of coupled cortical columns. Within a single column a population of pyramidal neurons feeds forward to the excitatory and inhibitory interneuron populations, which in turn feed back to the pyramidal neurons leading to positive and negative feedback, respectively. The excitatory and inhibitory interneurons receive only excitatory input from the pyramidal neurons residing in the same cortical column. The pyramidal population receives an inhibitory input from the inhibitory interneurons, and an excitatory input [denoted with p^pyr_{tot, exc}(t)] coming from both the excitatory interneurons and from sources external to the column. The latter part we indicate with p^pyr_{ext, exc}(t), which for a cortical column i is given by:

$\begin{array}{cc}{p}_{\text{ext},\text{exc}}^{\text{pyr},i}(t)=p_{\text{const}}+{\displaystyle \sum _j{p}_{\text{coup}}^{j,i}(t)+p_{\text{osc}}(t)}+{\xi }_{\text{ou}}^i(t).& \left(6\right)\end{array}$

FIGURE 1

Figure 1. Driving inputs and connectivity between neuronal populations within a single cortical column and between columns. The model is formed by a network of coupled cortical columns described by the Jansen-Rit model. In each column a population of pyramidal neurons (green triangle), feeds forward to a population of excitatory interneurons (blue hexagons) and inhibitory interneurons (red circles). The two interneuron populations feed back into the pyramidal neurons. Excitatory connections are marked with blue arrows, and inhibitory connection with red lines with circular endings. The model is homogeneous; all columns are identical and all inter-column connections have same strength. Columns are fed with a common deterministic input and independent realizations of an Ornstein-Uhlenbeck noise. The deterministic input denoted here by a sinusoid has a constant DC level.

Here p_const is a constant component and ξⁱ_ou(t) is a stochastic one. Those two components stand for a contribution from sensory input or brain areas that are not explicitly included in the model. The component ∑_j p^{j, i}_coup(t) is the summed contribution from other cortical columns (indexed with j) connected with column i, and the oscillatory component p_osc(t) = Ã sin (2πft) stands for a periodic stimulus.

The system of equations that describe the dynamics of a single cortical column i is given by:

$\{\begin{array}{l}\begin{array}{cc}{\ddot{y}}_0^i(t)+2a{\dot{y}}_0^i(t)+a^2{y}_0^i(t)=Aa\text{Sigm}[y_1^i(t)-y_2^i(t)]& \left(7\right)\end{array}\\ {\ddot{y}}_1^i(t)+2a{\dot{y}}_1^i(t)+a^2{y}_1^i(t)=Aa\{p_{\text{ext},\text{exc}}^{\text{pyr},i}(t)\\ \begin{array}{cc}+C_2\text{Sigm}[C_1{y}_0^i(t)]\}& \left(8\right)\end{array}\\ \begin{array}{cc}{\ddot{y}}_2^i(t)+2b{\dot{y}}_2^i(t)+b^2{y}_2^i(t)=Bb\{C_4\text{Sigm}[C_3{y}_0^i(t)]\}& \left(9\right)\end{array}\end{array}$

where i runs from 1 to the number of columns N. y₀ is the excitatory postsynaptic membrane potential that feeds into the two populations of interneurons and y₁ and y₂ are excitatory and inhibitory postsynaptic membrane potentials that enter into the pyramidal population, respectively. C_k (k = 1, 2, 3, 4) are constants representing the connection strengths between populations.

The pyramidal neurons are known to be the main source of the EEG signal, which locally is proportional to the difference between their excitatory and inhibitory potentials. In the notation introduced above this value is expressed as y₁(t) − y₂(t). Below we will analyze the model behavior in terms of this quantity, in order to compare it with experimental data. Model parameters are listed in Table 1.

TABLE 1

Table 1. Parameter values of our neural mass model.

The Jansen-Rit model is capable of exhibiting two different types of periodic behavior: alpha oscillations and spiky oscillations (Grimbert and Faugeras, 2006; Spiegler et al., 2010). The latter are characterized by a lower frequency and higher amplitude than the alpha dynamics. The dynamical regime in which the system operates depends on the external excitatory input p^pyr_{ext, exc}. For the parameters used here and in the sole presence of the constant input component (p^pyr_{ext, exc} = p_const), the model is known (Grimbert and Faugeras, 2006) to undergo a Hopf bifurcation at an input value p_const = 89.83 Hz, above which a limit cycle appears that corresponds to the alpha oscillations. For 113.58 > p_const > 137.38 Hz, the alpha regime coexists with spiky oscillations, whereas for p_const > 137.38 Hz the alpha regime becomes the only attractor again. Finally, for p_const = 315.70 Hz the system undergoes another Hopf bifurcation and the limit cycle collapses back to a fixed point.

Experimental observations reveal a large autocorrelation time in EEG signals (see Aburn et al., 2012, and references therein). A number of studies interpret this observation as an instance of critical behavior such as the one found in the proximity of second-order phase transitions (Chialvo, 2010; Deco et al., 2013). This suggestion is, however, still under debate; for example Bédard et al. (2006) explained one of the putative signatures of criticality, 1/f scaling of the EEG power spectrum, in a way that does not rely on critical phenomena, but on filtering properties of the brain's tissue. In our simplified Jansen-Rit model description, we reproduce this spectral feature by operating close to one of the two Hopf bifurcations delimiting the alpha regime (Aburn et al., 2012). This allows us to see an emergent, but not fully developed, alpha resonance, characteristic of the proximity of a transition to an oscillatory regime (Kang et al., 2010; Battaglia and Hansel, 2011). We set the external excitatory input parameters in such a way that its value averaged over time (and over columns in case of the network of columns), 〈p^pyr_{ext, exc}〉 ≃ 90 Hz, is located close to the first Hopf bifurcation point. The external input p^pyr_{ext, exc} delivered to each cortical column contains constant, stochastic and periodic components, as well as an input coming from coupled columns. Only the constant component p_const and the contribution from the afferent columns ∑_j p^{j, i}_coup(t) have non-zero mean, and therefore determine 〈p^pyr_{ext, exc}〉. The columns receive a common periodic signal mimicking a sensory stimulus received by cortical areas from the thalamus. Stochastic and periodic components have zero means and thus, even though they do not affect the average of the total input 〈p^pyr_{ext, exc}〉, they do contribute to its variance.

As mentioned in the Introduction, the brain contains sources of noise originating from different mechanisms (Faisal et al., 2008). Noise has been taken into account in past studies of the Jansen-Rit model (Jansen and Rit, 1995; Pons et al., 2010; Aburn et al., 2012) usually in the form of white noise. Here we use Ornstein-Uhlenbeck noise, which has a finite correlation time and is a more realistic representation of background synaptic noise. Destexhe and Rudolph (2004), for instance, showed that Poisson spike trains acting upon a neuron lead to a temporal correlation in membrane conductivity fluctuations, which under certain conditions can be modeled with Ornstein-Uhlenbeck noise. This noise, corresponding to ξ_ou variable in Equation (6), is generated by the following linear stochastic differential equation:

$\begin{array}{cc}\begin{array}{ll}\frac{d{\xi }_{\text{ou}}}{dt}\hfill & =-\frac{{\xi }_{\text{ou}}}{\tau }+\frac{\sqrt{2D}}{\tau }{\xi }_w(t)\hfill \end{array}& \left(10\right)\end{array}$

where ξ_w(t) is a random variable representing Gaussian white noise with zero mean, 2D defines the amplitude of the stochastic component and τ is the correlation time of the Ornstein-Uhlenbeck noise. Each column in the model was fed with an independent realization of the Ornstein-Uhlenbeck noise.

2.2. Network of Coupled Cortical Columns

The full model is composed of a number of cortical columns, each modeled in the way described above and connected with each other via the pyramidal neurons. For simplicity we choose all-to-all bidirectional connectivity and do not consider delay in the coupling (see Figure 1). The latter is justified by the fact that neighboring cortical columns are separated by a distance smaller than 1 mm and signal speed propagation in axons has a lower limit 0.1 m/s (Segev and Schneidman, 1999). Therefore, the delay between neighboring columns is not larger than milliseconds, which is at least one order of magnitude smaller than the characteristic timescales of the system. We performed simulation with non-zero time delays that reproduced qualitatively the results reported below, thereby validating this approach.

The input component ∑_j p^{j, i}_coup(t) in Equation (6) stands for the excitatory contribution to column i from its neighboring columns, and is given by:

$\begin{array}{cc}\begin{array}{l}{\displaystyle \sum _j{p}_{\text{coup}}^{j,i}(t)}=\frac{1}{N-1}{\displaystyle \sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N{\mathcal{K}}_{j,i}}\text{Sigm}[y_1^j(t)-y_2^j(t)],\hfill \end{array}& \left(11\right)\end{array}$

where the matrix stands for the excitatory connectivity strengths between columns. The diagonal elements of are equal to zero (no self-connectivity), and all other elements are equal to each other, taking a value that from now on will be referred to as the connectivity strength K. The total external contribution to a single cortical column is normalized by the number of its afferent connections. The contribution from one cortical column to the input of other columns is determined by the activity of the former, which depends on its own total input. So, the values of p_const and ∑_j p^{j, i}_coup (through K) are chosen adequately in order to establish self consistent conditions that yield an effective input 〈p^pyr_{ext, exc}〉 ≃ 90 Hz.

2.3. Numerical Methods

The model was integrated using the stochastic Heun integration method (Toral and Colet, 2014) with a time step equal to 0.001 s. We validated the integration method by running computations with decreased integration step. In each run we simulated 1010 s of activity, discarding the first 10 s. We computed the power spectra by applying the Welch algorithm from the Matplotlib Python module, using a Hanning window. The length of each time segment was chosen to be 20 s, with an overlap between segments equal to 10 s.

3. Results

3.1. Spectral and Temporal Properties of a Single Column

We first studied the behavior of a single cortical column receiving an input p_const = 90 Hz, systematically varying the parameters of the Ornstein-Uhlenbeck noise, namely its noise intensity D (varied in a range from 0.1 Hz to 1000 Hz) and correlation time τ (varied in a range from 0.001 s to 10 s). We obtained the power spectrum in each case and compared it with the one reported in the experiments of Bayer et al. (2011) (see Figure 2A). Our goal here was to choose the noise parameters for which the computational result reproduced the experimental characteristics, namely an 1/f^b shape with an embedded peak in the alpha band. Figure 2B shows three power spectra obtained from the model driven by noisy inputs with the same variance but different intensities and correlation times. We found that experimental characteristics were qualitatively reproduced for τ = 0.15 s and noise intensity D = 350 Hz (see Figures 2A,B). Our result is robust for a range of D and τ values, provided that τ ≲ 0.2 s. Beyond that region, the power spectrum at low frequencies (≲ 2 Hz) becomes noticeably steeper than at higher frequencies (green trace in Figure 2B), which is not the case for experimental data. Moreover, for these large τ values the alpha peak becomes too prominent, whereas in the opposite limit it decreases as the correlation time τ is reduced (blue trace in Figure 2B), becoming for τ < 0.15 s significantly smaller (with respect to the 1/f^b background) than in experimental data. This dependence of the signal spectrum on the noise characteristics can be explained in the following way: the system operates on average close to the Hopf bifurcation, where the limit cycle regime begins. This regime is explored transiently by the system due to the stochastic driving. The duration of the episodes in which the system stays in the oscillatory regime is dictated by the correlation time of the noise. Small τ implies rapid changes of the input to the system, which does not have time to relax in the limit cycle regime and alpha oscillations do not occur. In contrast, for relatively large τ, the input changes in a more smooth manner and the system has time to relax and exhibit alpha oscillations, which contribute to the alpha peak in the power spectrum. The noise intensity D plays a role too, because it dictates how deep the system can go into the limit cycle regime. The broadband shape of the power spectrum roughly follows the shape of the spectrum of the noise, which depends on the control parameter τ. This effect is noticeable specially for low frequencies, and originates in the regime which in deterministic conditions corresponds to a fixed point (referred to as a “random fixed point” from now on), where the system follows the noisy driving, and thus yields power spectra similar to that of the Ornstein-Uhlenbeck noise. In this way, the combined effect of different dynamics gives rise to a realistic power spectrum. This effect does not rely on critical behavior, but requires that the system explores different dynamical regimes.

FIGURE 2

Figure 2. Comparison of experimental and computational power spectra. (A) Shows experimental data extracted from Figure 1D of the article by Bayer et al. (2011). (B) Shows power spectra obtained for three different noise parameter sets. Variance of the noise is the same for each case, but noise intensity D and correlation time τ change. Black line corresponds to time series shown in (C), green line to conditions when the correlation time τ is 10 times increased and the blue one when it is 10 times decreased. (C) Shows the time trace obtained from a Jansen-Rit model of a single cortical column with an input consisting of two parts: a constant component equal to 90 Hz and a stochastic one determined by Ornstein-Uhlenbeck noise with intensity D = 350 Hz and correlation time τ = 0.15 s. Three distinct types of dynamics are apparent: time trace begins with a noisy behavior based on the fixed point, then spiky dynamics show up and finally time trace ends with alpha oscillations.

The correlation time τ dictated by synaptic effects is conjectured to be of the order of 10 ms (Mazzoni et al., 2008; Sancristóbal et al., 2013) rather than 100 ms. In our case, however, noise stands for background activity arising from collective effects at the mesoscopic scale. In this study we focused on alpha rhythm, which has a characteristic time scale ~ 0.1 s, therefore choosing a noise correlation time τ equal to 0.15 s is reasonable. Note also that in the computational results shown in panel B of Figure 2 the alpha peak is shifted toward lower frequencies with respect to experimental results shown in panel A. The location of the peak could have been shifted by changing parameters of the model, however we chose to perform the analysis with the original set of parameters proposed by Jansen and Rit (1995), in order to maintain coherence with studies that adopted that set of parameters.

Although the power spectrum obtained with these noise parameters (Figure 2B) reproduces qualitatively the experimental results (Figure 2A), its corresponding temporal evolution exhibits a strong spiky behavior (Figure 2C), which is far from what is typically observed in experimental EEG recordings of healthy subjects. An adequate change of noise parameters would suppress the spiky dynamics in favor of alpha oscillations, but then the alpha peak in the power spectrum would be much higher than the 1/f^b background, in contrast with the experimental observations. Therefore, we conclude that the behavior of a single column is not able to recapitulate realistically both the temporal and spectral characteristics of the experimental observations at the same time. For this reason we extended our model to several coupled columns.

3.2. Coupled Cortical Columns

The signals measured in experimental EEG recordings do not arise from a single cortical column, but from an aggregate of columns. In order to take this into account we extended our model to represent multiple coupled columns. As a simplifying assumption, we consider that the signal registered by an electrode is an average of the signals generated by individual columns in the probed area. The model considers only excitatory connections between populations of pyramidal neurons residing in different cortical columns. There are two new parameters with respect to the single-column case: the number of columns N and the coupling strength K. Our aim was to test our cross-frequency transfer hypothesis in a simple model, therefore we considered only N = 4 cortical columns coupled in a simple all-to-all manner with equally strong connections. Already for 4 coupled columns, individual spikes in the temporal domain are substantially attenuated due to averaging (Figure 3A), rendering time traces that qualitatively resemble EEG signals. On the other hand, the power spectrum of the averaged signal still resembles the experimental one (Figure 3B). This approach finds support in experiments; so-called “microseizures”—spiky, epileptic-like activity that may be detected only in a very fine spatial scale (~1 mm electrode array resolution, 40 μm electrode size)—have been observed experimentally by Stead et al. (2010), not only in epileptic subjects, but also sporadically in healthy ones.

FIGURE 3

Figure 3. Effect of signal averaging. (A, B) Show a typical time trace and power spectrum, respectively, analogous to (B, C) of Figure 2, but obtained from the model of four columns coupled in an all-to-all bidirectional manner. The coupling strength constant K was set to 15. Besides the input from the other columns, each column received the same external input as in the single column case, described above. (A) Shows time traces of two individual columns (dashed lines) and the average signal of all four columns (solid line). Averaging attenuates individual spikes, but does not affect the power spectrum substantially (as shown in B).

The coupling contribution to each column was normalized by the number of afferent columns [N − 1 for our all-to-all connectivity topology, as expressed in Equation (11)]. This normalization allowed us to study the dependence of the results on the system size N for a constant value of the coupling strength K. Without that normalization, increasing the number of columns N in all-to-all topology would ultimately lead to saturation due to excessive external driving. Increasing the number of columns to N = 100 showed that the spectrum remains qualitatively consistent with the experimental one independently of the system size.

We next studied how robust the behavior shown in Figure 3 above is with respect to changes in the constant input and coupling strength. To vary those parameters it is necessary to take into account the fact that the input into a cortical column from other columns [∑jp^{j, i}_coup(t)] is implicitly dependent on p_const, which has the same value for all columns and determines the dynamical regime in which the columns operate. In order to keep the effective total input 〈p^{pyr, i}_{ext, exc}〉 close to 90 Hz, the coupling strength K needs to be compensated by reduction in the input constant component p_const below 90 Hz. Taking this into consideration, we ran a series of simulations for N = 4 columns varying both the constant input component p_const and the coupling strength K. For each condition we averaged the inputs coming from the coupling terms over time and over all columns. In this manner we obtained the mean contribution of inter-column coupling to the input of a column. Adding this value to the constant input p_const gives the average total input acting upon a column, 〈p^pyr_{ext, exc}〉. We varied the coupling strength in the range 0 < K < 70, and in each case we chose the constant input component within the range 50 Hz < p_const < 90 Hz in such a way that the average external input to each column 〈p^pyr_{ext, exc}〉 was close to 90 Hz, in accordance with our assumption regarding proximity to the bifurcation point. In these conditions we found a linear dependence of the average of the coupling input on coupling strength K. Moreover, we found that the variation coefficient of this input was close to unity, which indicates that not only its average, but also its standard deviation grows linearly with K. This can be explained by two effects. First, higher synchronization causes in-sync spiking that weakens the effect of averaging between the columns. Second, these periods of high activity alternate with periods of low activity when due to the lower value of p_const the columns operate effectively in the fixed point regime.

In the whole parameter range that fulfills the condition 〈p^pyr_{ext, exc}〉 ~ 90 Hz, the temporal behavior and power spectrum resembled the experimentally observed one (results not shown). With that caveat, our results are robust not only with respect to the parameters of the Ornstein-Uhlenbeck noise and the system size, but also with respect to the constant input and coupling strengths. The simulations described below correspond to K = 15 and p_const = 75 Hz.

3.3. Effect of an Oscillatory Input

The experimental study that we are interested in was performed by Bayer et al. (2011) on human subjects. Volunteers were placed on a bed, which was swung at a frequency 0.25 Hz. In our model we represent the stimulus associated with this movement as a harmonic driving. EEG data was recorded from the Fz electrode during the N2 sleep phase of the subjects for two conditions: swinging (bed in motion) and stationary (bed still). The experiment showed that swinging facilitates the transition from the awake state to sleep, and that it enhances the EEG power of both slow and alpha oscillations (see Figure 4A). According to the experimental setup the bed motion is harmonic, thus we started by mimicking swinging applying to each column in our model an oscillatory component à sin (2πft), with f = 0.25 Hz. We set the driving amplitude à to 45 Hz, and left all other parameters unchanged with respect to the stationary conditions described in the previous section.

FIGURE 4

Figure 4. Effect of low-frequency driving. (A) Shows experimental EEG power spectra published by Bayer et al. (2011) (Figure 1D of that article), recorded from the Fz electrode in N2 sleep phase for both stationary (black line) and swinging (red line) conditions. Bayer et al. (2011) performed paired 2-tailed t-tests and found statistically significant increase in power in the ranges denoted here with horizontal arrows. In these ranges they found p-value to be p < 0.05, except for the frequency range 1 Hz–2 Hz, where p < 0.005. (B) Shows a typical time trace obtained from the model for two full periods of the driving sine signal, for which (C) shows the corresponding power spectra obtained in the periodically driven case (blue line) compared with the absence of driving (black line) and driving with the composed signal (red). This signal comprises sinusoid ingredients with varying frequencies, amplitudes and random phases. See text for details. Gray marks one standard deviation of distribution of power spectra obtained for different values of phases in the composed signal. In all cases four all-to-all connected columns subject to Ornstein-Uhlenbeck noise with intensity D = 350 Hz and correlation time τ = 0.15 s were used. In the case of sinusoidal driving, the input had amplitude à = 45 Hz and frequency f = 0.25 Hz. The columns were coupled with coupling strength K = 15.

A typical time trace and power spectrum of the signal obtained from the model with harmonic driving are presented in Figures 4B,C. Panel B shows that both the instantaneous amplitude of the alpha oscillations and the average value of the signal are modulated by driving signal. The mechanism underlying both these effects originates in the bifurcation structure of the model (Grimbert and Faugeras, 2006; Spiegler et al., 2010), and was first reported by Tsodyks et al. (1997) for the case of gamma-theta coupling in a Wilson-Cowan model. Recently this effect was more generally discussed by Deco et al. (2008). We now examine in detail this mechanism for our case. In stationary conditions, the system explores four different dynamical regimes: two random fixed points, a regime of spiky behavior and a regime of alpha oscillations. This exploration arises as a consequence of the conjunction of two factors. The first factor is the proximity and coexistence of different dynamical regimes in the vicinity of the chosen input value 〈p^pyr_{ext, exc}〉 = 90 Hz. The second factor is stochastic driving, which enforces alternations between these regimes. For the noise parameters chosen here, D = 350 Hz and τ = 0.15 s, the variance of the Ornstein-Uhlenbeck noise is ≃ 48 Hz, which means that the system may explore two regimes separated by the Hopf bifurcation at 89.83 Hz: a regime of alpha oscillations and a random fixed point. Moreover, the system can enter the spiky dynamics regime, which starts at a saddle-node bifurcation at 113.58 Hz. Below this point the system may enter another random fixed point regime from which it may undergo noisy excitation and also exhibit spiky dynamics.

In the presence of oscillatory driving this situation changes and a number of factors contribute to an overall increase of power in the alpha band. Firstly, the amplitude of alpha oscillations is smaller in the direct vicinity of the Hopf bifurcation (p^pyr_{ext, exc} ≳ 90 Hz), than for greater p^pyr_{ext, exc} values, determined by the driving signal amplitude. Secondly, for p^pyr_{ext, exc} > 137.38 Hz the alpha oscillations become the only allowed dynamics, so flipping between different regimes ceases to occur. This results not only in an increased alpha activity of each individual column, but also in an increase of synchronization between the columns. For sufficiently high driving amplitude columns go through transient in-phase synchronization periods, where the averaged amplitude of their alpha oscillations is greater than in the periods of unsynchronized alpha activity. These periods may occur for only some, or for all columns in the system, they may be terminated by noise and then may reappear due to coupling between the columns. A few peaks, which emerge due to synchronization between the columns during the driving sine upswing are shown in Figure 4B. They appear only for single or few oscillation cycles. It is due to the fact that for this case the driving amplitude is à = 45 Hz and thereby sets the maximal value of the deterministic part of the input to 135 Hz, which coincides with the starting point of the purely alpha regime at 137.38 Hz. This regime may still be explored due to noise, but under these conditions the system operates at best on the edge of resonance (Kang et al., 2010; Battaglia and Hansel, 2011). To the contrary, for higher driving amplitudes synchronization develops fully, enhancing further power increase in the alpha band, as shown in Figure 5A.

FIGURE 5

Figure 5. Relative change of the power spectrum with respect to stationary conditions as the function of the amplitude (A) and frequency (B) of the driving oscillatory signal, and the number of columns in the system (C). Color represents the relative change of the power spectrum expressed in dB, as defined by 10 log₁₀(PSD_driven/PSD_stat), where PSD_driven is the power spectrum in the presence of driving and PSD_stat is the power spectrum in its absence. In (A, B) the number of columns was N = 4, in (B, C) the driving amplitude à was fixed to 45 Hz and in (A, C) the driving frequency was fixed to 0.25 Hz. Again, the analysis was performed on the system of all-to-all connected Jansen-Rit models of cortical columns in the presence of Ornstein-Uhlenbeck noise with intensity D = 350 Hz and correlation time τ = 0.15 s. The coupling strength between the modules K was set to 15.

During the positive half of the driving cycle the alpha oscillations are superimposed to the oscillatory signal (see Figure 4B). During the negative half of the cycle the system moves further away from the oscillatory regimes and it may dwell in one of the random fixed point regimes. Consequently, less spiky behavior is observed, which results in a slight decrease of power in 0 Hz–5 Hz frequency band that corresponds to spiking (Faugeras et al., 2009). This decrease is shown in Figure 4C and in all panels of Figure 5. In the case of alpha oscillations, however, the system may remain in the regime of the random fixed point associated with the alpha limit cycle via the Hopf bifurcation. Noisy perturbations alter the system which may oscillate with alpha frequency. For the noise parameters chosen here (D = 350 Hz and τ = 0.15 s) the latter effect is minor. Nevertheless, all the discussed effects together lead to an overall increase of power in the alpha band in comparison to stationary conditions, as shown in Figure 4C, thereby reproducing qualitatively the experimental observation shown in Figure 4A.

The experiment also showed a statistically significant increase of power for low frequencies. The increase is also observed in the model, but is much more centered (peaked) at the driving frequency (0.25 Hz) than in the experiment, which is much smoother, probably due to reshaping of the low-frequency harmonic signal by sensory, thalamic and/or thalamocortical processing. In order to test this assumption, in the next step we drove the model with a reshaped signal of the form:

$\begin{array}{cc}f(t)=\widetilde{A^{\prime }}{\displaystyle \sum _{n=n_{\text{min}}}^{n_{\text{max}}}1}{0}^{-\frac{n·f_{\text{step}}-f_{\text{min}}}{f_{\text{max}}-f_{\text{min}}}}\sin [2\pi (n{f}_{\text{step}}t+X_n)]& \left(12\right)\end{array}$

where n_min = f_min/f_step, n_max = f_max/f_step, f_min = 0.05 Hz is the minimal frequency of the driving, f_max = 4 Hz is the maximal one, f_step was set to 0.05 Hz and X_n is a random number in the range [0,1). This formula describes a signal composed of a sum of sines with frequencies from f_min to f_max taken every f_step, with exponentially decaying amplitudes and with randomly distributed phases. The choice of f_max has been dictated by the upper limit of the frequency interval in which Bayer et al. (2011) observed a significant increase of power. The amplitude Ã′ = 10.76 Hz was set is such a way that this composed signal delivered the same power to the model as the previously used simple sine signal with amplitude à = 45 Hz. We performed 10 full simulations for different distributions of random phases and averaged the power spectra obtained. The resulting averaged spectrum along with one standard deviation of the power spectrum distribution is shown in Figure 4C. This figure shows that driving with the composed signal reproduces the experimental results better than driving with a simple sine: the increase in the alpha band is present regardless of the randomization of phases, and instead of a decrease of power for low frequencies (as observed in the simple sine driving case), an increase (similar to the experimental result) is observed.

Next, in order to examine systematically the response of this broadband oscillation to harmonic signals, we studied the response of our model for N = 4 columns to a variation in the amplitude and frequency of the driving. Figure 5A depicts the dependence of the power spectrum (in color code) on the driving signal amplitude, with red indicating an increase in power and blue representing a decrease with respect to the stationary conditions. This figure shows that the power increase in the alpha band is robust with respect to the driving amplitude, provided its value is large enough. The slight increase in the frequency that responds maximally, observed for large amplitude values, might be understood from the fact that the frequency of the limit cycle exhibited by the NMM increases slightly for increasing input to the columns (Spiegler et al., 2010).

More importantly, we examined the response of the model with N = 4 columns for a large range of signal frequencies, ranging from values much smaller than its intrinsic alpha frequency (as we have been discussing so far) all the way to much larger frequencies (up to 25 Hz). The results, shown in Figure 5B, reveal that an increase in alpha occurs only for low-enough driving frequencies (f ≲ 2 Hz). As the driving frequency increases, the initial response at alpha splits and leads to increase of power at frequencies smaller and larger than alpha. Interestingly, at this point the alpha band undergoes a decrease, rather than an increase, in power. The low-frequency power (≲ 0.5 Hz) is also reduced for a wide range of driving frequencies. The response is dominated by a straight diagonal line corresponding to 1:1 response to the driving frequency, and by its first harmonic. This strong 1:1 response means that in the case of harmonic driving every injected frequency is transferred by the system. This explains why driving of the form as in Equation (12) leads to a broad increase of power in low frequencies in the system's output. The same study performed for stronger driving showed that the effects discussed above are robust with respect to the driving amplitude (results not shown), although for stronger signals higher harmonics show up and relative changes in the power spectrum are enhanced and widened. Finally, we studied the impact of the system size on the observed effect. The result shown in Figure 5C indicates that under the chosen conditions (oscillatory driving with amplitude à = 45 Hz and frequency f = 0.25 Hz) the results are robust with respect to the system size.

4. Discussion

We have studied a minimal model that gives rise to broadband oscillations in the alpha regime. The model consists of a small number of cortical columns coupled in an all-to-all configuration. Similarly to what happens in microscopic models (Sancristóbal et al., 2013), taking into account that the neuronal population receives a background signal from the rest of the brain in the form of a temporally correlated Ornstein-Uhlenbeck noise (Destexhe and Rudolph, 2004) leads to the characteristic 1/f^b spectral profile of neural activity. Furthermore, selecting adequately the population dynamics of the NMM produces oscillations with a well-defined frequency (David and Friston, 2003), which appear superimposed to the 1/f^b profile. On the other hand, generating this profile entails that the NMM operates in a spiking regime, which differs from the characteristic dynamics observed via EEG. This type of macroscopic measurement, however, reflects the behavior of multiple coupled columns, and we showed that when this situation is considered in our model the spiking behavior disappears due to averaging, rendering signals which recapitulate the experimentally observed EEG while maintaining the broad power spectrum.

Within this dynamical regime, we have examined the effect of a low-frequency driving with a simple sine signal and a composed one, showing a cross-frequency transfer through which this driving signal increases the power not only of low-frequency rhythms, but also of the alpha activity. The result qualitatively reproduces the experimental observations of Bayer et al. (2011) on the effect of rocking on alpha activity and sleep and is robust with respect to the choice of model parameters. The increase of power in the alpha band results from both an enhancement of alpha activity of individual cortical columns and collective synchronization effects. Our results suggest that certain types of cross-frequency transfer in the brain can be simply the result of passive driving of a broadband neuronal oscillator, which brings this effect close to the vast body of work dealing with the driving and synchronization of chaotic oscillators (Pikovsky et al., 2001; Boccaletti et al., 2002; Anishchenko et al., 2007) Interestingly, a systematic analysis shows that the frequency transfer only occurs toward the intrinsic frequency of the oscillator (alpha) when the driving frequency is low; as it increases the response shifts to both lower and higher frequencies, and the power in the alpha band decreases instead of increasing. Taken together, our results indicate that a relatively simple oscillation generation mechanism in neuronal populations has a strongly nontrivial response to periodic driving, providing a rich scenario to interpret a variety of cross-frequency phenomena in the brain.

Conflict of Interest Statement

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

Acknowledgments

We thank the reviewers of the paper for very insightful comments. This work was supported by the European Commission through the FP7 Marie Curie Initial Training Network 289146 (NETT: Neural Engineering Transformative Technologies), and the Ministerio de Economia y Competividad (Spain, project FIS2012-37655). JGO acknowledges support from the ICREA Academia programme.

References

Aburn, M. J., Holmes, C. A., Roberts, J. A., Boonstra, T. W., and Breakspear, M. (2012). Critical fluctuations in cortical models near instability. Front. Physiol. 3:331. doi: 10.3389/fphys.2012.00331

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Amari, S. (1974). Homogeneous nets of neuron-like elements. Biol. Cybern. 17, 211–220. doi: 10.1007/BF00339367

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Anishchenko, V. S., Astakhov, V., Vadivasova, T., Neiman, A., and Schimansky-Geier, L. (2007). Nonlinear Dynamics of Chaotic and Stochastic Systems: Tutorial and Modern Developments. Berlin: Springer.

Google Scholar

Aru, J., Aru, J., Priesemann, V., Wibral, M., Lana, L., Pipa, G., et al. (2015). Untangling cross-frequency coupling in neuroscience. Curr. Opin. Neurobiol. 31C, 51–61. doi: 10.1016/j.conb.2014.08.002

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Battaglia, D., and Hansel, D. (2011). Synchronous chaos and broad band gamma rhythm in a minimal multi-layer model of primary visual cortex. PLoS Comput. Biol. 7:e1002176. doi: 10.1371/journal.pcbi.1002176

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bayer, L., Constantinescu, I., Perrig, S., Vienne, J., Vidal, P.-P., Mühlethaler, M., et al. (2011). Rocking synchronizes brain waves during a short nap. Curr. Biol. 21, R461–R462. doi: 10.1016/j.cub.2011.05.012

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bédard, C., Kröger, H., and Destexhe, A. (2006). Does the 1/f frequency scaling of brain signals reflect self-organized critical states? Phys. Rev. Lett. 97:118102. doi: 10.1103/PhysRevLett.97.118102

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Boccaletti, S., Kurths, J., Osipov, G., Valladares, D. L., and Zhou, C. S. (2002). The synchronization of chaotic systems. Phys. Rep. 366, 1–101. doi: 10.1016/S0370-1573(02)00137-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bojak, I., and Liley, D. T. J. (2010). Axonal velocity distributions in neural field equations. PLoS Comput. Biol. 6:e1000653. doi: 10.1371/journal.pcbi.1000653

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bruns, A., and Eckhorn, R. (2004). Task-related coupling from high- to low-frequency signals among visual cortical areas in human subdural recordings. Int. J. Psychophysiol. 51, 97–116. doi: 10.1016/j.ijpsycho.2003.07.001

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Buzsáki, G., and Draguhn, A. (2004). Neuronal oscillations in cortical networks. Science 304, 1926–1929. doi: 10.1126/science.1099745

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Canolty, R. T., and Knight, R. T. (2010). The functional role of cross-frequency coupling. Trends Cogn. Sci. 14, 506–515. doi: 10.1016/j.tics.2010.09.001

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Chialvo, D. R. (2010). Emergent complex neural dynamics. Nat. Phys. 6, 744–750. doi: 10.1038/nphys1803

CrossRef Full Text | Google Scholar

Cohen, M. X., Axmacher, N., Lenartz, D., Elger, C. E., Sturm, V., and Schlaepfer, T. E. (2009a). Good vibrations: cross-frequency coupling in the human nucleus accumbens during reward processing. J. Cogn. Neurosci. 21, 875–889. doi: 10.1162/jocn.2009.21062

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Cohen, M. X., Elger, C. E., and Fell, J. (2009b). Oscillatory activity and phase-amplitude coupling in the human medial frontal cortex during decision making. J. Cogn. Neurosci. 21, 390–402. doi: 10.1162/jocn.2008.21020

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Cona, F., Zavaglia, M., Massimini, M., Rosanova, M., and Ursino, M. (2011). A neural mass model of interconnected regions simulates rhythm propagation observed via TMS-EEG. Neuroimage 57, 1045–1058. doi: 10.1016/j.neuroimage.2011.05.007

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Coombes, S. (2010). Large-scale neural dynamics: simple and complex. Neuroimage 52, 731–739. doi: 10.1016/j.neuroimage.2010.01.045

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Cox, R., van Driel, J., de Boer, M., and Talamini, L. M. (2014). Slow oscillations during sleep coordinate interregional communication in cortical networks. J. Neurosci. 34, 16890–16901. doi: 10.1523/JNEUROSCI.1953-14.2014

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

David, O., and Friston, K. J. (2003). A neural mass model for MEG/EEG: coupling and neuronal dynamics. Neuroimage 20, 1743–1755. doi: 10.1016/j.neuroimage.2003.07.015

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

de Lange, F. P., Jensen, O., Bauer, M., and Toni, I. (2008). Interactions between posterior gamma and frontal alpha/beta oscillations during imagined actions. Front. Hum. Neurosci. 2:7. doi: 10.3389/neuro.09.007.2008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Deco, G., Jirsa, V. K., and McIntosh, A. R. (2013). Resting brains never rest: computational insights into potential cognitive architectures. Trends Neurosci. 36, 268–274. doi: 10.1016/j.tins.2013.03.001

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Deco, G., Jirsa, V. K., Robinson, P. A., Breakspear, M., and Friston, K. (2008). The dynamic brain: from spiking neurons to neural masses and cortical fields. PLoS Comput. Biol. 4:e1000092. doi: 10.1371/journal.pcbi.1000092

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Destexhe, A., and Rudolph, M. (2004). Extracting information from the power spectrum of synaptic noise. J. Comput. Neurosci. 17, 327–345. doi: 10.1023/B:JCNS.0000044875.90630.88

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Erlhagen, W., and Bicho, E. (2006). The dynamic neural field approach to cognitive robotics. J. Neural Eng. 3, R36. doi: 10.1088/1741-2560/3/3/R02

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Faisal, A. A., Selen, L. P. J., and Wolpert, D. M. (2008). Noise in the nervous system. Nat. Rev. Neurosci. 9, 292–303. doi: 10.1038/nrn2258

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Faugeras, O., Touboul, J., and Cessac, B. (2009). A constructive mean-field analysis of multi-population neural networks with random synaptic weights and stochastic inputs. Front. Comput. Neurosci. 3:1. doi: 10.3389/neuro.10.001.2009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Freeman, W. (1972). Linear analysis of the dynamics of neural masses. Annu. Rev. Biophys. Bioeng. 1, 225–256. doi: 10.1146/annurev.bb.01.060172.001301

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Freeman, W. J., Rogers, L. J., Holmes, M. D., and Silbergeld, D. L. (2000). Spatial spectral analysis of human electrocorticograms including the alpha and gamma bands. J. Neurosci. Methods 95, 111–121. doi: 10.1016/S0165-0270(99)00160-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Freestone, D. R., Nesic, D., Jafarian, A., Cook, M. J., and Grayden, D. B. (2013). “A neural mass model of spontaneous burst suppression and epileptic seizures,” in Conference proceedings: 2013, 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 5942–5945.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Friston, K. J. (2000). The labile brain. I. Neuronal transients and nonlinear coupling. Philos. Trans. R. Soc. Lond. B Biol. Sci. 355, 215–236. doi: 10.1098/rstb.2000.0560

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Goodfellow, M., Schindler, K., and Baier, G. (2012). Self-organised transients in a neural mass model of epileptogenic tissue dynamics. Neuroimage 59, 2644–2660. doi: 10.1016/j.neuroimage.2011.08.060

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Grimbert, F., and Faugeras, O. (2006). Bifurcation analysis of Jansen's neural mass model. Neural Comput. 18, 3052–3068. doi: 10.1162/neco.2006.18.12.3052

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Helmstaedter, M., de Kock, C. P. J., Feldmeyer, D., Bruno, R. M., and Sakmann, B. (2007). Reconstruction of an average cortical column in silico. Brain Res. Rev. 55, 193–203. doi: 10.1016/j.brainresrev.2007.07.011

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hubel, D. H., and Wiesel, T. N. (1977). Functional architecture of macaque monkey visual cortex. Proc. R. Soc. Lond. B Biol. Sci. 198, 1–59. doi: 10.1098/rspb.1977.0085

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hutt, A., and Frank, T. (2005). Critical fluctuations and 1/fα-activity of neural fields involving transmission delays. Acta Phys. Pol. A 108, 1021–1040.

Google Scholar

Hutt, A., and Atay, F. (2006). Effects of distributed transmission speeds on propagating activity in neural populations. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 73:021906. doi: 10.1103/PhysRevE.73.021906

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Isler, J. R., Grieve, P. G., Czernochowski, D., Stark, R. I., and Friedman, D. (2008). Cross-frequency phase coupling of brain rhythms during the orienting response. Brain Res. 1232, 163–172. doi: 10.1016/j.brainres.2008.07.030

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Jansen, B. H., and Rit, V. G. (1995). Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol. Cybern. 73, 357–366. doi: 10.1007/BF00199471

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Jansen, B. H., Zouridakis, G., and Brandt, M. E. (1993). A neurophysiologically-based mathematical model of flash visual evoked potentials. Biol. Cybern. 68, 275–283. doi: 10.1007/BF00224863

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Jensen, O., and Colgin, L. L. (2007). Cross-frequency coupling between neuronal oscillations. Trends Cogn. Sci. 11, 267–269. doi: 10.1016/j.tics.2007.05.003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Jirsa, V., and Haken, H. (1996). Field theory of electromagnetic brain activity. Phys. Rev. Lett. 77, 960–963. doi: 10.1103/PhysRevLett.77.960

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Jirsa, V., and Müller, V. (2013). Cross-frequency coupling in real and virtual brain networks. Front. Comput. Neurosci. 7:78. doi: 10.3389/fncom.2013.00078

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Kang, K., Shelley, M., Henrie, J. A., and Shapley, R. (2010). LFP spectral peaks in V1 cortex: network resonance and cortico-cortical feedback. J. Comput. Neurosci. 29, 495–507. doi: 10.1007/s10827-009-0190-2

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Lopes da Silva, F. H., Hoeks, A., Smits, H., and Zetterberg, L. H. (1974). Model of brain rhythmic activity. Kybernetik 15, 27–37. doi: 10.1007/BF00270757

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Marshall, L., Helgadóttir, H., Mölle, M., and Born, J. (2006). Boosting slow oscillations during sleep potentiates memory. Nature 444, 610–613. doi: 10.1038/nature05278

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Massimini, M., Ferrarelli, F., Esser, S. K., Riedner, B. A., Huber, R., Murphy, M., et al. (2007). Triggering sleep slow waves by transcranial magnetic stimulation. Proc. Natl. Acad. Sci. U.S.A. 104, 8496–8501. doi: 10.1073/pnas.0702495104

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Mazzoni, A., Brunel, N., Cavallari, S., Logothetis, N. K., and Panzeri, S. (2011). Cortical dynamics during naturalistic sensory stimulations: experiments and models. J. Physiol. Paris 105, 2–15. doi: 10.1016/j.jphysparis.2011.07.014

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Mazzoni, A., Panzeri, S., Logothetis, N. K., and Brunel, N. (2008). Encoding of naturalistic stimuli by local field potential spectra in networks of excitatory and inhibitory neurons. PLoS Comput. Biol. 4:e1000239. doi: 10.1371/journal.pcbi.1000239

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Mazzoni, A., Whittingstall, K., Brunel, N., Logothetis, N. K., and Panzeri, S. (2010). Understanding the relationships between spike rate and delta/gamma frequency bands of LFPs and EEGs using a local cortical network model. Neuroimage 52, 956–972. doi: 10.1016/j.neuroimage.2009.12.040

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

McDonnell, M., and Ward, L. (2011). The benefits of noise in neural systems: bridging theory and experiment. Nat. Rev. Neurosci. 12, 415–425. doi: 10.1038/nrn3061

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Nguyen Trong, M., Bojak, I., and Knösche, T. R. (2012). Associating spontaneous with evoked activity in a neural mass model of visual cortex. Neuroimage 66C, 80–87. doi: 10.1016/j.neuroimage.2012.10.024

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Nunez, P. (1974). The brain wave equation: a model for the EEG. Math. Biosci. 291, 279–297. doi: 10.1016/0025-5564(74)90020-0

CrossRef Full Text | Google Scholar

Osipova, D., Hermes, D., and Jensen, O. (2008). Gamma power is phase-locked to posterior alpha activity. PLoS ONE 3:e3990. doi: 10.1371/journal.pone.0003990

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Pikovsky, A., Rosenblum, M., and Kurths, J. (2001). Synchronization: A Universal Concept in Nonlinear Sciences. Cambridge, UK: Cambridge University Press. doi: 10.1017/CBO9780511755743

CrossRef Full Text

Pons, A. J., Cantero, J. L., Atienza, M., and Garcia-Ojalvo, J. (2010). Relating structural and functional anomalous connectivity in the aging brain via neural mass modeling. Neuroimage 52, 848–861. doi: 10.1016/j.neuroimage.2009.12.105

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Ponten, S. C., Daffertshofer, A., Hillebrand, A., and Stam, C. J. (2010). The relationship between structural and functional connectivity: graph theoretical analysis of an EEG neural mass model. Neuroimage 52, 985–994. doi: 10.1016/j.neuroimage.2009.10.049

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Sancristóbal, B., Vicente, R., Sancho, J. M., and Garcia-Ojalvo, J. (2013). Emergent bimodal firing patterns implement different encoding strategies during gamma-band oscillations. Front. Comput. Neurosci. 7:18. doi: 10.3389/fncom.2013.00018

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Segev, I., and Schneidman, E. (1999). Axons as computing devices: basic insights gained from models. J. Physiol. Paris 93, 263–270. doi: 10.1016/S0928-4257(00)80055-8

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Shirvalkar, P. R., Rapp, P. R., and Shapiro, M. L. (2010). Bidirectional changes to hippocampal theta-gamma comodulation predict memory for recent spatial episodes. Proc. Natl. Acad. Sci. U.S.A. 107, 7054–7059. doi: 10.1073/pnas.0911184107

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Spiegler, A., Kiebel, S. J., Atay, F. M., and Knösche, T. R. (2010). Bifurcation analysis of neural mass models: impact of extrinsic inputs and dendritic time constants. Neuroimage 52, 1041–1058. doi: 10.1016/j.neuroimage.2009.12.081

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Spiegler, A., Knösche, T. R., Schwab, K., Haueisen, J., and Atay, F. M. (2011). Modeling brain resonance phenomena using a neural mass model. PLoS Comput. Biol. 7:e1002298. doi: 10.1371/journal.pcbi.1002298

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Stead, M., Bower, M., Brinkmann, B. H., Lee, K., Marsh, W. R., Meyer, F. B., et al. (2010). Microseizures and the spatiotemporal scales of human partial epilepsy. Brain 133, 2789–2797. doi: 10.1093/brain/awq190

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Toral, R., and Colet, P. (2014). Stochastic Numerical Methods. Weinheim: Wiley-VCH. doi: 10.1002/9783527683147

CrossRef Full Text

Tort, A. B. L., Komorowski, R. W., Manns, J. R., Kopell, N. J., and Eichenbaum, H. (2009). Theta-gamma coupling increases during the learning of item-context associations. Proc. Natl. Acad. Sci. U.S.A. 106, 20942–20947. doi: 10.1073/pnas.0911331106

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Ts'o, D. Y., Zarella, M., and Burkitt, G. (2009). Whither the hypercolumn? J. Physiol. 587(Pt 12), 2791–2805. doi: 10.1113/jphysiol.2009.171082

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Tsodyks, M. V., Skaggs, W. E., Sejnowski, T. J., and McNaughton, B. L. (1997). Paradoxical effects of external modulation of inhibitory interneurons. J. Neurosci. 17, 4382–4388.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Ursino, M., Cona, F., and Zavaglia, M. (2010). The generation of rhythms within a cortical region: analysis of a neural mass model. Neuroimage 52, 1080–1094. doi: 10.1016/j.neuroimage.2009.12.084

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Voytek, B., Canolty, R. T., Shestyuk, A., Crone, N. E., Parvizi, J., and Knight, R. T. (2010). Shifts in gamma phase-amplitude coupling frequency from theta to alpha over posterior cortex during visual tasks. Front. Hum. Neurosci. 4:191. doi: 10.3389/fnhum.2010.00191

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wang, P., and Knösche, T. R. (2013). A realistic neural mass model of the cortex with laminar-specific connections and synaptic plasticity - evaluation with auditory habituation. PLoS ONE 8:e77876. doi: 10.1371/journal.pone.0077876

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wendling, F., Bellanger, J. J., Bartolomei, F., and Chauvel, P. (2000). Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals. Biol. Cybern. 83, 367–378. doi: 10.1007/s004220000160

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wilson, H. R., and Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1–24. doi: 10.1016/S0006-3495(72)86068-5

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar