Higher-Order Interactions Characterized in Cortical Activity

Coherence potential analysis suggests an intrinsic threshold in nLFP pattern formation. A, Coherence potential analysis utilizes the full baseline-to-baseline excursion of the nLFP waveform (red). B, The nLFP waveform at a source electrode is compared with the simultaneously recorded LFP at other electrodes (black rectangle). Similarities for each pairwise comparison are quantified by the Pearson correlation coefficient, r. Comparisons for randomly chosen, length-matched LFP segments from the source electrode serve as control (green rectangle). C, nLFP amplitude correlates with transient increase of local unit firing. Normalized increase of firing (see Materials and Methods) plotted as a function of nLFP amplitude (mean ± SD; n = 56 units). For all measured nLFP amplitudes, the increase of firing is significant (p ≪ 10⁻⁶; t test). D, Similarity in nLFP waveforms at distant sites increases nonlinearly with nLFP peak amplitude. Plotted are the fractions of electrodes on the array with high similarity r > 0.8, i.e., R_min = 0.8, with the source electrode nLFP as a function of nLFP amplitude. Average increase in nLFP waveform similarity over all source channels (black; n = 91). Expectation in similarity for random comparisons (green). Threshold (−2.5 SD; arrow) used for nLFP detection in the current neuronal avalanche analysis. E, The nonlinear coherence potential function is revealed for high-similarity requirements, i.e., R_min > 0.6. Coherence potential probability plotted as a function of minimal similarity r > R_min (color coded). F, Difference between functions in E and expectation from random controls.

DG model. A, The distribution of a 3D Gaussian and its projections at individual 2D subspaces. Marginal distributions at a 2D space are represented by probability density contours. Λ, pairwise covariance. B, Converting continuous Gaussian variables to binary variables (events) by thresholding. Three Gaussian variables with pairwise covariance as specified in A and without higher-order interactions (left). Individual thresholds (thr_i, thr_j, and thr_k; broken lines) applied to each variable. The resulting trains of events, i.e., nLFPs or spikes, after thresholding (right). To fit the DG model to the data, the pairwise covariance (Λ) of the multidimensional continuous Gaussian and each threshold need to be adjusted such that the resulting binary variables have event rates and pairwise correlations identical to the data. We note that there is no direct correspondence between the threshold to extract nLFP patterns and the individual thresholds to adjust the rates of the binary variables.

The DG model accurately reconstructs neuronal avalanche dynamics and predicts second- and third-order interactions in neuronal avalanches significantly better than the Ising model. A, The DG model (blue) reconstructs the power law in size distribution for avalanche patterns (dots) more accurately than the Ising model (red). Results for data, Ising model, and Independent (Ind) model replotted from Figure 1G for comparison. B, The DG model is also superior in predicting the probabilities of avalanche patterns (corresponding data from A). Observed pattern probability P is plotted against model predictions. Solid line indicates equality. The most common pattern (all zeros; inactive) is not shown for visual clarity. C, Quantification of model prediction demonstrates 1–2 orders of magnitude in improvement of the DG model over the Ising model. Cumulative distribution of JS divergence between the observed and predicted size distribution of avalanche patterns for n = 30 randomly chosen, spatially compact 10-electrode groups (left; monkey A; compare Fig. 1H). Corresponding analysis for avalanche pattern probabilities (right). One-half of randomly chosen bins was used to predict the other half of data (half-data; see Materials and Methods) D, The DG model accurately predicts second- (θ_ij) and third-order (θ_ijk) interactions. In comparison, the Ising model is less accurate for second-order interactions (right arrow) and fails by design to predict third-order interactions (ellipsoid, arrow). Measured interactions of first to third order (θ_i, θ_ij, and θ_ijk, respectively) are plotted against model predictions for n = 98 randomly chosen 3-electrode groups (monkey A). E, Corresponding measured avalanche pattern probability plotted against model predictions (over all 3-electrode groups). Solid line indicates equality. F, The DG model accurately predicts the power law in the size distribution for avalanche patterns for systems much larger than 10 sites. Prediction for n = 24, 47, and 91 electrode sites are shown (compare Fig. 1F).

Demonstration of significant temporal correlations in avalanche dynamics. A–B, Size distribution for spatiotemporal clusters from bin-shuffled nLFPs (gray broken line) compared with the original avalanche size distribution (black solid line). For bin-shuffled data, temporal correlations beyond Δt, i.e., the duration of a single bin, are removed, which destroys the power law in cluster sizes. A, 10-electrode group. B, 91-electrode group. Insets, areas used for analysis.

The DG model accurately predicts instantaneous nLFP patterns. The dataset used in A–E are the same as in Figure 4A–F, and all analysis was performed in a corresponding way, except that activity patterns were based on instantaneous nLFP patterns found within a single bin of duration Δt, without further temporal concatenation into spatiotemporal clusters. For details, see Figure 4 legend.

Summary of the comparison between the DG and Ising models in approximating nLFP patterns for both monkeys. The DG model provides much more accurate predictions for both avalanche patterns (A) and instantaneous nLFP patterns (B) (monkey A, gray; monkey B, black). JS divergence of Ising model in predicting pattern probabilities is plotted against that of the DG model (10-electrode groups; spatially compact: filled circle, 30 groups; random chosen: unfilled circle, 30 groups). Solid line indicates equality. Insets: Results for approximating the size distribution of avalanche patterns (A) and instantaneous nLFP patterns (B).

Relation between LFP, instantaneous LFP patterns, and variables of the DG model. A, LFP amplitude distribution fitted by a Gaussian function. Whereas the fit seems good at linear scale (left), a clear deviation is visible for amplitudes larger than ±2–3 SD in semilogarithmic coordinates (right). −2.5 SD threshold (broken line). Example from single electrode. B, The pairwise correlation in the continuous LFP correlates with the pairwise correlation found for the corresponding hidden variables in the DG model (n = 91 electrodes; −2.5 SD nLFP threshold). Note the systematic, absolute difference in correlation between the hidden variables and corresponding LFP sites. Solid line, identity. C, Performance of the DG model when the covariances estimated at nLFP threshold T are used to predict nLFP patterns at threshold T′. Here T was varied from −1.5 to −3.5 SD to predict probability of nLFP patterns at T′ = −2.5 SD (an example 10-electrode group; solid line indicates equality). The predictions made by the Ising model (red) are provided for comparison. D, Quantification of results in C. JS divergence between the actual pattern probability distribution and that predicted by the DG model is plotted as a function of the threshold T to predict patterns at nLFP threshold T′ = −2.5 SD (mean ± SD; 30 10-electrode groups analyzed). Both matched (T = T′) and mismatched (T ≠ T′) DG models significantly outperform the (matched) Ising model. However, the mismatched DG models are significantly less accurate than the matched DG model. E, Covariance of the DG model constructed at threshold T = −2.5 SD is used to predict nLFP patterns obtained at T′ = [−1.5, −2, −2.5, −3 SD]. Results for the corresponding matched Ising and DG models are provided for comparison.

The magnitude of higher-order interactions introduced by thresholding depends on event rate and strength of pairwise correlations. A, Third-order interactions in the binary variables cannot be neglected if, for the Gaussian variable, the covariance is strong and the mean is far below or above 0. Change in the magnitude of third-order interaction (θ_ijk) is shown as a function of the mean (γ) and the covariance (Λ) of an underlying 3D Gaussian, u ∼ N (γ, Λ). θ_ijk was calculated for binary variables obtained by applying the threshold u_i > 0. The white dot marks the average mean and covariance of the hidden Gaussians estimated in the DG model for nLFPs (Fig. 4D). B, Differences in pattern probabilities between the Ising and the DG models, here quantified as entropy difference, are most pronounced when, for the binary variables, the rate deviates far away from 0.5 (per bin) and pairwise correlations are strong, as is the case for avalanche patterns (high-rate regime produces a symmetric plot; data not shown). Entropy difference is plotted as a function of average rate and pairwise correlation. Dot (red) marks the average rate and pairwise correlation for avalanche patterns of monkey A (average of 91 channels and all pairs). For both A and B, simulations assume homogenous rate and pairwise correlations for simplicity.

Summary of the performance of the DG and the Ising models in approximating the spiking activity for all five datasets. A, The DG model provides more accurate predictions than the Ising model for both pattern probabilities and pattern sizes, i.e., neuronal synchrony defined as the number of concurrent active units (inset). JS divergence of the Ising model in predicting pattern probabilities and number of concurrent spikes plotted against those of the DG model (15–30 randomly chosen, 10-unit groups for individual datasets; M_A, monkey A; M_B, monkey B; C_A_P1, cat A, probe one; C_A_P2, cat A, probe two; C_B, cat B). B, Similar to A but for neuronal groups with strong pairwise correlations (see Materials and Methods for details about selection). Fifteen to 30 10-unit groups were analyzed for individual datasets. For both A and B, the performance differences between the two models were statistically significant (p ≪ 10⁻⁶ for all comparisons; signed rank test; data combined across all datasets).

For strongly coupled subgroups, the DG model predicts spiking activities in ongoing and stimulus-evoked conditions significantly better than the Ising model. A–C, Ongoing spiking activity during avalanche dynamics in monkey A. A, Predictions in pattern probability by DG and Ising models (30 10-unit groups with strongest pairwise correlation). Zoomed-in view, showing that the DG model more accurately predicted high probability patterns (inset). B, Cumulative distributions of JS divergence for the groups shown in A. C, Measured and predicted interactions for >200 three-unit groups with strongest pairwise correlation. Measured (x-axis) and predicted (y-axis) pattern probabilities (pp) for the same three-unit groups (inset). For more details, see legend of Figure 4B–D. D, Distribution of concurrent spikes for all 56 neurons recorded in monkey A (circles; ongoing activity) and prediction made by the DG model (line). The distribution for the shuffled data was plotted for comparison (dotted line). The shuffling was made for individual neurons and the total number of spikes as well as the distribution of interspike intervals was preserved. E–H, Corresponding results for spiking activity evoked by drifting gratings recorded in area 17 of an anesthetized cat (cat A, probe 1).