  A central challenge in the computational modeling of neural dynamics is the
trade-off between accuracy and simplicity. At the level of individual neurons,
nonlinear dynamics are both experimentally established and essential for
neuronal functioning. An implicit assumption has thus formed that an accurate
computational model of whole-brain dynamics must also be highly nonlinear,
whereas linear models may provide a first-order approximation. Here, we provide
a rigorous and data-driven investigation of this hypothesis at the level of
whole-brain blood-oxygen-level-dependent (BOLD) and macroscopic field potential
dynamics by leveraging the theory of system identification. Using functional
MRI (fMRI) and intracranial EEG (iEEG), we model the resting state activity of
700 subjects in the Human Connectome Project (HCP) and 122 subjects from the
Restoring Active Memory (RAM) project using state-of-the-art linear and
nonlinear model families. We assess relative model fit using predictive power,
computational complexity, and the extent of residual dynamics unexplained by
the model. Contrary to our expectations, linear auto-regressive models achieve
the best measures across all three metrics, eliminating the trade-off between
accuracy and simplicity. To understand and explain this linearity, we highlight
four properties of macroscopic neurodynamics which can counteract or mask
microscopic nonlinear dynamics: averaging over space, averaging over time,
observation noise, and limited data samples. Whereas the latter two are
technological limitations and can improve in the future, the former two are
inherent to aggregated macroscopic brain activity. Our results, together with
the unparalleled interpretability of linear models, can greatly facilitate our
understanding of macroscopic neural dynamics and the principled design of
model-based interventions for the treatment of neuropsychiatric disorders.
