The first component of the TVB-AdEx model is at the cellular level, and consists of networks of integrate-and-fire adaptive exponential (AdEx) neurons. As shown in previous studies (Destexhe, 2009; Zerlaut et al., 2018; di Volo et al., 2019), networks of AdEx neurons with adaptation can display asynchronous, irregular (AI) states, as well as synchronous, regular slow-wave dynamics that alternate between periods of high activity (Up) and periods of near silence (Down). The necessary mechanistic ingredients needed to obtain both dynamical regimes include leak conductance and conductance-based synaptic inputs. Each neuron's input is comprised by the firing rates of synaptically connected neurons, weighted by synaptic strengths, as well as stochastic noise (hereafter called “drive”; see Materials and Methods), related biologically to miniature postsynaptic currents. AdEx neurons have the ability to integrate synaptic inputs and fire action potentials, followed by a refractory period (Brette and Gerstner, 2005). AdEx networks with conductance-based synapses can capture features offered by more detailed and computationally expensive models, including AI states and slow-wave dynamics. Figure 1 shows an example of such AI states (Figure 1A) and Up-Down dynamics (Figure 1B) simulated by the same AdEx network, changing only the level of spike-frequency adaptation current (parameter b in the equations, see Material and Methods). In AI states, the firing of individual units remains irregular, but sustained (Figure 1A), whereas in slow-wave states the dynamics alternate between depolarized Up states with asynchronous dynamics and hyperpolarized Down states of near silence (Figure 1B). As such, changes in spike-frequency adaptation lead to differences in cellular kinetics between sleep and wake states. Biologically, spike-frequency adaptation is suppressed by enhanced concentrations of neuromodulators such as acetylcholine during active, conscious brain states that tends to close K⁺ leak channels, resulting in sustained depolarization of neurons (McCormick, 1992) which promotes the emergence of asynchronous, irregular (AI) action potential firing and fluctuation-driven regimes associated with waking states. In contrast, low levels of neuromodulation during unconscious brain states leave leak K⁺ channels open, leading to waves of synchronous depolarisation and hyperpolarization due to the buildup and decay of spike-frequency adaptation, accounting for the emergence of slow-wave dynamics as observed in previous modeling work (Jercog et al., 2017; Nghiem et al., 2020).

The second component of the TVB-AdEx model is a mean-field equation derived from spiking-neuron network simulations, capturing the typical dynamics of a neuron in response to inputs and hence able to describe the mean behavior of a neuronal population, Zerlaut et al. (2018) and di Volo et al. (2019) using a Master-Equation formalism (El Boustani and Destexhe, 2009). This formalism allows one to derive mean fields from conductance-based integrate-and-fire models. It has been shown that—using numerical fits of the transfer function (Zerlaut et al., 2016), an analytical expression for the relationship between a neuron's input and output rates - one can describe complex neuronal models, such as AdEx neurons, and even Hodgkin-Huxley type biophysical models (Carlu et al., 2020). In Figures 1C, D, excitatory and inhibitory firing rates are compared between mean-field simulations using the Master Equation formalism and spiking neural network simulations (time binned population spike counts divided by time bin length T = 0.5 ms). The average adaptation value is also shown for this network (Figures 1C, D, orange curves). These population variables are suitably captured by the mean-field model including adaptation (Capone et al., 2019; di Volo et al., 2019). This mean-field model can exhibit both AI (Figure 1E) and Up-Down dynamics (Figure 1F). Like in the AdEx spiking network model, the transition between two states can be obtained by changing the adaptation parameter called b in both cases (di Volo et al., 2019). With no adaptation, the dynamics are fluctuation-driven around a fixed point exhibiting nonzero firing rates. With adaptation, as the neurons self-inhibit due to adaptation, the nonzero rate fixed point is progressively destabilized by adaptation buildup, driving the dynamics back to the near-zero firing rate fixed point until adaptation wears off and noise drives the system back to the vicinity of the higher-rate fixed point. Thus, with adaptation, the system displays noise-driven alternation between the two fixed-points to generate slow waves (Figure 1G). The regimes are achieved using the mean-field model, which describes excitatory (RS) and inhibitory (FS) population firing rates as well as the mean adaptation level of excitatory populations (Figure 1H).

We have used the simulation engines of the Human Brain Project's (HBP's) EBRAINS neuroscience research infrastructure (https://ebrains.eu and The Virtual Brain https://ebrains.eu/service/the-virtual-brain) to make access as wide as possible. Replication of the TVB-AdEx findings can be done here, with a free EBRAINS account, and users can clone the repositories to further test or extend the present capacities. The models can also be downloaded from Github at https://gitlab.ebrains.eu/kancourt/tvb-adex-showcase3-git to run locally.

Here, the connection of mean-field models was defined by human tractography data (https://zenodo.org/record/4263723, Berlin subjects/QL_20120814) from the Berlin empirical data processing pipeline (Schirner et al., 2015) (Figure 2A). A parcellation of 68 regions was used to place localized mean-field models, with long-range excitatory connections (Figure 2B) and delays (Figure 2C) defined by tract length and weight estimates in human diffusion tensor imaging (DTI) data (Sanz-Leon et al., 2015). Now it becomes possible to simulate brain-scale networks using AdEx-based mean-field models in TVB, hence the name “TVB-AdEx” model.

Having coupled AdEx mean-field models that capture average microscopic characteristics of neural activity, we sought to ascertain if hallmarks of brain-scale (macroscopic) spontaneous activity resembling human brain states were reproduced, as well as whether increases in adaptation strength account for transitions between wake-like and sleep-like macroscopic dynamics. Characterizing temporal hallmarks of simulated neural activity (Figure 3), we find that asynchronous, wake-like dynamics across nodes are recovered in the absence (b = 0 pA, Figure 3A), but not in the presence (b = 60 pA, Figure 3B) of adaptation. Power spectral analysis reveals a peak in the delta range (0.5 − 5 Hz) for the high-adaptation condition (Figure 3C) consistent with empirical data from deeply sleeping individuals. Further, the power spectrum in the low-adaptation condition shows a maximum near 10 Hz (alpha range), consistent with empirically observed dynamics during resting wakefulness (Figure 3C). Therefore, changes in a simulated microscopic process (spike-frequency adaption) influence spectral features of macroscopic brain states, with low-adaptation regimes resembling waking states and high-adaptation regimes reminiscent of slow-wave sleep.

Increasing adaptation can also tune the spatial correlation structure of neural activity across brain states. Indeed, as shown in Figure 4, Pearson correlations across nodes are enhanced in the presence of adaptation, consistent with asynchronous dynamics seen during wakefulness vs. synchronous slow waves seen during deep sleep (Figures 4A, E). We also observe the correlation matrix is structured into two clusters corresponding to the two hemispheres in the slow-wave condition (b = 60 pA, Figure 4E). In addition, increased adaptation strength also causes the emergence of significantly larger correlations between inhibitory than excitatory firing rates across nodes during sleep-like dynamics (Figures 4B, F). This reveals that microscopic variation in adaptation strength alone can account for empirical reports of increased correlations between inhibitory neurons across long distances and even different cortical regions specifically for inhibitory (Peyrache et al., 2012; Le Van Quyen et al., 2016; Olcese et al., 2016). This is due to different effects of adaptation on excitatory regular-spiking neurons and inhibitory fast-spiking neurons, key to reproducing empirical dynamics in unconscious states (Jercog et al., 2017; Nghiem et al., 2020). Moreover, the Phase Lag Index (PLI) is increased during sleep-like dynamics (Figures 4C, G), suggesting systematic phase relations between nodes consistent with traveling slow waves empirically observed during spontaneous unconscious dynamics (Destexhe et al., 1999; Steriade, 2003). Such phase relations, evidenced by a significantly larger PLI, are more pronounced for inhibitory than excitatory neurons in sleep-like dynamics, reminiscent of the key role of inhibitory neurons in organizing the emergence of synchronous dynamics during sleep (Nghiem et al., 2018).

Next, we investigate how the connectome's structure shapes the landscape of synchrony and phase coherence across brain regions, alongside adaptation. In particular, how do the Pearson correlation and PLI scale with spatial distance between nodes? In both b = 0 pA and b = 60 pA conditions (Figures 5A, D), the Pearson correlation between excitatory firing rates significantly decreases with Euclidean distance between regions, corresponding to tract-length-related delays in our model. A steeper negative slope is observed in the awake-like (Figure 5A) than in the slow-wave regime (Figure 5D), suggesting that the spatial extent of synchrony between regions is enhanced in the presence of high-adaptation, sleep-like dynamics.

In Figures 5B, C, E, F, we show a scatter plot and box plot of the Phase-Lag Index (PLI) as a function of distance between regions. In both b = 0 pA and b = 60 pA cases, significant differences are observed in the PLI across regions (Kruskal-Wallis test, ***p < 0.001), suggesting systematic phase relations consistent with the propagation of traveling waves are particularly predominant at intermediate scales. Specifically, in the slow-wave condition (b = 60 pA), we observe that the PLI between regions approximately 65 mm apart is significantly enhanced (Mann-Whitney U test, **p < 0.01) in comparison to the PLI at both shorter and longer distances.

Further, we test whether the predominance of slow waves at an intermediate spatial scale is an emergent property from the structure of the human connectome. To this purpose, we shuffle the connection weights successively for each region and every other region, to generate a connectome with the same distribution of connection weights between regions but retaining none of the graph structure of the empirical connectome (Figure 5G). The tract lengths are not modified. Repeating simulations with shuffled connectomes, we again compute the PLI as a function of distance. With shuffled connectomes, we find that the PLI no longer varies significantly as a function of distance in wake-like dynamics (Figure 5H). As well, the intermediate peak in PLI as a function of distance—denoting an intermediate spatial scale for traveling waves (Figure 5I)—is lost. These results suggest that the non-trivial organization of phase coherence phenomena across cortex is an emergent property of both high adaptation and the weighted graph structure of the human connectome.

Finally, the transition between wake and sleep-like dynamics when changing the level of adaptation is robust for different combinations of parameters of the model (Figure 6). With a high-density scan of the parameter space (see Methods), we find that, for lower values of spike-frequency adaptation (b = 0 pA), AI states are present independently of the timescale (T) of the AdEx mean-field model network nodes and the coupling strength between nodes (S). Consequently, when increasing adaptation, there is a sharp transition from wake-like dynamics to slow-wave activity captured by a marked increase in firing rate standard deviation, again robustly across all values of T and S in the explored range.

After reproducing features of spontaneous dynamics between brain states, we test the hypothesis that changing adaptation in the TVB-AdEx model can also explain differences in empirically observed stimulus-evoked brain responses, with stimuli encoded in more sustained, widespread, reliable, and complex patterns during conscious states (Massimini et al., 2005). To this end, a square wave of 0.1 Hz, matching the magnitude of stochastic drive, with 50 ms duration was input to the firing rates of the transfer function for excitatory populations in the right premotor cortex of the TVB-AdEx simulation during awake-like and slow wave sleep-like conditions, as in previously published empirical studies (Casali et al., 2013).

Figure 7 illustrates the effect of perturbing the large-scale network defined by the TVB-AdEx models. The effect of an external stimulus is apparent for both deep sleep-like and wake-like states (Figure 7A). The average traces of the stimulated region are shown in black, and take into account the 40 realizations shown in gray. To examine the spread of activity following perturbations, the time at which the excitatory firing rate of each region becomes significantly different from the unstimulated baseline (prior to perturbation) is plotted using a color map showing earlier significant changes in brighter colors. Here, we find that responses are more widespread across time and space across brain regions in wake-like (Figure 7B) than sleep-like dynamics (Figure 7C), corresponding to experimental observations in response to Transcranial Magnetic Stimulation (TMS) (Massimini et al., 2005).

To better characterize the complexity of stimulus-evoked responses, the perturbational complexity index (PCI), used in previous experimental works involving TMS, was computed. The PCI is the ratio between Lempel-Ziv complexity, which captures the number of all possible different binary words that can be extracted from binarized responses to stimuli across time and regions, and the entropy of the binarized response that describes how often the response is above the pre-stimulus baseline (see Methods for binarization procedure). A low PCI value indicates a “simple” response to stimulus, while a high PCI value indicates more “complex” response, typically propagating more effectively to different brain areas (Casali et al., 2013). The models were perturbed with stimuli smaller (0.01 Hz) than the stochastic drive, comparable to the noise in amplitude (0.1 Hz), and larger than the drive (1.0 Hz) for simulations in which the value of spike frequency adaptation (b_e) was varied between 0 and 60 pA. As shown in Figures 7D–F, computing the PCI from the TVB-AdEx model shows that PCI values are typically higher for lower-adaptation, wake-like regimes than for higher-adaptation, slow-wave sleep like regimes. In particular, a sharp drop in PCI values is observed between b = 40 pA and b = 60 pA, suggesting an abrupt transition between highly responsive asynchronous and less responsive slow-wave dynamics as adaptation increases. For each value of noise, a one-way ANOVA revealed significant differences between PCI distributions across b-values (p < 0.0001), with multiple comparisons highlighting that the PCI was significantly larger for wake-like than sleep-like conditions, in particular for lower-amplitude stimuli. The same behavior was observed when comparing awake subjects with subjects in slow-wave sleep (Massimini et al., 2005; Casali et al., 2013). One may note a sharp drop of the PCI after b = 40 pA, for which enhanced PCI is observed, reminiscent of a transition from conscious to unconscious response dynamics. As concentrations of neuromodulators such as acetylcholine are known to increase with attention and vigilance, b = 40 pA, which is at the higher-adaptation, lower-neuromodulation end of the spectrum of high-PCI states, could be reminiscent of states that lie between waking vigilance and deep sleep, such as resting wakefulness. Also note the wider distribution of PCI values in sleep-like simulations, suggesting more variable responses for each realization of the same stimulus and therefore less reliable stimulus encoding.

We considered networks of integrate-and-fire neuron models displaying spike-frequency adaptation, based on two previous papers (Destexhe, 2009; Zerlaut et al., 2018). We used the Adaptive Exponential (AdEx) integrate-and-fire model (Brette and Gerstner, 2005). We considered a population of N = 10⁴ neurons randomly connected with a connection probability of p = 5%. We considered excitatory and inhibitory neurons, with 20% inhibitory neurons. The AdEx model permits to define two cell types, “regular-spiking” (RS) excitatory cells, displaying spike-frequency adaptation, and “fast spiking” (FS) inhibitory cells, with no adaptation. The dynamics of these neurons is given by the following equations: (1)cmdvkdt=gL(EL-vk)+gLΔevk-vthrΔ-wk+Isyn (2)uwdwkdt=-wk+b∑tsp(k)δ(t-tsp(k))+a(vk-EL), where c_m = 200 pF is the membrane capacitance, v_k is the voltage of neuron k and, whenever v_k > v_peak = −47.5 mV for inhibitory neurons and v_k > v_peak = −40.0 mV for excitatory at time t_sp(k), v_k is reset to the resting voltage v_reset = −65 mV and fixed to that value for a refractory time T_refr = 5 ms. The voltage threshold v_thr is –50 mV. The leak term g_L had a fixed conductance of g_L = 10 nS and the leakage reversal E_L was of −65 mV for inhibitory and −63 for excitatory. The exponential term had a different strength for RS and FS cells, i.e., Δ = 2 mV (Δ = 0.5 mV) for excitatory (inhibitory) cells. Inhibitory neurons were modeled as fast spiking FS neurons with no adaptation (a = b = 0 for all inhibitory neurons) while excitatory regular spiking RS neurons had a lower level of excitability due to the presence of adaptation (while b varied in our simulations we fixed subthreshold adaptation a = 0 nS and u_w = 500 ms).

The synaptic current I_syn received by neuron i is the result of the spiking activity of all neurons j ∈ pre(i) pre-synaptic to neuron i. This current can be decomposed in the synaptic conductances evoked by excitatory E and inhibitory I pre-synaptic spikes Isyn=Gsyne(Ee-vk)+Gsyni(Ei-vk), Where E_e = 0 mV (E_i = −80 mV) is the excitatory (inhibitory) reversal potential. Excitatory synaptic conductances were modeled by a decaying exponential function that sharply increases by a fixed amount Q_E at each pre-synaptic spike, i.e.,: Gsyne(t)=Qe∑exc.preΘ(t-tspe(k)) e-(t-tspe(k))/ue, Where Θ is the Heaviside function, u_e = u_i = 5 ms is the characteristic decay time of excitatory and inhibitory synaptic conductances, and Q_e = 1.5 nS (Q_i = 5 nS) the excitatory (inhibitory) quantal conductance. Inhibitory synaptic conductances are modeled using the same equation with e → i. This network displays two different states according to the level of adaptation, b = 0 pA for asynchronous irregular states, and b = 60 pA for Up-Down states (see Zerlaut et al., 2018 for details).

We considered a population model of a network of AdEx neurons, using a Master Equation formalism originally developed for balanced networks of integrate-and-fire neurons (El Boustani and Destexhe, 2009). This model was adapted to AdEx networks of RS and FS neurons (Zerlaut et al., 2018), and later modified to include adaptation (di Volo et al., 2019). The latter version is used here, which corresponds to the following equations using Einstein's index summation convention where sum signs are omitted and repeated indices are summed over: (3)T∂νμ∂t=(Fμ-νμ)+12cλη∂2Fμ∂νλ∂νη (4)T∂cλη∂t=δληFλ(1/T-Fη)Nλ+(Fλ-νλ)(Fη-νη)               +∂Fλ∂νμcημ+∂Fη∂νμcλμ-2cλη (5)∂W∂t=-W/uw+bνe+a(μV(νe,νi,W)-EL) , where μ = {e,i} is the population index (excitatory or inhibitory), ν_μ the population firing rate and c_λη the covariance between populations λ and η. W is a population adaptation variable (di Volo et al., 2019). The function F_{μ = {e,i}} = F_{μ = {e,i}}(ν_e, ν_i, W) is the transfer function which describes the firing rate of population μ as a function of excitatory and inhibitory inputs (with rates ν_e and ν_i) and adaptation level W. These functions were estimated previously for RS and FS cells and in the presence of adaptation (di Volo et al., 2019).

At the first order, i.e., neglecting the dynamics of the covariance terms c_λη, this model can be written simply as: (6)Tdνμdt=(Fμ-νμ) , Together with Equation (5). This system is equivalent to the well-known Wilson-Cowan model (Wilson and Cowan, 1972), with the specificity that the functions F need to be obtained according to the specific single neuron model under consideration. These functions were obtained previously for AdEx models of RS and FS cells (Zerlaut et al., 2018; di Volo et al., 2019) and the same are used here.

For a cortical volume modeled as two populations of excitatory and inhibitory neurons, the equations can be written as: (7)Tdνedt=Fe(νe+νaff+νdrive,νi)-νe (8)Tdνidt=Fi(νe+νaff,νi)-νi (9)dWdt=-W/uw+bνe+a(μV(νe,νi,W)-EL), where ν_aff is the afferent thalamic input to the population of excitatory and inhibitory neurons and ν_drive is an external noisy drive simulated by an Ornstein-Uhlenbeck process. The function μ_V is the average membrane potential of the population and is given by μV=μGeEe+μGiE+i+gLEL-WμGe+μGi+gL, where the mean excitatory conductance is μ_Ge = ν_eK_eu_eQ_e and similarly for inhibition.

This system describes the population dynamics of a single region, and was shown to closely match the dynamics of the spiking network (di Volo et al., 2019).

Extending our previous work at the mesoscale (Chemla et al., 2019; di Volo et al., 2019) to model large brain regions, we define networks of mean-field models, representing interconnected brain regions (each described by a mean-field model). We considered interactions between cortical regions as excitatory, while inhibitory connections remain local to each region. The equations of such a network, expanding the two-population mean-field (Equation 7), are given by: (10)Tdνe(k)dt=Fe[νeinput(k)+νaff(k),νi(k)]-νe(k)Tdνi(k)dt=Fi[νeinput(k)+νaff(k),νi(k)]-νi(k) (11)dW(k)dt=-W(k)/uw+bνe(k)               +a(μV(νe(k),νi(k),W(k))-EL) , where ν_e(k) and ν_i(k) are the excitatory and inhibitory population firing rates at site k, respectively, W(k) the level of adapation of the population, and νeinput(k) is the excitatory synaptic input. The latter is given by: (12)νeinput(k)=νdrive(k)+∑jCjk νe(j,t-||j-k||/vc) where the sum runs over all nodes j sending excitatory connections to node k, and C_jk is the strength of the connection from j to k (and is equal to 1 for j = k). Note that ν_e(j, t − ||j − k||/v_c) is the activity of the excitatory population at node k at time t − ||j − k||/v_c to account for the delay of axonal propagation. Here, ||j − k|| is the distance between nodes j and k and v_c is the axonal propagation speed.

The Phase-Lag Index (PLI) was computed for each pair of nodes, averaged over simulation time. The Hilbert transform is employed to extract the phase ψ(t) of the time series. From there, the PLI, given by (13)PLI≡|<sign(ψi(t)-ψj(t))>|, is computed for nodes i and j, where < · > denotes averaging over time (Silva Pereira et al., 2017). One may note that the PLI takes values between 0 (random phase relations or perfect synchrony) and 1 (perfect phase locking). In this work we report the mean PLI over all time epochs for excitatory and inhibitory firing rates of each region pair for each adaptation value.

A model such as the TVB-AdEx contains many parameters whose impact on the dynamics needs to be understood. Additionally, it is necessary to have reasonable, physiological ranges determined for them. As described above, most of them have were already fixed via biological or mathematical arguments, but there is still a subset of parameters whose impact needed to be studied to have a deeper and general understanding of the model. In Table 1, one can find the characteristics of the parameters chosen and the reason for their choice. For each parameter, 16 evenly spaced values were obtained inside the described range. Preliminary results allowed to reduce the explored parameter space to a total of 675,840 different configurations to be analyzed. Using High Performance Computing, the simulation of each parameter combination was parallelized. For each configuration, a seven second simulation was run and, afterwards, multiple features were extracted (mean and standard deviation of the excitatory and inhibitory firing rates, mean value of functional connectivity, etc.). By plotting the value of these features as a function of the parameter values, one can observe the influence of the latter on the model's dynamics (as is shown in Figure 6).

We computed the Perturbational Complexity Index (PCI) in response to a localized square wave stimulus, over the firing rates in a given brain region of a TVB-AdEx simulation, following the method proposed by Casali et al. (2013). This stimulus was applied by augmenting the firing rate of the excitatory population during the pulse. This is done for multiple trials with the same stimulus delivered to the same node at a random point in time and with different realizations of noise. The PCI is the ratio of two quantities: the Lempel-Ziv algorithmic complexity and the source entropy (Casali et al., 2013). To compute both quantities, firing rates ν(t) must be binarized to produce significance vectors s(t). First, the trials are aligned to stimulation time, considering only the 300 ms before and after stimulus onset. Then, each node's firing rate is re-scaled and mean and standard deviation given by pre-stimulus activity averaged over nodes. Afterwards, all pre-stimulus firing rates are randomized across time bins, this procedure being repeated 500 times. The threshold for significance T is then given by the one-tail percentile of the maximum absolute value over all repetitions within a series of 20 trials. For each trial of those 20 trials, we can then write s(t) = 1 whenever post-stimulus ν(t) > T and s(t) = 0 otherwise. For what follows, we concatenate all s(t) vectors from all simulation nodes into one single significance vector S(t) per trial.

The Lempel-Ziv complexity LZ(S) is the length of the “zipped” vector S(t), i.e., the number of possible binary “words” that make up the binary vector S(t). Briefly, S(t) is sectioned successively into consecutive words of between one and N_t characters where N_t is the total length of S(t). Scanning sequentially through all words, each new encountered word is added to a “dictionary,” and LZ(S) is the total number of words in the dictionary at the end of the procedure.

The spatial source entropy H(S) is given by: (14)H(S)=-p(S=0)log2(p(S=0))-p(S=1)log2(p(S=1)), where log₂ denotes the base-two logarithm.

The PCI can then be expressed as PCI(S)=LZ(S)H(S).