Critical fluctuations in cortical models near instability

Abstract

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.

Keywords: Hopf bifurcation; autocorrelation; critical fluctuations; neural mass model.

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.

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.

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.

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.

Figure 5

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

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

(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.

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.

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.

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²).

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

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

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

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.