Critical Fluctuations in Cortical Models Near Instability

Computational studies often proceed from the premise that cortical dynamics operate in a linearly stable domain, where fluctuations dissipate quickly and show only short memory. Studies of human electroencephalography (EEG), however, have shown significant autocorrelation at time lags on the scale of minutes, indicating the need to consider regimes where non-linearities influence the dynamics. Statistical properties such as increased autocorrelation length, increased variance, power law scaling, and bistable switching have been suggested as generic indicators of the approach to bifurcation in non-linear dynamical systems. We study temporal fluctuations in a widely-employed computational model (the Jansen–Rit model) of cortical activity, examining the statistical signatures that accompany bifurcations. Approaching supercritical Hopf bifurcations through tuning of the background excitatory input, we find a dramatic increase in the autocorrelation length that depends sensitively on the direction in phase space of the input fluctuations and hence on which neuronal subpopulation is stochastically perturbed. Similar dependence on the input direction is found in the distribution of fluctuation size and duration, which show power law scaling that extends over four orders of magnitude at the Hopf bifurcation. We conjecture that the alignment in phase space between the input noise vector and the center manifold of the Hopf bifurcation is directly linked to these changes. These results are consistent with the possibility of statistical indicators of linear instability being detectable in real EEG time series. However, even in a simple cortical model, we find that these indicators may not necessarily be visible even when bifurcations are present because their expression can depend sensitively on the neuronal pathway of incoming fluctuations.

Introduction

Computational models of neocortex and other brain structures have proved very useful for a range of research problems in neuroscience (Braun and Mattia, 2010; Friston and Dolan, 2010). Interpreting empirical data using dynamical models is particularly fruitful in neuroimaging, where underlying processes are obscured by the low temporal resolution of fMRI or the coarse spatial source resolution of electroencephalography (EEG)/MEG. This allows testing of hypotheses about internal dynamical mechanisms (e.g., Freyer et al., 2012) and, through model inversion, the estimation of neural and connectivity parameters that cannot be observed directly (Friston et al., 2003). In contrast to modeling at the microscopic scale, where the range of dynamics of healthy neurons is known to include non-linear behavior such as limit cycles, modeling at the larger scale of mesoscopic neural masses, or neural fields often assumes that the dynamics at this scale operate close to a stable fixed point where input fluctuations result in only small and brief perturbations of the population state. This premise is predicated on the diffusion approximation that states that correlations amongst neuronal inputs are reduced as the size of the population increases (for review, see Deco et al., 2008). This approach enables the calculation of spectra from the composition of transfer functions, a powerful technique that allows physiological parameters to be estimated from non-invasive functional neuroimaging (Friston et al., 2003) and neurophysiological (van Albada et al., 2010) data.

Dynamic instabilities in models at the larger scale of neural masses have typically been associated with the pathological activity of epileptic seizures (Wendling et al., 2000; Robinson et al., 2002; Breakspear et al., 2006). However, empirical data shows that such instabilities may also underlie healthy neural activity (Freyer et al., 2009, 2011). Indeed, the Jansen–Rit neural mass model (Jansen and Rit, 1995) and its derivatives (Wendling et al., 2002; David and Friston, 2003; Zavaglia et al., 2006; Moran et al., 2007; Sotero et al., 2007; Spiegler et al., 2010) reach bifurcations where fixed points become linearly unstable while still within the healthy physiological range of parameters. In fact, oscillations in the model output that have been identified with normal cortical alpha activity have been shown to arise from limit cycle activity following a supercritical Hopf bifurcation (Grimbert and Faugeras, 2006; Spiegler et al., 2010).

The term “linear instability” here does not necessarily imply that the dynamics of the system as a whole lose stability. Indeed, in the case of the supercritical Hopf bifurcation, stability of the attractor is maintained as it deforms continuously from a stable fixed point to a stable limit cycle, which then increases in size in the phase space. Hence, there is no discontinuous transition. The distinction is that the dominant dynamics in the system are no longer linear. The presence of quadratic and higher order flow terms that become significant in the neighborhood of a bifurcating fixed point have a profound influence on the system’s statistical properties and its response to stochastic perturbations.

The putative presence of linear instabilities in healthy, mesoscopic cortical activity is ultimately an empirical question that must be answered with reference to the theory of non-linear stochastic dynamical systems. For a wide range of systems, statistical measures such as increased autocorrelation, increased variance, and bistable switching have been proposed as generic indicators that the system is losing linear stability on approaching a bifurcation (Scheffer et al., 2009; Kelso, 2010). Increased autocorrelation length is a direct consequence of critical slowing down, which occurs as the strength of attraction to a stable fixed point becomes weaker before changing to equally weak repulsion. Long-range correlations may also reveal a transition from exponential to power law relaxation in the vicinity of linear instabilities as a result of the higher order (non-linear) flow terms.

Within neuroscience, statistical indicators of bifurcations have been studied at a range of scales, in both computational models and empirical analyses. In the context of single neuron models, increase of variance close to a bifurcation and the spectral peak near a Hopf bifurcation have been examined (Steyn-Ross et al., 2006). Spectral features and variance close to instability have been explored in large-scale mean field corticothalamic models (Robinson et al., 1997, 2002; Roberts and Robinson, 2012) and mean field models of the brainstem and hypothalamus (Robinson et al., 2010). Slowing down, instability, and bifurcations have also been studied at the highest level of brain function, particularly in human movement. For example, increased variance and critical slowing have been observed in human bimanual motor control (Kelso et al., 1986; Scholz et al., 1987) and are explained by a simplified phenomenological model of coordination (Haken et al., 1985).

In addition to the analyses of empirical data contained within these computational studies, signatures of transitions in neuroimaging data have been the subject of a number of predominantly empirical studies. Amplitude fluctuations of human brain oscillations have been shown to have long time autocorrelations with power law decay in EEG (Linkenkaer-Hansen et al., 2001), consistent with effects expected near linear instability. Scale-free cortical activity has also been reported in surface electrocorticogram (ECoG) activity, although the significance, scaling coefficient, and likely mechanisms remain contested (Bedard et al., 2006; Miller et al., 2009; He et al., 2010). Similarly, Stam and de Bruin (2004) reported scale-free fluctuations in the degree of synchronization between surface EEG recordings. These findings are consistent with prior reports of intermittent non-linear structure within (Stam et al., 1999) and between (Breakspear and Terry, 2002) surface EEG channels. More recently, Freyer et al. (2009) found that 10 Hz oscillations showed intermittent switching between two distinct bistable modes, although the dwell times within each mode followed a stretched exponential, not a power law decay.

The objectives of the present study are to examine linear instabilities in the Jansen–Rit model, a closed set of equations describing the activity of a small cortical region and one of the simplest cortical neural mass models. At the same time it is a base upon which many extensions and derivative models have been built (Wendling et al., 2002; David and Friston, 2003; Zavaglia et al., 2006; Moran et al., 2007; Sotero et al., 2007; Spiegler et al., 2010). The phenomena which we report in this simple model therefore highlight the possibility of similar behavior in a wider class of models. We focus on one key indicator of linear instability (autocorrelation length) and one important bifurcation (supercritical Hopf). Time series for each neural population in the model are generated for sets of parameters approaching a bifurcation. We then test whether the autocorrelation indicator of proximity to bifurcation is reliably detectable in the time series of the pyramidal population and also examine scaling properties of fluctuations in this time series. In this way we explore whether simple bifurcations at the population scale have the potential to contribute to indicators such as lengthened autocorrelation times and power law scaling of fluctuations reported in human EEG data.

Materials and Methods

Jansen–Rit Neural Mass Model

Building on the earlier work of Lopes da Silva et al. (1974) and Wilson and Cowan (1972), Jansen and Rit developed a simple computational model of a small cortical region (Jansen et al., 1993; Jansen and Rit, 1995). The model produces an output signal similar to spontaneous EEG alpha oscillations, and also shows responses similar to evoked potentials following pulsatile input. The Jansen–Rit model is a closed set of differential equations that describe the local average states of three interconnected neural populations (Figure 1), excitatory interneurons, pyramidal cells, and inhibitory interneurons. Here we follow David and Friston (2003) in identifying the excitatory interneurons in the model with layer IV spiny stellate cells. The spiny stellate and pyramidal neurons are both excitatory and both populations receive external input, although only the pyramidal cells project out of the local region.

FIGURE 1

Figure 1. Schematic connectivity of the Jansen–Rit model. (A) Basic connectivity diagram showing the three neuronal populations, their excitatory (arrows) and inhibitory (circle) connections, and inputs from outside the local cortical region (u and p). (B) A block diagram then summarizes how this translates directly to a mathematical model: linear filter boxes labeled h_e(t) and h_i(t) model the mean response of excitatory and inhibitory synapse populations respectively, including their postsynaptic dendritic filtering. Sigmoid boxes (denoted S) represent conversion of mean summed soma membrane potential to mean output firing rate. Connectivity constants γ₁ to γ₄ model the number and strength of connections between populations.

Each second order equation in the model corresponds to a population of synapses and their postsynaptic dendritic processes (Freeman, 1992; Deco et al., 2008). Critically damped second order linear filters describe the time course of the population mean of postsynaptic potentials, further dispersed due to variability of parameters within each population. This mean behavior summarizes both synaptic and dendritic dynamics of the individual neurons. For excitatory and inhibitory synapses respectively these filters are expressed by the differential operators

$L_{\text{e}}=\frac{1}{H_{\text{e}}{\kappa }_{\text{e}}}\left(\frac{d^2}{d{t}^2}+2{\kappa }_{\text{e}}\frac{d}{dt}+{\kappa }_{\text{e}}^2\right),$ (1)

$L_{\text{i}}=\frac{1}{H_{\text{i}}{\kappa }_{\text{i}}}\left(\frac{d^2}{d{t}^2}+2{\kappa }_{\text{i}}\frac{d}{dt}+{\kappa }_{\text{i}}^2\right),$ (2)

where the scalar parameters H_e and H_i determine the maximum amplitude of the postsynaptic population response to excitatory and inhibitory inputs, respectively. The rate constants κ_e and κ_i determine the time scale of these population responses. As the synaptic filters are linear, synapses with different source or target populations can be merged where synapses are assumed to have the same aggregate properties. For example Jansen and Rit (1995) consolidated their original model to just three second order equations, as it was implicitly assumed that excitatory and inhibitory interneuron populations would always have identical states up to a scaling constant. Following synaptodendritic filtering, fluctuations in membrane potential sum in the cell soma and lead to changes in the average population firing rate. The sigmoid function

$S\left(v\right)=\frac{2e_0}{1+exp\left[{\rho }_1\left({\rho }_2-v\right)\right]},$ (3)

describes how the mean firing rate of a neural population depends on the mean soma membrane potential v, incorporating the dispersion of responses due to variability in the parameters and underlying neuronal states (Marreiros et al., 2008). Parameters e₀, ρ₂, and ρ₁ determine the maximum firing rate, threshold potential, and sensitivity, respectively.

We express the Jansen–Rit model as a set of four second order differential equations, thus allowing both pyramidal and spiny stellate populations separately to receive extrinsic input. We follow the variable and parameter names of Moran et al. (2007). The dynamical variables v₁, v₂, and v₄ represent the positive contributions to population mean soma potentials by excitatory synapses targeting spiny stellate, pyramidal, and inhibitory interneuron populations, respectively. Variable v₃ represents the negative contribution to the mean soma potential of the pyramidal population originating from inhibitory synapses. Thus the resulting mean soma potential of the pyramidal population is v₂ − v₃. This is taken as the main output of the model (Jansen and Rit, 1995; David et al., 2005) because the size and orientation of the apical dendrites of pyramidal neurons mean that pyramidal activity is most closely associated with EEG signals. These equations are given by

$\begin{array}{rcll}{L}_e{v}_1& =& {\gamma }_{1}S\left(v_2-v_3\right)+⟨u⟩+{\sigma }_u{\xi }_u\left(t\right),& \text{(4)}\\ L_e{v}_2& =& {\gamma }_{2}S\left(v_1\right)+⟨p⟩+{\sigma }_p{\xi }_p\left(t\right),& \text{(5)}\\ L_i{v}_3& =& {\gamma }_{4}S\left(v_4\right),& \text{(6)}\\ L_e{v}_4& =& {\gamma }_{3}S\left(v_2-v_3\right).& \text{(7)}\end{array}$

Equation 4 describes excitatory synaptic input targeting the spiny stellate population. Eqs 5 and 6 describe excitatory and inhibitory synaptic input to the pyramidal population, respectively. Equation 7 describes excitatory synaptic input to the inhibitory interneuron population. Parameters 〈u〉 and 〈p〉 are the mean per neuron external input firing rates to the cortical region, targeting spiny stellate, and pyramidal populations, respectively. Langevin white noise terms ξ_u(t) and ξ_p(t) in the extrinsic input represent the fluctuations in the input firing rates, with σ_u and σ_p denoting their standard deviations. Scalar connectivity constants γ₁ to γ₄ represent at the population scale the number and strength of connections between the three neural populations.

This system of equations is equivalent to a single eight-dimensional stochastic first order differential system:

$dv=f(v)dt+GdW(t),$ (8)

where matrix elements of G determine the cross-correlation of noise inputs to the pyramidal and spiny stellate populations. This is the equation that we integrate numerically.

Table 1 lists the values of parameters used for all simulations in this study; they are the standard parameter values introduced by Jansen and Rit (1995).

TABLE 1

Table 1. Jansen–Rit standard parameter values.

Jansen and Rit (1995) themselves focused on numerical simulations of this non-linear model. Through a survey of the simulated behavior with physiologically realistic parameters, they observed a variety of noise-driven rhythmic behaviors consistent with human alpha and beta rhythms. Wendling et al. (2000) studied the emergence of “spike-wave” oscillations resembling epileptic activity when the ratio of excitation to inhibition was increased. Bifurcations in this model where subsequently examined by Grimbert and Faugeras (2006) who treated the input p as the bifurcation parameter in order to understand better the original simulation results of Jansen and Rit (limit cycle beyond a Hopf bifurcation causing alpha oscillations) and Wendling et al. (the emergence of a large amplitude non-harmonic oscillator near a sniper bifurcation). More recently, Spiegler et al. (2010) performed a more general bifurcation analysis that included time scale parameters and analyzed the presence of qualitatively different oscillatory regimes.

Bifurcation Parameters

In the original model (Jansen et al., 1993), both pyramidal and excitatory interneuron populations were the targets of extrinsic inputs, with the two inputs being always proportional (fully correlated). In the model of Jansen and Rit (1995), all extrinsic input was delivered to the pyramidal neurons only, with external stimulation of the other population dropped.

David and Friston (2003) revisited the Jansen–Rit model, in particular explicitly identifying the “excitatory interneuron” population of the original model with spiny stellate cells in layer IV of the neocortex. Their motivation was to send extrinsic input to the layer IV spiny stellate cells in the model rather than to the pyramidal cells. This was arguably a more realistic model of connectivity for input representing thalamocortical sensory afferents. However the equations as published retained the pyramidal-only input of the original Jansen–Rit model.

Moran et al. (2007), in the context of Dynamic Causal Modeling (DCM, a framework for model selection and parameter estimation), extended the Jansen–Rit model with several innovations, including firing rate adaptation, recurrent inhibition, and a differently shaped sigmoid function. In particular Moran et al. did change the target of the extrinsic input to be the spiny stellate population, as foreshadowed by David and Friston (2003). We refer to this model as the Moran–Friston model hereafter.

For the present study we minimally extend the Jansen–Rit model, so that extrinsic input can be delivered either to the pyramidal population (as in Jansen and Rit, 1995), the spiny stellate excitatory population (as in Moran et al., 2007) or more realistically a combination of the two. In this way, the system input can be varied continuously from the Jansen–Rit design to the Moran–Friston design or anywhere in between. In addition, for the case of input to both populations, these two inputs can be chosen as uncorrelated, fully correlated, or partially correlated in their fluctuations. Hence we study the bifurcations of this model as input is varied in the combined (u, p) plane. This subset of parameter space includes a one-dimensional space explored by Jansen and Rit containing a supercritical Hopf bifurcation studied by Grimbert and Faugeras (2006) that is within the physiological range of parameters. We present the bifurcation analysis in Section “Bifurcation diagram.”

Numerical Simulation and Analysis

The model is a system of stochastic differential equations (SDEs) with additive noise. Equation 8 is integrated numerically using the Heun algorithm, which is applicable to SDEs in Stratonovich form (Rümelin, 1982). This ensures that noise amplitude is scaled in appropriate proportion to the square root of the integration time step. We use an integration time step of 0.2 ms. The first transient 5 s of each simulation is discarded from further analysis.

As reviewed in the introduction, the approach to linear instabilities in systems of equations such as Eqs 4–7 is widely assumed to cause changes in the autocorrelation length and/or a peak in the power spectral density function (in the case of a Hopf bifurcation). This is because it is often assumed that the linear treatment of these systems – which predicts both an enhancement of spectral peaks and a lengthening of the autocorrelation time – can be extrapolated from the setting when the system is linearly stable to when it is in the neighborhood of a bifurcation.

To estimate the normalized autocorrelation function of the resulting time series we first normalize each time series to a mean of 0 and standard deviation of 1, and then compute the cross-correlation of the series with itself applying unbiased normalization,

${\widehat{R}}_{\text{xx,unbiased}}\left(m\right)=\frac{1}{N-m}\sum _{n=0}^{N-m-1}{x}_{n+m}{x}_{\text{n}}\left(m\ge 0\right),$ (9)

where m is the lag expressed as number of samples (Orfanidis, 1996). In each case we compute autocorrelation at lag times from 0 to 1/4 of the total time series length for further analysis. Since time series are generated in the vicinity of Hopf bifurcations with natural frequency about 11 Hz, the autocorrelation functions all have a strong 11 Hz component. Because we are primarily interested in the decay of the amplitude of this autocorrelation over a longer time scale, each simulation is repeated 16 times with identical parameters to generate 16 sample paths each of duration 600 or 1,800 s. The autocorrelation function is calculated as described above for each sample path separately. The decay is then quantified by calculating the modulus of the Hilbert transform of the normalized autocorrelation functions computed above. The pointwise mean and standard deviation of this autocorrelation amplitude across 16 sample paths are then plotted. Power spectra are estimated using the Welch algorithm with Hamming window and a segment length of 80,000 samples or 16 s.

The full MATLAB code implementing the model, integration, and time series analyses is available from the authors on request.

Results

Bifurcation Diagram

From the earlier bifurcation analysis of Grimbert and Faugeras (2006) the model is known to have a supercritical Hopf bifurcation when the pyramidal input p = 89.8 s⁻¹ and the other parameters are set to the values used by Jansen and Rit (1995). This assumed no input to the excitatory (spiny stellate) population. We label this bifurcation point H1; it has mean input firing rate 〈p〉 = 89.8 s⁻¹ to the pyramidal population and zero input to the spiny stellate population (i.e., 〈u〉 = 0, σ_u = 0). This maps directly back to the original Jansen–Rit model with pyramidal-only input. Matching the effective noise level used by Jansen and Rit (correcting a scaling error in the original paper) is achieved by allowing p to fluctuate with standard deviation σ_p = 0.5390 s⁻¹.

To examine the difference between cases where input is provided in different ratios to the spiny stellate population and pyramidal population, we continue the bifurcation point H1 in the (u, p) plane in parameter space, using the numerical continuation package MATCONT (Dhooge et al., 2003).

Figure 2 shows the bifurcation diagram in the (u, p) plane. This plane is a two-dimensional slice through the larger parameter space of the model, so that a curve in this plane corresponds to a surface in parameter space. The Hopf curve is almost horizontal for u > 0, implying that the level of pyramidal cell stimulation required to reach the supercritical Hopf bifurcation in the model (p ∼ 75–90 s⁻¹) is roughly independent of the level of spiny stellate cell stimulation for u ≳ 0. For comparison with H1, we select point H2 on this same surface of supercritical Hopf points, but this time with greater mean input to the spiny stellate population (〈u〉 = 270 s⁻¹) than to the pyramidal population (〈p〉 = 73 s⁻¹). The magnitude of fluctuations $\sqrt{{\sigma }_u^2+{\sigma }_p^2}$ in the input is kept the same as at H1, with a standard deviation of σ_u = 0.5203 s⁻¹ in the spiny stellate input and σ_p = 0.1407 s⁻¹ in the pyramidal input.

FIGURE 2

Figure 2. Bifurcation diagram in the (u, p) parameter plane showing Hopf curve (thick, solid curve) and the location of the chosen Hopf bifurcation points H1, H2, and H3 on this curve. A generalized Hopf bifurcation (GH) marks the transition from subcritical Hopf points (the curve below GH) to supercritical Hopf points (the curve continuing beyond GH). Regions where p < 0 or u < 0 are non-physical. Below the Hopf curve (regions I and II) a stable fixed point exists, which gradually loses linear stability as the curve is approached. Above the Hopf curve (region III) this point has lost linear stability and become a stable limit cycle. The dashed line is a curve of fold bifurcation points. In region II a single stable fixed point exists. In region I the system is bistable with a second stable fixed point also existing, at lower excitation.

Autocorrelation Indicator Behaves Differently at H1 and H2

For each of the bifurcation points we simulate the dynamics at four locations in parameter space: the approach to the bifurcation from the linearly stable side (two points), at the bifurcation point (one point), and beyond the bifurcation (one point). In each case the output pyramidal time series (v₂ − v₃) is the focus of our analysis.

Each simulation is performed separately with parameter values fixed at these different values, rather than performing a single dynamic simulation with sliding parameters. This approach allows the time series analyzed at a fixed parameter value to be approximately stationary (provided the total time simulated is long enough) so that statistics for the process at that parameter point can be estimated from a finite time series. Where variance is considered as an indicator of instability this approach also avoids any spurious short-time increases in variance due to the dynamically shifting range of the system in phase space, as distinct from increased noise-driven variance at the new parameter values (Kuehn, 2011).

To determine the effect of proximity to a bifurcation on the fluctuation statistics, we analyze the approach and passage through bifurcations H1 and H2. Figure 3 shows the results for bifurcation H1. An exemplar pyramidal time series (Figure 3A) reveals a fluctuating oscillatory system, whose power spectrum (Figure 3B) peaks at the frequency of the Hopf instability, namely 11 Hz. The series of panels in Figure 3C shows that when approaching and passing point H1 (from left-to-right), the autocorrelation time stays approximately constant.

FIGURE 3

Figure 3. (A) Example section of simulated time series at the bifurcation point H1. (B) Power spectrum at H1. (C) Autocorrelation amplitude at points in parameter space approaching the bifurcation point H1 (<p> = 74.8 s⁻¹, <p> = 84.8 s⁻¹), at the bifurcation point H1 (<p> = 89.8 s⁻¹) and beyond the bifurcation point (<p> = 94.8 s⁻¹). The line indicates the mean over 16 trials and the gray area indicates one standard deviation.

For comparison, the corresponding analyses for bifurcation H2 are shown in Figure 4. By eye, the fluctuation envelope of the amplitude appears smoother. As is evident in Figure 4C, the autocorrelation amplitude decays much more slowly as the system approaches the bifurcation. At the bifurcation point H2 the autocorrelation remains above 20% of its zero-lag value at a lag of 15 s.

FIGURE 4

Figure 4. (A) Example section of simulated time series at the bifurcation point H2. (B) Power spectrum at H2. (C) Autocorrelation amplitude at points in parameter space approaching the bifurcation point H2 (<p> = 58.0 s⁻¹, <p> = 68.0 s⁻¹), at the bifurcation point H2 (<p> = 73.0 s⁻¹) and beyond the bifurcation point (<p> = 78.0 s⁻¹). The line indicates the mean over 16 trials and the gray area indicates one standard deviation.

The variance of the output pyramidal time series increases as the bifurcation H1 is approached, with standard deviations of 0.4550, 0.5344, 0.5630, and 0.6110 mV at the four parameter points respectively. Approaching H2 this also occurs, with standard deviations of 0.1454, 0.2160, 0.2582, and 0.3449 mV respectively for the output time series. It is notable that in the vicinity of point H2, the standard deviation of the simulated pyramidal output time series is between 1.8 and 3.1 times smaller than in the vicinity of point H1, while autocorrelation times are roughly 7 times longer than at H1.

As expected for a Hopf bifurcation in a stochastic system the dynamics change gradually and continuously through the bifurcation (Rowat and Greenwood, 2011). The amplitude of oscillations increases when moving toward and beyond the bifurcation point as revealed by increased variance of the output time series. Close to the bifurcation point this reflects weakening of the stability of the (fixed point) attractor while beyond the bifurcation point it reflects increasing size of the (limit cycle) attractor. The increase in power at 11 Hz is visible in the power spectrum (Figure 3B).

The comparison of Figures 3C and 4C shows that autocorrelation is a useful and clearly visible indicator of linear instability in the vicinity of point H2, but not for point H1. This is despite these being points on the same surface of bifurcations with the same variance of input fluctuations.

We conjecture that the key difference between H1 and H2 is the orientation of the input fluctuations in phase space with respect to the two-dimensional center manifold of the bifurcation, which determines the specific directions in which linear stability is weakening. When close to the equilibrium point, the center manifold surface can be approximated by the center eigenspace of the bifurcation. Since the eigenvectors of the linearized system are far from orthogonal, the relevant reference plane to determine the noise component projected into the center eigenspace is the plane that is perpendicular to the stable eigenspace. For H1 the resulting projection of the noise vector onto this reference plane is cosα = 0.0031. For H2 the projection is cosα = 0.0008; i.e., the noise input has a projection onto that plane that is four times larger in the case of H1 than in the case of H2.

However, the comparison between points H1 and H2 does not by itself give strong support for this hypothesis, because there are several other factors that are significantly different between H1 and H2. In particular H2 has 3.1 times the total input firing rate of H1, so that on this basis the difference in autocorrelation could simply be due to greater level of excitation for point H2. This motivates the comparison constructed below.

Autocorrelation Depends on Orientation of Input Fluctuations

In order to separate the effect of different mean firing rates from the effect of different noise orientation, we construct two new scenarios H3p and H3u, where the only difference between them is the noise orientation; all other parameters are kept identical. We choose point H3 on the same bifurcation line of supercritical Hopf points, but with equal mean input firing rates to pyramidal and spiny stellate populations (mean input firing rate per neuron of 〈u〉 = 〈p〉 = 80.35 s⁻¹). We simulate two scenarios at point H3 to test the conjecture, with both scenarios using the same values for all model parameters, and in particular with both scenarios using the same mean input firing rates, as illustrated in Figure 5.

FIGURE 5

Figure 5. One second of sample external inputs u and p, (A) for scenario H3p and (B) for scenario H3u.

We define the scenario H3p as the case where only the pyramidal input is allowed to fluctuate about its mean, while spiny stellate input is held steady at its mean value, corresponding to the parameters 〈p〉 = 80.35 s⁻¹, σ_p = 0.5390 s⁻¹, 〈u〉 = 80.35 s⁻¹, and σ_u = 0 s⁻¹.

Scenario H3u is defined as the case where only the spiny stellate input is allowed to fluctuate, while pyramidal input is held steady, corresponding to the parameters 〈p〉 = 80.35 s⁻¹, σ_p = 0 s⁻¹, 〈u〉 = 80.35 s⁻¹, and σ_u = 0.5390 s⁻¹.

By using these two constructed scenarios, all parameters in the simulation are kept identical between scenarios H3p and H3u except for the direction of the fluctuations of input in phase space, which is rotated in phase space from the pyramidal direction to the spiny stellate direction. Rotating the vector of fluctuations independently from the vector of mean inputs is non-physiological. The simulated results of the non-physiological scenarios H3p and H3u are used to shed light on the reason for different autocorrelation in the original realistic scenarios H1 and H2.

Comparable analyses of these two scenarios are presented in Figures 6 and 7. The contrast between scenarios H3p and H3u is clear in Figures 6C and 7C. When the fluctuations are in the input to the pyramidal population (scenario H3p), the decay of autocorrelation amplitude changes little as the bifurcation point is approached. By contrast, when the fluctuations are in the input to the spiny stellate population (scenario H3u) the indicator of increased autocorrelation length is very prominent. A large increase in autocorrelation heralds the transition to linear instability in scenario H3u with significant autocorrelation at lags of up to 450 s. This indicator is much less evident in scenario H3p.

FIGURE 6

Figure 6. (A) Example section of simulated time series for scenario H3p. (B) Power spectrum at H3. (C) Autocorrelation amplitude at points in parameter space approaching the bifurcation point H3 (<p> = 65.3 s⁻¹, <p> = 75.3 s⁻¹) at the bifurcation point H3 (<p> = 80.3 s⁻¹) and beyond the bifurcation point (<p> = 85.3 s⁻¹). The line indicates the mean over 16 trials and the gray area indicates one standard deviation.

FIGURE 7

Figure 7. (A) Example section of simulated time series for scenario H3u. Note the y-axis scale is much smaller than that of Figure 6A, reflecting much smaller output variance in this case. (B) Power spectrum at H3. (C) Autocorrelation amplitude at points in parameter space approaching the bifurcation point H3 (<p> = 65.3 s⁻¹, <p> = 75.3 s⁻¹), at the bifurcation point H3 (<p> = 80.3 s⁻¹) and beyond the bifurcation point (<p> = 85.3 s⁻¹). The line indicates the mean over 16 trials and the gray area indicates one standard deviation.

It is also instructive to view the autocorrelation amplitude with log scaling of the delay axes. The results for the four scenarios (H1, H2, H3p, H3u) we have thus far considered are shown in Figure 8. Whereas the autocorrelation length stays almost invariant across the bifurcation in scenario H1 (Figure 8A), a clear increase is seen in scenario H2 (Figure 8B). Where scenario H3p shows a small, but systematic lengthening (Figure 8C), a progression through the same points in parameter space – but now with input fluctuations aligned with the stable eigenspace – can again be seen to lead to a dramatic increase (Figure 8D).

FIGURE 8

Figure 8. Autocorrelation amplitude in log(delay)-linear(correlation) coordinates. Each panel shows step immediately before (blue), at (green), and beyond (orange) Hopf bifurcation. (A) Scenario H1, (B) H2, (C) H3p, (D) H3u.

In both scenarios H3p and H3u the variance of the output pyramidal time series increases as the bifurcation point H3 is approached, with standard deviations of 0.4884, 0.4730, 0.5134, and 0.5339 mV for H3p and 0.0188, 0.0295, 0.0873, and 0.3033 mV for H3u. The standard deviation is starkly different between these two scenarios, with standard deviation between 2 and 22 times smaller in scenario H3u than in scenario H3p. Thus changing the noise input direction results in both reduced variance and increased autocorrelation length. Variance also increases more rapidly in scenario H3u than H3p as the bifurcation is approached.

Relating this to the orientation of input noise, the contrast in alignment is even greater between scenarios H3p and H3u than in the comparison of H1 and H2. For H3p the projection of the input noise onto the reference plane perpendicular to the stable eigenspace is cosα = 0.003, i.e., noise input has a non-negligible component perpendicular to the stable eigenspace of the bifurcation near the equilibrium point, whereas for H3u the projection is cosα = 0.00006, i.e., the projection of the noise input onto the reference plane is 50 times smaller in the case of H3u than H3p.

Input Correlation and Output Variable not Important

The results presented above are calculated from the pyramidal time series. Applying the same process to time series for the other two populations in the model (spiny stellate and inhibitory) shows that in each case, the results for autocorrelation decay and variance show the same behavior as the pyramidal time series. This is important, as it rules out the possibility that the autocorrelation difference results from different amounts of filtering between the noise input and the measured output.

In the case of point H2, both pyramidal input and spiny stellate input have fluctuations. To check whether correlations between the two input fluctuations are important to the results we examine the two extreme cases of independent and perfectly correlated inputs. Cross-correlation of input fluctuations does not affect the results: autocorrelation amplitude of the pyramidal output decays over a similar time scale whether the inputs to the two populations are independent or perfectly correlated.

Scaling Properties of Output Fluctuations

As reviewed earlier, long-tailed fluctuation distributions have been observed in the amplitude fluctuations of alpha (Freyer et al., 2009) and beta oscillations in scalp EEG data (Linkenkaer-Hansen et al., 2004). Therefore we study the statistical properties of fluctuations at and near the bifurcation points in scenarios H3u and H3p. In particular, we characterize fluctuations by the distributions of sizes and durations of excursions in the amplitude envelope of the detrended pyramidal time series. More specifically, we analyze the squared Hilbert amplitude, which is a measure of instantaneous power. We extract excursions above a threshold (sometimes termed “avalanches” in the literature) where each excursion is delineated by the time points at which the instantaneous power crosses the threshold from below and the next crossing from above. Fluctuation duration is thus the length of the time interval for which the power is above threshold, and we define fluctuation size to be the time integral of the instantaneous power over this interval (i.e., the area under the curve, a measure of energy in the fluctuation). We choose the threshold for each time series such that it approximately maximizes the number of identified fluctuations and falls in a regime where the fluctuation statistics are relatively insensitive to small changes in this value.

We analyze the fluctuation size and duration distributions following the methods of Clauset et al. (2009). For each set of fluctuation statistics we calculate the inverse cumulative distribution function and fit candidate distributions to the tail using the method of maximum likelihood: power law (the Pareto distribution), power law with exponential cutoff, lognormal, and exponential. Here the tail is all the data above a lower bound that we identify as the value that minimizes the Kolmogorov–Smirnov goodness-of-fit statistic between the power law model and the data (Clauset et al., 2009). This method of determining the range of the fit from the data strikes a balance between fitting too wide a range (i.e., outside the power law regime) and too narrow a range (i.e., throwing away data unnecessarily). We use the same fitting range for all four candidate distributions. We estimate a p-value for the fitted power law by comparing the data to 1,000 synthetic data sets drawn from a true power law, which accounts for whether the deviation between the data and the fitted power law is within the range expected for finite sampling of a true power law. The p-value is taken to be the fraction of synthetic data sets that deviate from the power law by at least as much as the data, and p > 0.1 indicates plausibility of the power law hypothesis (Clauset et al., 2009). We compare the fitted power law with alternative distributions using likelihood ratio tests. Significant deviation of the likelihood ratio from zero is tested using Vuong’s methods (Vuong, 1989). For the nested hypothesis of power law versus power law with cutoff (the latter family includes the former), the null hypothesis is that the power law is best-fitting distribution. For all other tests, the null hypothesis is that both distributions are equally far from the true distribution.

Figure 9 shows the fluctuation distributions for H3u. The empirical distributions for both duration (Figure 9A) and area (Figure 9B) exhibit a scaling regime over approximately four orders of magnitude. The power law fits for duration and area have exponents 1.56 and 1.51, and p-values p = 0.27 and p = 0.72, respectively, and are thus consistent with the hypothesis that the true distribution is a power law. The fitted exponents depend weakly on the threshold value but the main finding of a broad scaling regime is unchanged. The lognormal and power law with exponential cutoff are also consistent with the data: the likelihood ratio tests do not distinguish between the lognormal and power law fits (duration: p = 0.15; area: p = 0.26), but favor the power law with exponential cutoff over both power law (duration: p = 0.016; area: p = 0.047) and lognormal (duration: p = 0.003; area: p = 0.004). The pure exponential distribution is strongly ruled out in all cases (p ≪ 0.001) and so is not shown.

FIGURE 9

Figure 9. Upper cumulative distributions of fluctuation statistics at the bifurcation point H3u, using squared Hilbert amplitude thresholded at 0.008 mV², with power law (red), power law with exponential cutoff (green), and lognormal (blue) fits plotted for the fitted range of the tail. (A) Fluctuation duration. (B) Fluctuation size as given by area under the curve.

Approach to this bifurcation, shown in Figure 10, reiterates the autocorrelation results of Section “Autocorrelation Depends on Orientation of Input Fluctuations.” Near H3u (Figure 10A), the long scaling regime of Figure 9A (black) is significantly diminished away from the bifurcation (red), with few fluctuations having durations >10 s. Here, the pure power law is ruled out (p < 0.001), and the power law with exponential cutoff is strongly favored over all the alternatives tested. For comparison, Figure 9B shows fluctuations at the same bifurcation when noise enters almost perpendicular to the center eigenspace (scenario H3p). At the bifurcation (black), there is no clear scaling regime, and the distribution is essentially unchanged by moving to a more stable point in parameter space (red). The power law fit is ruled out for both points (p < 0.001), and again the power law with cutoff is strongly favored. Thus, as in the autocorrelation cases (Figures 7C and 8C), the fluctuation statistics clearly herald the approach to the bifurcation for H3u but only negligibly for H3p.

FIGURE 10

Figure 10. Comparison of fluctuation duration distributions between points approaching (p = 65.3 s⁻¹) and at the bifurcation point (p = 80.3 s⁻¹) for the two noise input directions. (A) Point H3u (black, threshold = 0.008 mV²) and nearby more stable point (red, threshold = 0.0002 mV²). (B) Point H3p (black, threshold = 0.4 mV²) and nearby more stable point (red, threshold = 0.18 mV²).

Discussion

The relevance of these results to physiology is twofold. Firstly, we have demonstrated a fundamental limitation in the use of autocorrelation as an indicator of the loss of linear stability, a limitation which will apply when attempting to detect bifurcations from actual human EEG, EMG, and MEG time series. Secondly, the demonstration of both long autocorrelation times and scale-free temporal fluctuations in a simple, low-dimensional stochastic model informs the debate about whether the brain exhibits self-organized criticality, because it shows that these features can also arise from mechanisms other than a multi-scale critical phase transition.

Close to the supercritical Hopf bifurcation in the Jansen–Rit model, we have shown that when lengthened autocorrelation times and scale-free fluctuations manifest in any one cell population as indicators of approach to the bifurcation, then they are indeed detectable in the pyramidal time series that is most closely associated with EEG signals. The standard parameters of the model provide sufficiently large coupling between the three neural populations that lengthened autocorrelation is evident in all three populations when it is present in any of them.

When considering long time autocorrelation and scale-free fluctuations that are present in human EEG time series this suggests that in addition to the possibility that these could arise in the brain at the point of phase transition in a complex, multi-scale system (Linkenkaer-Hansen et al., 2004; Stam and de Bruin, 2004), there may also be a role for low-dimensional stochastic dynamics at the population scale in generating these indicators.

More importantly, we have shown that even in a very simplistic cortical model, these indicators can already be subtle in their dependency on neuronal inputs. Longer autocorrelation times are not guaranteed to be evident in the output just because there is a bifurcation where linear stability is lost. In particular we have shown that a change of the orientation in phase space of small fluctuations in the input can be sufficient to enhance or almost completely remove this indicator.

Jansen and Rit (1995) suggested that input to excitatory interneurons could be removed from the model, as input to the pyramidal population from coupled columns was expected to have the same effect. Our results show that when fluctuations in the input are taken into account, the statistical properties of the model output are sensitive to the choice of which neural population receives the extrinsic input.

The Jansen–Rit model is representative of a broad class of models that mathematically can be expressed as a composition of sigmoid functions and second order linear filters. It is worth noting that neural field models (such as Jirsa and Haken, 1996; Robinson et al., 1997), when restricted to spatially uniform solutions, can also be expressed in this mathematical form, with an additional critically damped second order linear filter capturing the time characteristics of local axonal propagation with a population spread of sources and axon parameters (Robinson et al., 1997). There are no particularities of the current model that suggest that the phenomena which we describe will be limited to this setting. The present results regarding fluctuation orientation hence speak broadly to the commonly employed neural mass and neural field models of large-scale neuronal activity.

Opposite Effect on Autocorrelation and Variance

Autocorrelation and variance of the output signal have been suggested as generic indicators of the approach to local bifurcation, as standard linear analysis shows they are both expected to increase as the bifurcation is approached and the real part of bifurcating eigenvalues approaches zero. We also observed that changing the orientation of input fluctuations can result in autocorrelation increasing at the same time as variance is decreased. Insight into these phenomena can be gained by considering the behavior of a simple low-dimensional linear stochastic system. In the one-dimensional linear case of an Ornstein–Uhlenbeck process, dx = −axdt + bdW, the (normalized) autocorrelation is given by exp(−aτ) and variance by 1/2 b²/a, so both increase as the size of the eigenvalue a approaches zero. In particular, variance increases linearly with the variance of noise input b² (Gardiner, 2010). The same is true for a linearized two-dimensional system near a Hopf bifurcation (Steyn-Ross et al., 2006). If we naively assume that aligning input noise with the center eigenspace increases the amount of noise affecting the slow dynamical system of the center manifold, we would expect variance to be greater when the angle with the center eigenspace plane was smaller, which was not the case. From consideration of the normal form transformation (Roberts, 2008) it is rather the plane perpendicular to the stable eigenspace that should be relevant in determining the magnitude of noise driving the slow dynamics. Because the eigenvectors of the linearized Jansen–Rit system are far from orthogonal, that reference plane is almost normal to the center eigenspace plane, resulting in the observed reversal of the expected relationship between noise orientation and variance.

It may be possible to study these bifurcation indicators more specifically in a normal form model by considering a full center manifold reduction. Close to the bifurcation non-linear terms can result in multiplicative noise in the slow dynamical system of the center manifold (Roberts, 2008). These occur in addition to the simpler additive noise that results directly from linear transformation of the input noise terms but so far we have yet to calculate the magnitude or importance of these multiplicative noise terms in the present system. Furthermore, any local analysis of the behavior close to the equilibrium point is valid only for the case of small noise, so that the state of the system remains local to the equilibrium point. That is not necessarily the case for this system, as suggested by larger output standard deviation near the bifurcation seen in the cases of H1 and H3p, which is comparable to the amplitude of the subsequent limit cycles. This implies that the system is exploring a wider region of phase space compared to the cases with high autocorrelation (H2 and H3u). Thus the structure of flow in the phase space further from the equilibrium point may be directly responsible for the quickly decaying autocorrelation in those cases. In particular if the center manifold curves away from the center eigenspace, then at a sufficient distance from the equilibrium point the directions of noise input which are “well aligned” and “poorly aligned” with the manifold may be reversed.

Power Law Scaling of Output Fluctuations

Analyzing the distributions of fluctuation sizes and durations, we observe the presence of a long power law scaling regime that extends over four orders of magnitude with a rapid truncation at the far right hand tail at the bifurcation when input fluctuations are normal to the reference plane. This power law scaling is not observed when input fluctuations have a significant projection onto the reference plane. Further away from the bifurcation, the power law regime extends for less than one order of magnitude so that the lengthy power law tail provides a signature of proximity to the bifurcation in that scenario.

A range of simple dynamical mechanisms are known to permit production of scale-free fluctuation structure of this kind. A relaxation process with a fractional operator formally yields a power law (Pareto) probability distribution of fluctuation durations (Sokolov and Klafter, 2005). Multiplicative noise (which arises when reducing oscillation dynamics of the model to two dimensions) can also in specific cases result in power law probability distributions (Anteneodo and Riera, 2005). However, the cause of the power law scaling of the distributions of fluctuations in our system is not yet determined.

Future Work

This study considered autocorrelation in the output of a single Jansen–Rit model region, representing a small area of cortex of the order of 2–3 mm². For the question of potential detectability in EEG it remains to examine the effect on autocorrelation of combining the output of a large number of cortical regions, whose oscillations may be synchronized to a greater or lesser degree and where the output measurement function relating EEG to the combination of sources plays an important role.

Within the Jansen–Rit model we also observed indicators close to other bifurcation types, including switching between attractors in a bistable region near a cusp bifurcation and “flickering” or intermittent switching away from a stable fixed point in a monostable region near a sniper bifurcation, which are not explored further in this paper. Therefore it remains to examine the sensitivity of these and other indicators, such as mean switching times as bifurcations are approached, to noise orientation.

It is hoped that normal form analyses near the bifurcation will shed some light on the mechanism by which the input noise affects autocorrelation. A first step will be to examine a simpler normal form system displaying the same behavior, where exact control over the shape of the center manifold can be afforded, initially targeting the limiting case of small fluctuations. Such an analysis will serve to separate the generic local effects of the Hopf bifurcation from global behavior due to excursion of the state further from the equilibrium point.

Examination of a normal form system will also be key to determining the reason for the power law scaling of fluctuation statistics. The results presented in this paper show that some of the indicators of instability reported in human EEG also arise in the output of a simple neural mass model near linear instability.

While similar indicators can also emerge from a critical phase transition in a complex, multi-scale system, we have shown in the present study that some of the same indicators can arise in a very different way, from the low-dimensional stochastic dynamics at a single scale: the mesoscopic scale of interacting populations. As the field advances, it will become increasingly important to move away from a single umbrella notion of “criticality” in brain dynamics toward defining a number of exact, and possibly distinct, mechanisms responsible for correlations and scale-free fluctuations in the time and/or spatial domains. It is certainly possible at this stage that multiple mechanisms play a role.

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

Michael Breakspear acknowledges the support of the Australian Research Council, the National Health and Medical Research Council and the Brain Network Recovery Group Grant JSMF22002082. Tjeerd W. Boonstra acknowledges support from the Netherlands Organization for Scientific Research (NWO #451–10-030). The authors thank Saeid Mehrkanoon for helpful technical advice.

References

Anteneodo, C., and Riera, R. (2005). Additive-multiplicative stochastic models of financial mean-reverting processes. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 72, 026106.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bedard, C., Kroeger, H., and Destexhe, A. (2006). Does the 1/f frequency scaling of brain signals reflect self-organized critical states? Phys. Rev. Lett. 97, 118102.

Braun, J., and Mattia, M. (2010). Attractors and noise: twin drivers of decisions and multistability. Neuroimage 52, 740–751.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Breakspear, M., Roberts, J. A., Terry, J. R., Rodrigues, S., Mahant, N., and Robinson, P. A. (2006). A unifying explanation of primary generalized seizures through nonlinear brain modeling and bifurcation analysis. Cereb. Cortex 16, 1296.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Breakspear, M., and Terry, J. R. (2002). Detection and description of non-linear interdependence in normal multichannel human EEG data. Clin. Neurophysiol. 113, 735–753.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Clauset, A., Shalizi, C. R., and Newman, M. E. J. (2009). Power-law distributions in empirical data. SIAM Rev. Soc. Ind. Appl. Math. 51, 661–703.

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

David, O., Harrison, L., and Friston, K. J. (2005). Modelling event-related responses in the brain. Neuroimage 25, 756–770.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

CrossRef Full Text

Dhooge, A., Govaerts, W., and Kuznetsov, Y. A. (2003). MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs. ACM Trans. Math. Softw. 29, 141–164.

CrossRef Full Text

Freeman, W. (1992). “Predictions on neocortical dynamics derived from studies in paleocortex,” in Induced Rhythms in the Brain, Brain Dynamics Series, Chap. 9, eds E. Basar and T. H. Bullock (Boston: Birkhäuser), 183–199.

Freyer, F., Aquino, K., Robinson, P. A., Ritter, P., and Breakspear, M. (2009). Bistability and non-Gaussian fluctuations in spontaneous cortical activity. J. Neurosci. 29, 8512.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Freyer, F., Roberts, J. A., Becker, R., Robinson, P. A., Ritter, P., and Breakspear, M. (2011). Biophysical mechanisms of multistability in resting-state cortical rhythms. J. Neurosci. 31, 6353.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Freyer, F., Roberts, J. A., Ritter, P., and Breakspear, M. (2012). A canonical model of multistability and scale-invariance in biological systems. PLoS Comput. Biol. doi:10.1371/journal.pcbi.1002634

CrossRef Full Text

Friston, K. J., and Dolan, R. J. (2010). Computational and dynamic models in neuroimaging. Neuroimage 52, 752–765.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Friston, K. J., Harrison, L., and Penny, W. (2003). Dynamic causal modeling. Neuroimage 19, 1273–1302.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gardiner, C. (2010). Stochastic Methods: A Handbook for the Natural and Social Sciences, Springer Series in Synergetics Series. New York: Springer.

Grimbert, F., and Faugeras, O. (2006). Bifurcation analysis of Jansen’s neural mass model. Neural Comput. 18, 3052–3068.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Haken, H., Kelso, J. A. S., and Bunz, H. (1985). A theoretical model of phase transitions in human hand movements. Biol. Cybern. 51, 347–356.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

He, B. J., Zempel, J. M., Snyder, A. Z., and Raichle, M. E. (2010). The temporal structures and functional significance of scale-free brain activity. Neuron 66, 353–369.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jirsa, V. K., and Haken, H. (1996). Field theory of electromagnetic brain activity. Phys. Rev. Lett. 77, 960–963.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kelso, J. A. S. (2010). Instabilities and phase transitions in human brain and behavior. Front. Hum. Neurosci. 4:23. doi:10.3389/fnhum.2010.00023

CrossRef Full Text

Kelso, J. A. S., Scholz, J. P., and Schöner, G. (1986). Nonequilibrium phase transitions in coordinated biological motion: critical fluctuations. Phys. Lett. A 118, 279–284.

CrossRef Full Text

Kuehn, C. (2011). A mathematical framework for critical transitions: bifurcations, fast-slow systems and stochastic dynamics. Physica D 240, 1020–1035.

CrossRef Full Text

Linkenkaer-Hansen, K., Nikouline, V. V., Palva, J. M., and Ilmoniemi, R. J. (2001). Long-range temporal correlations and scaling behavior in human brain oscillations. J. Neurosci. 21, 1370.

Pubmed Abstract | Pubmed Full Text

Linkenkaer-Hansen, K., Nikulin, V. V., Palva, J. M., Kaila, K., and Ilmoniemi, R. J. (2004). Stimulus-induced change in long-range temporal correlations and scaling behaviour of sensorimotor oscillations. Eur. J. Neurosci. 19, 203–218.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lopes da Silva, F. H., Hoeks, A., Smits, H., and Zetterberg, L. H. (1974). Model of brain rhythmic activity. Biol. Cybern. 15, 27–37.

Marreiros, A. C., Daunizeau, J., Kiebel, S. J., and Friston, K. J. (2008). Population dynamics: variance and the sigmoid activation function. Neuroimage 42, 147–157.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Miller, K. J., Sorensen, L. B., Ojemann, J. G., and den Nijs, M. (2009). Power-law scaling in the brain surface electric potential. PLoS Comput. Biol. 5, e1000609. doi:10.1371/journal.pcbi.1000609

CrossRef Full Text

Moran, R. J., Kiebel, S. J., Stephan, K. E., Reilly, R. B., Daunizeau, J., and Friston, K. J. (2007). A neural mass model of spectral responses in electrophysiology. Neuroimage 37, 706–720.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Orfanidis, S. J. (1996). Optimum Signal Processing, 2nd Edn. New York: McGraw-Hill.

Roberts, A. J. (2008). Normal form transforms separate slow and fast modes in stochastic dynamical systems. Physica A 387, 12–38.

CrossRef Full Text

Roberts, J. A., and Robinson, P. A. (2012). Corticothalamic dynamics: structure of parameter space, spectra, instabilities, and reduced model. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 85, 011910.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Robinson, P. A., Rennie, C. J., Phillips, A. J. K., Kim, J. W., and Roberts, J. A. (2010). Phase Transitions in Physiologically-Based Multiscale Mean-Field Brain Models, Modeling Phase Transitions in the Brain. New York: Springer, 179–201.

Robinson, P. A., Rennie, C. J., and Rowe, D. L. (2002). Dynamics of large-scale brain activity in normal arousal states and epileptic seizures. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 65, 041924.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Robinson, P. A., Rennie, C. J., and Wright, J. J. (1997). Propagation and stability of waves of electrical activity in the cerebral cortex. Phys. Rev. E Stat. Phys. Plasmas. Fluids Relat. Interdiscip. Topics 56, 826.

CrossRef Full Text

Rowat, P. F., and Greenwood, P. E. (2011). Identification and continuity of the distributions of burst-length and interspike intervals in the stochastic Morris-Lecar neuron. Neural Comput. 23, 1–31.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rümelin, W. (1982). Numerical treatment of stochastic differential equations. SIAM J. Numer. Anal. 19, 604–613.

CrossRef Full Text

Scheffer, M., Bascompte, J., Brock, W. A., Brovkin, V., Carpenter, S. R., Dakos, V., Held, H., van Nes, E. H., Rietkerk, M., and Sugihara, G. (2009). Early-warning signals for critical transitions. Nature 461, 53–59.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Scholz, J. P., Kelso, J. A. S., and Schöner, G. (1987). Nonequilibrium phase transitions in coordinated biological motion: critical slowing down and switching time. Phys. Lett. A 123, 390–394.

CrossRef Full Text

Sokolov, I. M., and Klafter, J. (2005). From diffusion to anomalous diffusion: a century after Einstein’s Brownian motion. Chaos 15, 026103–026103.

CrossRef Full Text

Sotero, R. C., Trujillo-Barreto, N. J., Iturria-Medina, Y., Carbonell, F., and Jimenez, J. C. (2007). Realistically coupled neural mass models can generate EEG rhythms. Neural Comput. 19, 478–512.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stam, C. J., and de Bruin, E. A. (2004). Scale-free dynamics of global functional connectivity in the human brain. Hum. Brain Mapp. 22, 97–109.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stam, C. J., Pijn, J. P. M., Suffczynski, P., and da Silva, F. H. L. (1999). Dynamics of the human alpha rhythm: evidence for non-linearity? Clin. Neurophysiol. 110, 1801–1813.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Steyn-Ross, D. A., Steyn-Ross, M. L., Wilson, M. T., and Sleigh, J. W. (2006). White-noise susceptibility and critical slowing in neurons near spiking threshold. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 74, 051920.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

van Albada, S. J., Kerr, C. C., Chiang, A. K. I., Rennie, C. J., and Robinson, P. A. (2010). Neurophysiological changes with age probed by inverse modeling of EEG spectra. Clin. Neurophysiol. 121, 21–38.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Vuong, Q. H. (1989). Likelihood ratio tests for model selection and non-nested hypotheses. Econometrica 57, 307–333.

CrossRef Full Text

Wendling, F., Bartolomei, F., Bellanger, J. J., and Chauvel, P. (2002). Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition. Eur. J. Neurosci. 15, 1499–1508.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wilson, H. R., and Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1–24.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Zavaglia, M., Astolfi, L., Babiloni, F., and Ursino, M. (2006). A neural mass model for the simulation of cortical activity estimated from high resolution EEG during cognitive or motor tasks. J. Neurosci. Methods 157, 317–329.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Appendix

FIGURE A1

Figure A1. Power spectrum at H1, using a larger window size of 150 s (750,000 samples) to show lower frequencies from 6.7 × 10⁻³ Hz.

FIGURE A2

Figure A2. Power spectrum at H2, using a larger window size of 150 s (750,000 samples) to show lower frequencies from 6.7 × 10⁻³ Hz.

FIGURE A3

Figure A3. Power spectrum for scenario H3p, using a larger window size of 450 s (2,250,000 samples) to show lower frequencies from 2.2 × 10⁻³ Hz.

FIGURE A4

Figure A4. Power spectrum for scenario H3u, using a larger window size of 450 s (2,250,000 samples) to show lower frequencies from 2.2 × 10⁻³ Hz.