Bridging multiple scales in the human brain using computational modelling

Brain dynamics span multiple spatial and temporal scales, from fast spiking neurons to slow fluctuations over distributed areas. No single experimental method links data across scales. Here, we bridge this gap using The Virtual Brain connectome-based modelling platform to integrate multimodal data with biophysical models and support neurophysiological inference. Simulated cell populations were linked with subject-specific white-matter connectivity estimates and driven by electroencephalography-derived electric source activity. The models were fit to subject-specific resting-state functional magnetic resonance imaging data, and overfitting was excluded using 5-fold cross-validation. Further evaluation of the models show how balancing excitation with feedback inhibition generates an inverse relationship between α-rhythms and population firing on a faster time scale and resting-state network oscillations on a slower time scale. Lastly, large-scale interactions in the model lead to the emergence of scale-free power-law spectra. Our novel findings underscore the integrative role for computational modelling to complement empirical studies.

RSNs. Several studies report negative correlations between cortical fMRI blood oxygen level-dependent (BOLD) contrast signals and EEG-derived α-rhythm power 11,12,31,32 . EEG α-rhythms have traditionally been interpreted as a sign of cortical "idling" 33 . This notion is now extended by a growing body of research that emphasizes their role for information processing 34 . Importantly, attention, perceptual awareness and cognitive performance are rhythmically modulated by α-power and phase [35][36][37][38] . The roles of α-rhythms for mediating topdown control, timing of oscillations and directing attention by blocking task-irrelevant pathways are central to prevailing hypotheses termed 'gating by inhibition' and 'pulsed inhibition' 14,15 . Intracellular recordings showed that inhibition is inseparable from excitatory events, resulting in ongoing excitation-inhibition balance (E/I balance) 39,40 . Finally, despite widespread interest in criticality and its ubiquity in nature 18 the origin of power-law scaling observed in both fMRI and EEG data is unclear 16,17 . A key point of consideration across the three phenomena we describe in the context of the present paper is that the computational model was not explicitly constructed to address each of them. Rather we used an empirically constrained connectome-based brain network model to provide plausible multiscale explanations that link these phenomena and uncover underlying mechanisms. Population models are injected with region-wise EEG source activity approximating locally generated synaptic input currents 9,10 . This enables person-specific simulations. Simulated whole-brain fMRI time series are fitted with empirical fMRI time series that were acquired simultaneously with EEG. 5-fold cross-validation was performed to guard against overfitting.

Prediction of individual fMRI time series
The underlying simulated neural activity (firing rates and synaptic activity) is analysed to reveal links between mesoscopic dynamics and neuroimaging signals. Excitatory postsynaptic currents (EPSCs) tend to be strongly correlated with the electric fields in their vicinity, while inhibitory postsynaptic currents (IPSCs) are anticorrelated with EPSCs 8,39,40 . Here, EPSCs are derived from EEG and IPSCs result from modelled inhibitory population activity.

Synaptic gating (the fraction of open channels of neural populations) was transformed to
region-wise fMRI signals using the Balloon-Windkessel forward model 41 . In contrast to noisedriven models that only capture stationary features, hybrid models as used here are able to predict the ongoing temporal dynamics of fMRI time series.
Using exhaustive searches, we fitted parameters for brain network models of 15 human adult subjects, who all had T1-weighted, diffusion-weighted and simultaneous resting state EEG-fMRI data (Fig. 1). For each brain network model, three parameters were varied to maximize the fit between empirical and simulated fMRI: the scaling of excitatory white-matter coupling and the strengths of the inputs injected into excitatory and inhibitory populations. Besides tuning these three global parameters to maximize fMRI fit, we also tuned all local inhibitory coupling strengths in order to obtain biologically plausible firing rates in excitatory populations. For this second form of tuning, termed feedback inhibition control (FIC), average population firing rates were the sole optimization criterion, without any consideration of prediction quality, which was only dependent on the three global parameters (see Methods). FIC is required to compensate for excess or lack of excitation resulting from the large variability in white-matter coupling strengths obtained by MRI tractography, which is a prerequisite to obtain plausible ranges of population activity that is relevant for some of the following results (Figs. 4 and 5). Prediction quality was measured as the average correlation coefficient between all simulated and empirical region-wise fMRI time series of a complete cortical parcellation over 20.7 minutes length (TR = 1.94s, 640 data points) thereby quantifying the ability of the model to predict the activity of 68 parcellated cortical regions.
Accounting for the large-scale nature of fMRI resting-state networks, the chosen parcellation size provides a parsimonious trade-off between model complexity and the desired level of explanation. What this parcellation may lack in spatial detail, it gains in providing a full-brain coverage that can reliably reproduce ubiquitous large-scale features of empirical data, which we further present below. To exclude overfitting and limited generalizability, we performed 5-fold cross-validation. For each subset (fold), the fMRI data were randomly divided into an 80% sample as training segment and a 20% sample as testing segment. Furthermore, despite that large range of possible parameters, the search converged to a global maximum ( Supplementary Fig. 1). Therefore, we ensured that when the model has been fit to a subset of empirical data, that it was able to generalize to new or unseen data. In contrast to model selection approaches, where the predictive power of different models and their complexity are compared against each other, we here use only a single type of model.
The hybrid model simulation results were compared with three control scenarios: (i) a model where local dynamics were driven by noise, (ii) a variant of the hybrid model that used random permutations of the region-wise EEG source activity time series and (iii) a statistical model where ongoing α-band power fluctuation of EEG source activity was convoluted with the canonical hemodynamic response function (henceforth called α-regressor). The first two controls are brain network models; the third is inspired by traditional analyses of empirical EEG-fMRI data. The combination of 5-fold cross validation and comparison against these control scenarios further reinforces the explanatory utility of the hybrid computational model.  Visual inspection of example time series showed good reproduction of characteristic slow RSN oscillations by the hybrid model and the α-regressor (albeit inverted for the latter), but poor reproduction of temporal dynamics in the case of noise and random permutations models (Fig. 2a). To quantify prediction quality, we computed the mean over all correlation coefficients between simulated and empirical fMRI time series for each parameter set and each of the four setups (i.e., hybrid model and the three control setups). Predictions from the hybrid model correlated significantly better with empirical fMRI time series than predictions from the two random models (Fig. 2b). Five-fold cross-validation showed no significant difference of prediction quality between training and validation data sets (p = .97, two-sample t-test) and between validation data sets and prediction quality for the full time series (p = .34, two-sample t-test). Predictions from the α-regressor were slightly less accurate (note that as α-power is inversely correlated with fMRI activity, correlation coefficients were inverted for comparison). We performed a group-level spatial independent component analysis (ICA) on the empirical fMRI data to compare simulation results with RSN activity (Fig. 2d). To account for temporal variation of prediction quality, we computed sliding-window correlations (100 fMRI scans = 194 s window length; one fMRI scan step width). Average prediction quality was largest for the hybrid model for eight out of nine RSNs (Fig. 2d). As with region-wise fMRI, correlation coefficients of the hybrid model and the α-regressor were comparably good, while the control network models performed poorly. We found that prediction quality varied over time, regions and subjects. Window-wise prediction quality was highly correlated with the standard deviation of RSN temporal modes (Fig. 2c). That is, the higher the temporal activation of a RSN (in a given subject, region and time window) the better the prediction of its activation time series. As a consequence, epochs in the upper quartile of RSN standard deviations were significantly better predicted than epochs in the lower quartile (Fig. 2d).
Similarly, the strength of functional networks (measured by average FC between all regionpairs) correlated with prediction quality (r=0.5 ± 0.32, mean ± s.d.), indicating that epochs of high functional network activity were better predicted than epochs that are characterized by uncorrelated activity.  (Fig. 2, 4-6 In contrast to time series prediction, the α-regressor performed poorly for estimating FC. Compared to the α-regressor, all three model-based approaches provided significantly better predictions of subject-specific long-term FC (computed over ~20 min between all 68 ROIs) and functional connectivity dynamics (Fig. 3). Furthermore, FC matrices obtained from hybrid model simulations correlated significantly higher with empirical activity than predictions obtained from the noise-driven model (Fig. 3a, b). Interestingly, correlations for hybrid and random permutation models were effectively the same, likely because the largescale network dynamics that drive the emergence of FC would be relatively preserved when permuting injected activity. Prediction of group-average FC (all pairwise FC values averaged over all subjects) was higher for the hybrid model compared to the α-regressor (Fig. 3c).

E/I balance generates α-phase related firing
After fitting the 15 person-specific hybrid models, we analysed the local population activity underlying predicted fMRI time series to infer the neurodynamic mechanisms. Simulated dynamics result from fitting three global model parameters: the scaling of global white-matter coupling and the strengths of inputs injected into excitatory and inhibitory populations. The optimized hybrid model reproduced the inverse relationship between α-phase and firing rates observed in invasive recordings 42 (Fig. 4a). To investigate these fast-acting dynamics related to α-phase, we computed grand average waveforms of modelled synaptic inputs, population firing rates and synaptic gating, all time-locked to the zero-crossings of α-cycles (Fig. 4b).
Resulting waveforms illustrate the relation of α-oscillations to population activity and how the ongoing balancing of rhythmic excitatory and inhibitory inputs generated pulsed inhibition.   8,43 . Accordingly, excitatory population input, firing rates and synaptic gating (red; columns IV, VI, VII) were inverted relative to the α-cycle.
As noted above, we injected region-wise EEG source activity into corresponding network nodes in order to simulate biologically realistic EPSCs (Fig. 1). In hybrid models, which were optimized for highest empirical fMRI time series fit and biological plausibility of average firing rates, EPSCs dominated inputs to inhibitory populations. Consequently, the sum of synaptic input currents to inhibitory populations closely followed the shape of EPSCs.
As a result of the monotonic relationship between input currents and output firing rates of populations (defined by Eqs. 3 and 4, see Methods), the waveform of inhibitory firing rates and synaptic gating also closely followed injected EPSCs. As increasing input to inhibitory populations leads to increasing inhibitory effect and vice versa, resulting feedback inhibition (i.e. IPSC) waveforms were inverted to injected EPSCs; that is, excitation and inhibition were balanced during each cycle, which is in accordance with published electrophysiology results 8,40 . Consequently, IPSCs peaked during the trough of the α-phase and were lowest during the peak of the α-phase. The previously fitted synaptic coupling strength parameters determined the relative contribution of excitation and inhibition. Fitting the models to fMRI activity resulted in a biologically plausible ratio 8,43 of EPSCs to IPSCs, with IPSC amplitudes being about three times larger than EPSC amplitudes (Fig. 4b). As per other empirical observations 8, 43 , we ensured the EPSC amplitudes at excitatory populations are smaller than IPSC amplitudes by constraining parameters such that the standard deviation of injected source activity was larger for inhibitory populations, i.e., ω BG (E) < ω BG (I) . Specifically, we scanned 12 different ratios of both parameters ω BG (I) / ω BG (E) with values between 5 and 200.
Because IPSCs have dominated excitatory population inputs, firing rates of excitatory populations showed a similar shape as feedback inhibition currents, i.e., they peaked during the trough of the α-cycle and fell to their minimum during the peak of the α-cycle, reproducing their empirical relationship 42 .
In summary, the mesoscopic activity underlying fMRI predictions showed a rhythmic modulation of firing rates on the fast time scale of individual α-cycles in accordance with empirical observations 42 . Model activity revealed that periodic microstates of excitation and inhibition resulted from the ongoing balancing of EPSCs by feedback inhibition currents.

'Pulsed inhibition' generates fMRI oscillations
We next focused on the interaction between α-power fluctuations and population dynamics to identify the mechanism underlying the inverse relationship between α-power and fMRI 12 . To isolate the effects of α-waves and separate it from other EEG rhythms, we replaced the EEGsource activity in the optimized model for each subject with artificial α-activity consisting of a 10 Hz sine wave that contained a single brief high power burst in its centre (using the single parameter set for which the highest average prediction quality over all subjects was obtained, Supplementary Fig. 1) and computed grand average waveforms over all resulting region time series (Fig. 5a). As for fMRI prediction, local inhibition strengths for each node were tuned using FIC to obtain biologically realistic firing rates. We analysed waveforms of model state variables and found, as expected, that average per-cycle firing rates and fMRI activity decreased in response to the α-burst (Fig. 5a).
Notably, this behaviour emerged despite the fact that injected activity was centred at zero, i.e., positive and negative modulations of input currents were balanced. The reason for the observed asymmetric response to increasing power levels originated from inhibitory population dynamics: while positive deflections of α-cycles generated large peaks in ongoing firing rates of inhibitory populations, negative deflections were bounded by 0 Hz (Fig. 5a).
Because of this rectification of high-amplitude negative half-cycles, average per-cycle firing rates of inhibitory populations increased with increasing α-power. As a result, also feedback inhibition had increased for increasing α-power, which in turn led to increased inhibition of excitatory populations, decreased average firing rates, synaptic gating variables and ultimately fMRI amplitudes.
We next analysed the relationship between fMRI oscillations and α-power. We generated artificial α-activity consisting of a 10 Hz sine wave that was amplitude modulated by slow oscillations (cycle frequencies between 0.01 and 0.03 Hz) and injected it into the models of all subjects using the same parameters as for the single α-burst (Fig. 5b). In summary, we found that increased α-power led to an increased activation of inhibitory populations. Resulting feedback inhibition generates fMRI oscillations, which can explain the empirically observed anticorrelation between α-power and fMRI.

Global coupling controls scale-freeness
Empirical fMRI power spectra follow a power-law distribution P ∝ f β , where P is power, f is frequency and β the power-law exponent, which is an indicator of criticality and scale-free dynamics 44 . In accordance with systematic analyses 44 , average power-law spectra of our fMRI data obeyed power-law distributions with exponents β emp = -0.97 (β emp = -0.82 when the frequency range was limited to frequencies < 0.1 Hz as in Ref [44]) for empirical and β sim = -0.76 for simulated data ( Fig. 6a and Supplementary Fig. 2); time series were tested for scale invariance using rigorous model selection criteria (see Methods). Our previous results (Fig. 5) associated resting-state fMRI oscillations with EEG by identifying a neural mechanism that transforms instantaneous EEG source power fluctuations into fMRI oscillations. Surprisingly, however, the spectrum of EEG source power had a considerably smaller negative exponent than fMRI (β α-band = -0.56 for α-power and β wide-band = -0.53 for wide-band power).
Comparison of power spectra indicated that power-law fits of simulated fMRI spectra had a higher negative exponent than power-law fits of instantaneous source power spectra because the power of slower oscillations increased relative to the power of faster oscillations ( Fig. 6a and Supplementary Fig. 2). Model dynamics transform input activity such that the amplitude of output oscillations increased inversely proportional to their frequency. The effect is visible in Fig. 5b, where fMRI, synaptic and firing rate amplitudes of slow oscillations were larger than amplitudes of fast oscillations, despite equally large amplitudes of input α-power oscillations. In comparison, simulation results obtained for deactivated large-scale coupling, but an otherwise identical setup, did not show this frequency-dependent amplification (Fig.   6b). Without large-scale coupling the power-law exponent of simulated fMRI (β sim_Gzero = -0.53) was close to the exponent of injected EEG source activity power (β wide-band = -0.53). We  Fig. 2 and 3). Screening of subject-specific parameter spaces showed that the power-law exponent of simulated fMRI depended on the balance of large-scale excitation and local inhibition: the 2D distribution of prediction quality (fMRI time series and functional connectivity) and critical exponent showed a characteristic diagonal pattern, which demonstrates the necessity of proper balance of large-scale excitation and local inhibition for optimal prediction (Supplementary Fig. 3).  Figs. 2 and 3). When large-scale coupling was absent or inadequately balanced, prediction quality decreased and models produced flatter spectra or no scale invariance at all. b, Grand average waveforms of simulation results using the same artificial α-rhythm input as in Fig. 5b

DISCUSSION
Connectome-based computational modelling has helped to provide new insights into largescale network dynamics related to phenomena such as resting-state networks 45 . We build from this foundation by more explicitly integrating empirical data to constrain the model dynamics.
Subject-specific connectomes and biophysically based brain models were driven by their own EEG-derived approximations of locally generated synaptic input currents and fit to their own fMRI time series. The optimized models enabled inferences about mesoscopic dynamics underlying fMRI resting state activity. We reinforce the validity of the model by showing that the underlying emerging processes, which were not explicitly incorporated in the model, reproduce neurophysiological results. The novel approach forms an integrative framework to construct hypotheses derived from the generated model activity, while remaining sensitive to constraints provided by empirical evidence. Thereby, it links data to theory by uniting models of neural dynamics with empirical observations across modalities and spatiotemporal scales.
We demonstrate how this framework can be used for the systematic study of individual brain activity by identifying neural mechanisms underlying person-specific resting-state fMRI activity, the inhibitory effect of EEG α-rhythms and the emergence of scale-free dynamics.
Power-law relationships between frequency and amplitude of oscillations are omnipresent in nature; the identified mechanism-increased amplification of slower processes through prolonged recurrent excitation-may have implications for other areas of science. The observed co-emergence of long-range correlations and power-law scaling may point to a unifying explanation within the theory of self-organized criticality, as previously proposed by others 46 .
Fitting models with subject-specific fMRI time series configured network interaction such that mesoscopic dynamics emerged that are consistent with other empirical data at that scale.
Estimated parameters, like the strengths of local and global coupling, created dynamical regimes that were characterized by ongoing E/I balance, which we also found to be a key mechanism underlying fMRI resting-state oscillations and α-inhibition. Our findings mechanistically link electrical activity to hemodynamic oscillations and thereby add to accumulating evidence suggesting that RSNs originate from neuronal activity 11-13, 25, 26, 29, 30 rather than being a purely hemodynamic phenomenon that is only correlated, but not caused by it 32,47,48 . Likewise, our results provide a mechanistic explanation for the emergence of power-law spectra that is based on large-scale interaction of brain regions. Our identified mechanism suggests a neural origin in line with empirical results that connect E/I balance to scale-free dynamics 49 and recurrent excitation to states of high or low activity in local circuits 50 . Our major results are integrated and visualized in Supplementary Video 1.
Source activity injection aims to increase the biological realism of brain simulation.
Computational models by necessity omit features of the real system for the sake of simplicity, generality and efficiency. Adding such features increases degrees of freedom, but may also render parameter spaces intractable and increases the risk of over-fitting. Injection of empirical activity is a way to systematically probe sufficiently abstract neural systems while maintaining biologically realistic behaviour. Thereby, the approach aims to balance a level of abstraction that is sufficient to provide insights with being detailed enough to guide subsequent empirical study.
The present approach proved promising for general decoding of mesoscopic neural information processing mechanisms. The ability to selectively target and perturb personspecific whole-brain models will contribute to our mechanistic understanding of neural information processing and its relation to brain function and dysfunction.   (7) where J i k denotes the value of feedback inhibition strength of node i and τ k denotes the adaptive tuning factor during the k-th iteration. In the first iteration, all J i values are initialized with 1 and τ k is initialized with 0.005. The adaptive tuning factor is dynamically changed during each iteration based on the result of the previous iteration:

METHODS
For the case that the result did not improve during the current iteration, i.e., After preprocessing, the cortical parcellation mask was registered to fMRI resting-state data (T2-weighted echo planar imaging with blood oxygen level-dependent contrast; TR = 1940 ms; voxel size = 3 mm isotropic; eyes-closed resting-state) of subjects and average fMRI signals for each region were extracted. The first five images of each scanning run were discarded to allow the MRI signal to reach steady state. To identify RSN activity a spatial Group ICA decomposition was performed for the fMRI data of all subjects using FSL MELODIC (MELODIC v4.0; FMRIB Oxford University, UK, Beckmann and Smith (61)) with the following parameters: high pass filter cut off: 100 s, MCFLIRT motion correction, BET brain extraction, spatial smoothing 5 mm FWHM, normalization to MNI152, temporal concatenation, dimensionality restriction to 30 output components. ICs that correspond to RSNs were automatically identified by spatial correlation with the nine out of the ten wellmatched pairs of networks of the 29,671-subject BrainMap activation database as described in Smith et al. 62 (excluding the cerebellum network).
EEG preprocessing. Details of EEG preprocessing are described in supplementary material of Schirner et al. 57 . First, to account for slow drifts in EEG channels all channels were high-pass filtered at 1.0 Hz (standard FIR filter). Imaging Acquisition Artefact (IAA) correction was performed using Analyser 2.0 (v2.0.2.5859, Brain Products, Gilching, Germany). The onset of each MRI scan interval was detected using a gradient trigger level of 300 µV/ms.
Incorrectly detected markers, e.g. due to shimming events or heavy movement, were manually rejected. To assure the correct detection of the resulting scan start markers each inter-scan interval was controlled for its precise length of 1940 ms (TR). For each channel a template of the IAA was computed using a sliding average approach (window length: 11 intervals) and subsequently subtracted from each scan interval. For further processing, the data was down   illustrate that when large-scale coupling was absent or inadequately balanced, models did not produce scale-free behaviour and prediction quality decreased.