State equation from the spectral structure of human brain activity††preprint: APS/123-QED††thanks: These authors contributed equally.††thanks: These authors contributed equally.††thanks: These authors contributed equally. ††thanks: These authors contributed equally.
Neural electromagnetic (EM) signals recorded non-invasively from individual human subjects vary in complexity and magnitude. Nonetheless, variation in neural activity has been difficult to quantify and interpret, due to complex, broad-band features in the frequency domain. Studying signals recorded with magnetoencephalography (MEG) from healthy young adult subjects while in resting and active states, a systematic framework inspired by thermodynamics is applied to neural EM signals. Despite considerable inter-subject variation in terms of spectral entropy and energy across time epochs, data support the existence of a robust and linear relationship defining an effective state equation, with higher energy and lower entropy in the resting state compared to active, consistently across subjects. Mechanisms underlying the emergence of relationships between empirically measured effective state functions are further investigated using a model network of coupled oscillators, suggesting an interplay between noise and coupling strength can account for coherent variation of empirically observed quantities. Taken together, the results show macroscopic neural observables follow a robust, non-trivial conservation rule for energy modulation and information generation.
Brain activity varies between states, e.g. in unconsciousness vs consciousness, but a formal macroscopic description remains poorly understood. Spectral analysis of neural electromagnetic (EM) signals indicates that fluctuations occur simultaneously at all measurable frequencies, along a 1/f distribution, above which resonant modes occur at different frequencies Nunez and Srinivasan (2006); Niedermeyer and da Silva (2005); Buzsaki (2006). While neural power spectra display broadband features, they are conventionally investigated by ad hoc partition of the frequency domain. Results of band-limited analyses have demonstrated that the power of bands shifts between brain states. The shift in power within a frequency band is thought at least in part to reflect changes in the level of synchronous neuronal activity at the time scale of that band. Indeed, vigilant behavioral states are commonly called ’desynchronized’ whereas unconscious states ’synchronized’ Niedermeyer and da Silva (2005); Nunez and Srinivasan (2006); Buzsaki (2006).
[rgb]0,0.,0. Macroscopic, EM brain signals measured outside the human head reflect the coordinated activity of tens of thousands of neurons Baillet (2017). Analogies to statistical physics have helped to understand how microscopic properties of neurons and their interactions lead to neural population behavior at the mesoscale Schneidman et al. (2006); Cocco et al. (2009); Marre et al. (2009); Capone et al. (2015); Ferrari (2016); Okun et al. (2015); Gardella et al. (2016); Capone et al. (2018) across brain states Tavoni et al. (2017); Nghiem et al. (2018). Since the use of methodology from statistical physics has uncovered descriptions of individual neurons’ activity and resulting population behaviour at mesoscopic scales, macroscopically observed changes in signal complexity and magnitude could be addressed in analogy to thermodynamics. In fact, for large-scale, whole-brain recordings, spectral entropy measures reportedly vary with levels of consciousness Sitt et al. (2014); Sarasso et al. (2015); McKay et al. (2006); Vanluchene et al. (2004). In the present paper, brain states are described as a function of global state variables, derived from the spectral structure of different human subjects’ brain activity.
Neural signals from each human subject in each brain state may be dissimilar, but the organization of spectra could follow general principles characterized by a state equation. To test this possibility, a statistical analysis of broadband signals is performed using magnetic fields obtained by magnetoencephalography (MEG) from healthy, adult, human subjects in two brain states - resting state with eyes open and fixated (rest) and while performing an N-back visual working memory task (active) Larson-Prior et al. (2013); Van Essen et al. (2012).
Specifically, the magnetic field X, measured for each subject , brain state , sensor , and time epoch where is the total number of epochs (Fig.1) is analyzed. In figure 1a, the time courses of fields for such in rest (, Fig.1a, top) and active (, Fig.1a, bottom) states is shown. The corresponding power spectra are computed for each epoch (Fig.1b).
Observables between resting and active states may present different statistical properties. In order to test this hypothesis, a spectral energy measure is introduced for each subject in each brain state \color[rgb]0,0.,0.and each time epoch, summing over sensors :
[rgb]0,0.,0.Similarly, a spectral entropy measure is defined, for each subject in each brain state and epoch, summing over sensors, as
where, \color[rgb]0,0.,0. describes the relative heights of power spectral frequency bins, from each subject in each brain state .
A negative relation between energy and entropy across time epochs is found (Fig.1c): when energy decreases, entropy increases, suggesting a linear relationship of the following form:
where and are respectively the slope and intercept of the entropy-energy relation across time epochs, for each subject in each brain state. and differ between brain states; in the resting state, energy is higher, entropy is lower, and their relationship follows a steeper (more negative slope) line with a higher-intercept than in the active state (Fig. 1c).
Such a clear linear relationship between the state variables and raises naturally the question of whether it could be a direct consequence of how the quantities are defined. This possibility is tested through manipulation of the signals where the entropy can be modified () by adding zero-mean noise to the power spectrum, with no consequent effect on energy as reported in the inset of Fig.1c. Conversely, multiplying the full spectrum by a constant results in a change of the energy () without altering the entropy. The results support the hypothesis that the observed linear relationship between energy and entropy could reflect fundamental rules that define the way neural activity is organized.
Despite important inter-subject variability in the means over epochs of and , as well as , and , it is verified that the active-rest differences are strongly consistent across subjects, evidence for robustness of the introduced framework (Fig. 2).
Toward a better understanding of mechanisms governing fluctuating oscillations in neural assemblies between brain states, a classical model of coupled oscillators (Kuramoto model) is next studied Kuramoto (1975); Breakspear et al. (2010). N oscillators of phase are arranged over a fully connected network, each one evolving according to the following equation:
where is the coupling strength. The bare frequencies are extracted from a bimodal distribution (sum of two Gaussians). A zero-mean white i.i.d. noise term of standard deviation is received by each oscillator.
From simulated network activity, a collective time-varying variable , to compare with , is defined:
The parameter captures the synchronization between the oscillators - indeed, when the phases are equally distributed over , . Conversly, when the system is completely synchronous, all oscillators have the same phase, so that . The variable represents the time evolution of the system and describes the frequency of synchronous events.
To investigate further the behavior of the effective energy-entropy relation found in human data, the dynamics of the Kuramoto model is investigated by changing the coupling parameter and the amplitude of noise , which are known to strongly modulate synchrony and complexity in coupled oscillator models Pikovsky et al. (2003). Varying , a transition from asynchronous active-like to synchronous rest-like dynamical regime (Fig.3a-c) occurs. To mimic the observed variation between time epochs, the amplitude of noise can be understood as varying in time (Fig.3d-f). Increasing noise amplitude results in less synchronous activity, consistent with previous reports, as well as decreased Pikovsky et al. (2003), reduced energy, and increased entropy (Fig. 3). Altogether, the results indicate that a change in brain state can be modeled by varying coupling strength, while variation across time can be simulated by fluctuations in noise amplitude.
If changes in coupling contribute to empirically observed shifts in neural state variables as predicted by the model, one may expect to find alterations in the phase synchrony in the data. In order to test this prediction, the Phase-Lag Index (PLI), a measure of synchronization across time series is applied to human data Silva Pereira et al. (2017). By definition, the PLI is larger in more phase-locked states. Indeed, in resting compared to active state data, a significant trend toward enhanced phase synchrony is found (Fig.4).
A potential interpretation of the robust relation between these effective state variables could reflect the fact that energy is used to generate information in brain activity. It is important to note that the energy quantity proposed here is an indirect measure of neural activity and may be non-trivially related to actual internal energy. We also find evidence for a change in coupling, due to the consistent change in intercept and slope of the energy-entropy relation. Coupling regulates the sensitivity of networks to noise, which can account for the slope of the aforementioned relation. Indeed, with strong coupling, perturbations propagate more efficiently throughout the network, causing changes in the collective activity and affecting the resulting energy. Conversely, in a weakly-coupled network, perturbations do not easily spread across the system. In other words, the ’gas-like’ regime found in active brains may be useful for the potential to generate information versus the ’liquid-like’ state in which strong coupling promotes a low entropy state with less informational capacity. This provides a new interpretation of previously reported changes in functional connectivity across brain states Biswal et al. (1995); Buckner et al. (2008); Raichle et al. (2001); Massimini et al. (2005); Bettinardi et al. (2015).
Lastly, in light of statistical physics inspired models of neuronal activity Schneidman et al. (2006); Cocco et al. (2009); Marre et al. (2009); Capone et al. (2015); Ferrari (2016); Okun et al. (2015); Gardella et al. (2016); Tavoni et al. (2017); Nghiem et al. (2018); Capone et al. (2018), our results suggest a strategy to bridge the gap from macroscopic system signals studied here to microscopic cellular level activity. Indeed, biophysical models of spiking neural networks could be used to investigate how the mechanisms uncovered with the Kuramoto model translate at the microscopic scale. For instance, biophysical processes that control the variability of noisy inputs or the network’s tendency to synchronize could be studied. This may also provide insight on how the spectral state variables defined here relate to entropy and energy of micro-states in the sense of statistical physics and spin glass models. While such an approach demands scale-integrated data from a sufficiently large population of single neurons simultaneous with global electromagnetic measurements, it would pave the way to more formal, scale-integrated descriptions of brain activity in healthy and pathological states.
Data were provided by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. This project was supported by the European Union (Human Brain Project H2020-720270 and H2020-785907). JSG was supported by the Ludmer Centre for Neuroinformatics and Mental Health, the Molson Neuroengineering Award, The NSERC-CREATE Training Program in Neuroengineering, and the Ann and Richard Sievers Foundation for Innovation in Neuroscience. JML was supported by an NSERC Discovery grant. The authors thank Gleb Bezgin, Bruno Cessac, Alexis Garcia, Peter Grutter, George Kostopoulos, Naj Mahani, Danko Nikolic, Bartosz Teleńczuk, and Núria Tort-Colet for useful discussion.
Ii Author contributions
T.-A. N. wrote data analysis codes, investigated the Kuramoto model, contributed to interpreting the results, and writing the manuscript. J.-M. L. contributed to writing the analysis codes and the manuscript. M. d. V. contributed to interpreting the results, investigating the Kuramoto model, and writing the manuscript. C. C. investigated the Kuramoto model, contributed to interpreting the results and editing the manuscript. A. C. E. provided resources and supervision of the work. A. D. designed research, provided supervision, and edited the manuscript. J. S. G. conceived of the question, designed research, analyzed data, interpreted results, and wrote the paper.
Iii Competing interests
The authors declare no competing interests.
- Nunez and Srinivasan (2006) P. L. Nunez and R. Srinivasan, Electric fields of the brain: the neurophysics of EEG (Oxford University Press, USA, 2006).
- Niedermeyer and da Silva (2005) E. Niedermeyer and F. L. da Silva, Electroencephalography: basic principles, clinical applications, and related fields (Lippincott Williams & Wilkins, 2005).
- Buzsaki (2006) G. Buzsaki, Rhythms of the Brain (Oxford University Press, 2006).
- Baillet (2017) S. Baillet, Nature Neuroscience 20, 327 (2017).
- Schneidman et al. (2006) E. Schneidman, M. J. Berry, R. Segev, and W. Bialek, Nature 440, 1007 (2006).
- Cocco et al. (2009) S. Cocco, S. Leibler, and R. Monasson, Proceedings of the National Academy of Sciences 106, 14058 (2009).
- Marre et al. (2009) O. Marre, S. El Boustani, Y. Frégnac, and A. Destexhe, Physical review letters 102, 138101 (2009).
- Capone et al. (2015) C. Capone, C. Filosa, G. Gigante, F. Ricci-Tersenghi, and P. Del Giudice, PloS one 10, e0118412 (2015).
- Ferrari (2016) U. Ferrari, Physical Review E 94, 023301 (2016).
- Okun et al. (2015) M. Okun, N. A. Steinmetz, L. Cossell, M. F. Iacaruso, H. Ko, P. Barthó, T. Moore, S. B. Hofer, T. D. Mrsic-Flogel, M. Carandini, et al., Nature 521, 511 (2015).
- Gardella et al. (2016) C. Gardella, O. Marre, and T. Mora, eneuro 3, ENEURO (2016).
- Capone et al. (2018) C. Capone, G. Gigante, and P. Del Giudice, arXiv preprint arXiv:1804.00153 (2018).
- Tavoni et al. (2017) G. Tavoni, U. Ferrari, F. P. Battaglia, S. Cocco, and R. Monasson, Network Neuroscience 1, 275 (2017).
- Nghiem et al. (2018) T.-A. Nghiem, B. Telenczuk, O. Marre, A. Destexhe, and U. Ferrari, arXiv preprint arXiv:1801.01853 (2018).
- Sitt et al. (2014) J. D. Sitt, J.-R. King, I. El Karoui, B. Rohaut, F. Faugeras, A. Gramfort, L. Cohen, M. Sigman, S. Dehaene, and L. Naccache, Brain 137, 2258 (2014).
- Sarasso et al. (2015) S. Sarasso, M. Boly, M. Napolitani, O. Gosseries, V. Charland-Verville, S. Casarotto, M. Rosanova, A. G. Casali, J.-F. Brichant, P. Boveroux, et al., Current Biology 25, 3099 (2015).
- McKay et al. (2006) I. D. McKay, L. J. Voss, J. W. Sleigh, J. P. Barnard, and E. K. Johannsen, Anesthesia & Analgesia 102, 91 (2006).
- Vanluchene et al. (2004) A. L. Vanluchene, H. Vereecke, O. Thas, E. P. Mortier, S. L. Shafer, and M. M. Struys, The Journal of the American Society of Anesthesiologists 101, 34 (2004).
- Larson-Prior et al. (2013) L. J. Larson-Prior, R. Oostenveld, S. Della Penna, G. Michalareas, F. Prior, A. Babajani-Feremi, J.-M. Schoffelen, L. Marzetti, F. de Pasquale, F. Di Pompeo, et al., Neuroimage 80, 190 (2013).
- Van Essen et al. (2012) D. C. Van Essen, K. Ugurbil, E. Auerbach, D. Barch, T. Behrens, R. Bucholz, A. Chang, L. Chen, M. Corbetta, S. W. Curtiss, et al., Neuroimage 62, 2222 (2012).
- Kuramoto (1975) Y. Kuramoto, in International symposium on mathematical problems in theoretical physics (Springer, 1975) pp. 420–422.
- Breakspear et al. (2010) M. Breakspear, S. Heitmann, and A. Daffertshofer, Frontiers in human neuroscience 4, 190 (2010).
- Pikovsky et al. (2003) A. Pikovsky, M. Rosenblum, J. Kurths, and J. Kurths, Synchronization: a universal concept in nonlinear sciences, Vol. 12 (Cambridge university press, 2003).
- Silva Pereira et al. (2017) S. Silva Pereira, R. Hindriks, S. Mühlberg, E. Maris, F. van Ede, A. Griffa, P. Hagmann, and G. Deco, Brain Connectivity 7, 541 (2017).
- Biswal et al. (1995) B. Biswal, F. Zerrin Yetkin, V. M. Haughton, and J. S. Hyde, Magnetic resonance in medicine 34, 537 (1995).
- Buckner et al. (2008) R. L. Buckner, J. R. Andrews-Hanna, and D. L. Schacter, Annals of the New York Academy of Sciences 1124, 1 (2008).
- Raichle et al. (2001) M. E. Raichle, A. M. MacLeod, A. Z. Snyder, W. J. Powers, D. A. Gusnard, and G. L. Shulman, Proceedings of the National Academy of Sciences 98, 676 (2001).
- Massimini et al. (2005) M. Massimini, F. Ferrarelli, R. Huber, S. K. Esser, H. Singh, and G. Tononi, Science 309, 2228 (2005).
- Bettinardi et al. (2015) R. G. Bettinardi, N. Tort-Colet, M. Ruiz-Mejias, M. V. Sanchez-Vives, and G. Deco, Neuroimage 114, 185 (2015).
Data and preprocessing
Magnetoencephalography (MEG) data were acquired, preprocessed, anonymized, and distributed by the Human Connectome Project (HCP). Healthy adult human subjects were imaged in a magnetically shielded room using a whole head MAGNES 3600 (4D Neuroimaging, San Diego, CA) system at the Saint Louis University (SLU) medical campus. HCP protocols collected resting state (rest) and visual N-back working memory (active) data from the same subjects. The MEG has 248 magnetometer channels. Electrooculography, electrocardiography, and electromyography recordings are all synchronized with the MEG. Data were acquired with a bandwidth of 400Hz, DC high pass filter, in a continuous acquisition mode, with a 2034.5 Hz sampling rate. MEG data were preprocessed by the HCP. Biological and environmental artifacts, including 60 Hz artifacts, bad channels, and bad segments were removed first automatically by an independent component analysis (ICA)-based, publically available, custom Fieldtrip-scripted pipeline and further quality controlled by human verification. Channel-level preprocessing pipelines also downsampled data to 508.63Hz. N-back data were split into groups of trials corresponding to task design, and divided into epochs time-locked to the onset of the stimulus (TIM) to yield Matlab structures containing epoched, preprocessed data. Further documentation describing data collection and preprocessing is available Larson-Prior et al. (2013) and .
Subjects with data in both N-back memory and resting states (N = 77) were included in the analyses. Per subject average recording time is: 14.060.55 min for resting state and 12.350.67 min for working memory. Respecting constraints introduced by epoch time and Nyquist artifacts, power spectra between [4,100] Hz were obtained with Welch windowing method per , indexing 2 sec MEG time series from each subject , brain state (active or rest), sensor , and epoch .
Spectral energy and entropy
Each power spectrum was used to construct a histogram with bin size 0.1 Hz. Summing over the histogram bins of the power spectra of all sensors, the energy is obtained (Eq. 1). By normalizing each histogram, one obtains the frequency probability distribution for each sensor, epoch, state, and subject. Once again, summing over all bins and all sensors, the entropy can be derived following Eq. 2.
In all our simulations, we consider an all-to-all network of oscillators. The bare frequencies of each oscillator are extracted from a distribution of the following form:
where Hz, Hz, and Hz. We verified that the general picture reported in the main text is not affected by the specific choice of these three parameters. The simulations were run for 15 minutes of simulated time (to match with the experimental data and to ensure a robust estimation of observables), using an Euler integration method with a time step of 5 ms. The energy and entropy are computed on power spectra estimated on windows of 2 s, as with the experiments. A transient time of 100 s was discarded in all simulations.
Phase Lag Index (PLI)
The PLI was computed for each subject and time window, summing over all sensors. The frequency range of interest [4,100] Hz was divided into ten equally wide bins. Within each frequency bin, a time series is obtained by filtering the signal with a Butterworth band-pass filter. The Hilbert transform is then employed to extract the phase of the time series in each frequency bin . From there, the PLI, given by
is computed, where denotes averaging over time Silva Pereira et al. (2017). One may note that the PLI takes values between 0 (random phase relations) and 1 (perfect phase locking). In this work we report the mean PLI over all time epochs for each subject in each brain state.