## Abstract

Brain function depends on segregation and integration of information processing in brain networks often separated by long-range anatomic connections. Neuronal oscillations orchestrate such distributed processing through transient amplitude and phase coupling, yet surprisingly, little is known about local network properties facilitating these functional connections. Here, we test whether criticality, a dynamical state characterized by scale-free oscillations, optimizes the capacity of neuronal networks to couple through amplitude or phase, and transfer information. We coupled *in silico* networks which exhibit oscillations in the α band (8–16 Hz), and varied excitatory and inhibitory connectivity. We found that phase coupling of oscillations emerges at criticality, and that amplitude coupling, as well as information transfer, are maximal when networks are critical. Importantly, regulating criticality through modulation of synaptic gain showed that critical dynamics, as opposed to a static ratio of excitatory and inhibitory connections, optimize network coupling and information transfer. Our data support the idea that criticality is important for local and global information processing and may help explain why brain disorders characterized by local alterations in criticality also exhibit impaired long-range synchrony, even before degeneration of axonal connections.

**SIGNIFICANCE STATEMENT** To perform adaptively in a changing environment, our brains dynamically coordinate activity across distant areas. Empirical evidence suggests that long-range amplitude and phase coupling of oscillations are systems-level mechanisms enabling transient formation of spatially distributed functional networks on the backbone of a relatively static structural connectome. However, surprisingly little is known about the local network properties that optimize coupling and information transfer. Here, we show that criticality, a dynamical state characterized by scale-free oscillations and a hallmark of neuronal network activity, optimizes the capacity of neuronal networks to couple through amplitude or phase and transfer information.

## Introduction

Segregation of function across cortical areas allows for economical wiring and reduces metabolic costs (Achard and Bullmore, 2007; Sporns, 2010; Wang and Clandinin, 2016). However, to function adaptively in a dynamic environment, the brain also has to flexibly integrate its activity across distant cortical regions. Neuronal oscillations are an important network-level mechanism driving coordination not only in local cortical assemblies but also through amplitude and phase coupling of distant brain regions (Varela et al., 2001; Fries, 2005). In fact, mapping functional connectivity in the brain often relies on quantifying such phase (Sauseng et al., 2005; Stam et al., 2007; Palva et al., 2010; Crespo-Garcia et al., 2013; Marzetti et al., 2019; Tewarie et al., 2019) and amplitude dependencies (Daffertshofer and van Wijk, 2011; Hipp et al., 2012; Vidal et al., 2012; Foster et al., 2015; O'Neill et al., 2015; Tewarie et al., 2019). However, although coupling of oscillations enables information exchange, a fully synchronized brain is also not advantageous, since such a stable state does not allow for creation of new information (Tagliazucchi and Chialvo, 2013; Tognoli and Kelso, 2014). For this reason, several investigators have argued that the brain resides optimally in a state that is neither fully segregated nor fully integrated, where transient changes in phase and amplitude coupling allow it to explore the largest repertoire of functional states supported by the static structural connectome (Tagliazucchi and Chialvo, 2013; Deco and Kringelbach, 2016; Li et al., 2019). The so-called “critical state” has often been suggested as a statistically well-defined dynamical state that could fulfill the contrasting demands of stability and susceptibility to inputs (Bak et al., 1987; Chialvo, 2010). Indeed, the hallmark of the critical state is spatiotemporal dynamics with power-law scaling behavior (Linkenkaer-Hansen et al., 2001; Beggs and Plenz, 2003), indicating that spatial and temporal patterns on a wide range of scales can occur.

In a computational model of critical oscillations (CROS; Poil et al., 2012; Dalla Porta and Copelli, 2019; Avramiea et al., 2020), it was shown that for any level of inhibitory connectivity, scale-free fluctuations in the amplitude of oscillations appeared when a suitable amount of excitatory connections were added to the network, indicating that critical dynamics of oscillations requires a specific balance of excitation and inhibition. Interestingly, at the same balance we find scale-free distributions of neuronal avalanches, another hallmark of critical brain dynamics (Beggs and Plenz, 2003). Both neuronal avalanches (Kinouchi and Copelli, 2006; Shew et al., 2009; Gautam et al., 2015) and critical oscillations (Avramiea et al., 2020) have been associated with the widest dynamic range to stimuli. This suggests that at the critical point, networks might be most responsive to inputs coming from distant brain areas. Consequently, alterations in the local criticality of networks might lead to impaired coupling. In agreement with this idea, brain disorders such as Alzheimer's (Montez et al., 2009) and schizophrenia (Nikulin et al., 2012) have been associated with alterations in criticality, excitation/inhibition (E/I) balance (Lisman, 2012; Busche and Konnerth, 2016), and reduced long-range phase coupling in fronto-parietal networks (Babiloni et al., 2004; Sheffield et al., 2015). Of note, there is considerable variation in critical brain dynamics among individuals (Linkenkaer-Hansen et al., 2007) and brain states such as focused attention (Irrmischer et al., 2018), suggesting that it is biologically relevant to understand the implications of near-critical brain dynamics and E/I balance on long-range coupling and information processing.

Here, we hypothesize that local criticality of neuronal networks is required for establishing long-range coupling in the amplitude and phase of neuronal network oscillations, as well as for transferring information. We tested this hypothesis in the critical oscillations model (Poil et al., 2012; Dalla Porta and Copelli, 2019; Avramiea et al., 2020). We connected two networks and varied their E/I ratio and, consequently, the level of criticality, to understand its impact on long-range coupling and information transfer. We found that information transfer is optimized in the critical regime, where networks show maximal coupling in the phase and amplitude of oscillations. Thus, regulation of criticality may be important for long-range communication, and alterations in criticality may explain disruptions in long-range synchrony in brain disorders, along with their functional consequences.

## Materials and Methods

#### CROS model

CROS is based on a probabilistic firing-rate model (Dayan and Abbott, 2001), which exhibits power-law avalanches in a network of excitatory neurons (Abbott and Rohrkemper, 2007). Poil et al. (2012) extended the model in Abbott and Rohrkemper (2007) to include inhibitory neurons, and showed that power-law avalanches as well as long-range temporal correlations (LRTCs) in the amplitude of oscillations co-occur at a balance between excitatory and inhibitory connectivity. Here, we used the model described in (Poil et al., 2012), with synaptic weights optimized for power-law avalanches and LRTCs. The model was implemented using the Brian2 simulator for spiking neural networks (Stimberg et al., 2014; RRID:SCR_002998). Specifically, we modeled networks of 75% excitatory and 25% inhibitory integrate-and-fire neurons arranged on a 50 × 50 open grid. Placing neurons on the grid using uniform sampling may result, by chance, in clusters of excitatory or inhibitory neurons. To avoid this, we first positioned inhibitory neurons on the grid according to Mitchell's best candidate algorithm (Mitchell, 1991), which results in a more even spatial distribution. The remaining spaces were filled with excitatory neurons. Networks differ in their two connectivity parameters, *P*, of a connection at a distance *r* was given by:
*i*, with connectivity probability

As such, we used the Nelder–Mead optimization algorithm to determine the value of α which minimizes the following function:

#### Neuron model

Neurons were modeled using a synaptic model integrating received spikes, and a probabilistic spiking model. Each time step (*dt*) of 1 ms starts with each neuron, *i*, updating the input

the weights *S* is a is a binary vector, with *j* fired in the previous time step, and

The activation of a neuron

The spiking probability *A* at the current timestep, as follows:

We determine whether the neuron spikes with the probability *A* is reset to the reset value,

#### Model parameters

All model parameters were the same as in the original paper (Poil et al., 2012) apart from the synaptic weights. Neuron model: (*A*_{0} = 0.000001, *A*_{r} = −2); inhibitory neurons (*A*_{0} =0, *A*_{r} = −20). To improve the range and stability of the LRTCs from the original model, an evolutionary algorithm (Smit and Eiben, 2011) was applied to the synaptic weights. The parameters that could vary were the two connectivity parameters (taking values between 0% and 100%), and the natural logarithm of the magnitude of the four synaptic weights, *W _{EE}*,

*W*,

_{IE}*W*and

_{EI}*W*(taking values between −5 and 1). For each run, a fitness value was calculated based on the avalanches size and duration distributions (see below, Neuronal avalanches), and the LRTCs [see below, Detrended fluctuation analysis (DFA) of LRTCs].

_{II}The optimum weights (*W _{ij,}* connecting the presynaptic neuron

*j*to the postsynaptic neuron

*i*) found by the algorithm were (

*W*= 0.0085,

_{EE}*W*= 0.0085,

_{IE}*W*= −0.569,

_{EI}*W*= −2). Importantly, the relationship between multi-level criticality and E/I balance which emerges in Figure 1

_{II}*F*,

*G*was not part of the optimization and, thus, is a generic feature of the model. More specifically, although the evolutionary algorithm tuned a combination of excitatory and inhibitory connectivity that produces critical dynamics, one can find for any other excitatory connectivity also an inhibitory connectivity that balances the excitation and produces critical dynamics.

The initial activation *A*_{0}.

#### Long-range connections

To assess how reciprocal long-range connections affect synchrony of oscillations and information transfer in our model, we generated two CROS model grids of 50 × 50 neurons. The placement of excitatory and inhibitory neurons on the two grids was exactly the same, to allow connection of neurons of the same type and at the same positions in the two networks. However, the local connectivity was allowed to vary between the two grids, according to the probabilistic connectivity rules described above. To create *N* long-range connections between the two networks (where *N* < number of excitatory neurons on the grid = 1875), we sampled *N* excitatory positions. Half of these (*N*/2) are long-range connections from the first to the second network, projecting to excitatory neurons at the same positions on the grid. The other half connect from the second to the first network. This procedure ensures that there are no connections that go both ways.

#### Conduction delay

Zero milliseconds for local connections within each grid; 25 ms on the long-range connections between the two grids, unless otherwise specified. The value of 25 ms for conduction delay was chosen because it allows for both amplitude coupling (Fig. 2*H*) and phase coupling (Fig. 3*G*) to occur between the two networks.

#### Network input

A Gaussian white noise signal S was generated, with a length and sampling frequency equal to that of the network signal. The signal was then convoluted with an exponential:

#### Network activity analysis

A network signal was created by summing the total number of neurons spiking at each time-step with a Gaussian white noise signal of the same length with mean = 0 and σ = 3. This level of white noise was set to allow all networks to achieve a time-varying phase, which is not the case without adding the noise, when there are silent periods in the network. Amplitude envelopes were extracted from this preprocessed signal, with Gaussian noise added. The analysis of neuronal avalanches was performed on the raw signal, consisting of the total number of neurons spiking at each time step.

#### Oscillation power

The power spectrum was computed using the Welch method with a Hamming window of size with 2^{11} fast Fourier transform (FFT) points.

#### Neuronal avalanches

A neuronal avalanche is defined as a period where neurons are spiking above an activity-dependent threshold (Beggs and Plenz, 2003; Shew et al., 2011). Networks in a strongly excitation-dominated regime exhibit relentless activity throughout the entire duration of the simulation, if the avalanche threshold was set to 0 for these networks, only one avalanche would be identified. To allow computation of avalanche statistics also for these networks, the avalanche threshold was set to be dependent on the overall activity of each network: in our case, it was set to half the median of activity. The size of the avalanche is the number of spikes during this period. We then computed the κ index (Shew et al., 2009; Poil et al., 2012), which calculates the difference between the distribution of our data and a power-law, by calculating the average difference of the cumulative distribution of a power-law function, *P* (with exponent −1.5 for size and −2.0 for duration; Zapperi et al., 1995) and that of our experimental data, *A*, at 10 equally spaced points on a logarithmic axis (β) and adding 1:

A subcritical distribution is characterized by κ < 1, and a supercritical distribution by κ > 1, whereas κ = 1 indicates a critical network.

#### Neuronal oscillations

Data were filtered using a one-way causal FIR filter (order 250) between 8 and 16 Hz, which was the most prominent oscillatory regime exhibited in our model (Poil et al., 2012). Oscillatory amplitude was computed as the absolute of the Hilbert transform.

#### Detrended fluctuation analysis (DFA) of LRTCs

DFA was used to analyze the scale-free decay of temporal (auto)correlations, also known as LRTCs. The DFA was introduced as a method to quantify correlations in complex data with less strict assumptions about the stationarity of the signal than the classical autocorrelation function or power spectral density (Peng et al., 1995; Linkenkaer-Hansen et al., 2001; Hardstone et al., 2012). In addition, DFA facilitates a reliable analysis of LRTC up to time scales of at least 10% of the duration of the signal (Chen et al., 2002; Gao et al., 2006), and is in this way superior to classical fluctuation analysis algorithms without a detrending step, which have low performance at large timescales (Shao et al., 2012). DFA exponents in the interval of 0.5–1.0 indicate scale-free temporal correlations (autocorrelations), whereas an exponent of 0.5 characterizes an uncorrelated signal. The main steps from the broadband signal to the quantification of LRTC using DFA have been explained in detail previously (Linkenkaer-Hansen et al., 2001; Hardstone et al., 2012). In brief, after applying the bandpass filter, the amplitude envelope was computed as the absolute of the Hilbert transform. The signal profile was computed as the cumulative sum of the demeaned amplitude envelope. After splitting the signal profile into overlapping windows (with 50% overlap) of a length *n*, the fluctuation function *F*(*t*) for each window was computed as the standard deviation of the linearly detrended signal profile. Figure 1*D* shows the average fluctuation function across all windows <*F*(*t*)> for a range of window sizes. The DFA exponent is the slope of <*F*(*t*)>, and can be interpreted as the strength of the autocorrelations in signals. For our analyses we computed DFA on the amplitude envelope of the signal filtered in the 8- to 16-Hz range, and the exponent was fit for window sizes between 5 and 30 s (Fig. 1*D*).

#### Simulations

To determine the relation between connectivity and criticality, as well as phase coupling and amplitude coupling, we set excitatory/inhibitory connectivity between 5% and 100% at 5% intervals, and five different network ensembles were created for each combination of *N* = 0 long-range connections (unconnected), and *N* = 400 long-range connections, with a conductance delay of 25 ms over the long-range connections. All simulations were run for 1000 s.

We then picked a diagonal across the critical line, where we varied excitatory connectivity between 5% and 100%, in steps of 2%. Inhibitory connectivity was computed as 105% minus excitatory connectivity. Additionally, we varied the number of long-range connections (between 0 and 1800 in steps of 100) for all networks across the diagonal, five runs per parameter combination. Separately, we also varied conductance delays (between 0 and 100 ms, in steps of 5 ms) for all networks across the diagonal, five runs per parameter combination. To analyze the effects of varying E/I balance between the two networks, we ran all combinations in one diagonal (first network) against all combinations in one diagonal (second network), five runs for each parameter combination. The same approach was used with the external input, where, in addition to conduction delays and the number of long-range connections, we also varied the time constant of the input signal,

#### Amplitude coupling

Amplitude coupling between the two networks was measured as the Spearman correlation between the amplitude envelopes of the signals filtered in the 8-16 Hz range. This value was compared with a surrogate where the causal relationship between the two signals was destroyed by swapping the two halves of one of the signals before filtering and computing the amplitude envelope.

#### Phase coupling

Phase coupling was calculated using intersite phase clustering (Lachaux et al., 2000; Cohen, 2014), i.e., the vector strength of phase differences between two signals. The advantage of this method, compared with magnitude-squared coherence, is that it does not include information about oscillation amplitude, and thereby allows us to isolate amplitude from phase coupling. First, we computed the amplitude envelope of the signals *x* and *y*, filtered in the 8-16 Hz range. Then we extracted the phase by taking the angle of the Hilbert transform

Phase coupling was then compared with a surrogate where the causal relationship between the two signals was destroyed by swapping the two halves of one of the signals before filtering and computing the phases.

#### Information transfer

Since the input connects to the first network unidirectionally, we can use direct mutual information to determine information transfer from the stimulus (Conrad and Jolivet, 2020). This approach enhances interpretability of the statistics compared with other information transfer methods. Additionally, we expect that the effect of the input on the amplitude of oscillations is a monotonous function, which allows us to further simplify the analysis, replacing mutual information with Spearman correlation. Consequently, we computed the Spearman correlation coefficient of the input signal *I* with the amplitude envelope of α oscillations (8–16 Hz) of the first network,

## Results

### Multilevel criticality is preserved after coupling networks

To model communication between two distant brain areas, e.g., frontal and parietal, we created two networks of excitatory and inhibitory neurons, and coupled them through *N* long-range connections (Fig. 1*A*). We first selected at random *N/2* excitatory positions from the first network and created synapses to excitatory neurons at the same positions in the second network, and repeated the operation for the second network, making sure that the reciprocal connections did not run through the same two neurons. We set a conduction delay of 25 ms across these long-range connections, which allows for both amplitude coupling (Fig. 2*H*) and phase coupling (Fig. 3*G*) to occur between the two networks, and is comparable with delays across long-range myelinated axons in the brain (Ringo et al., 1994). The neuronal model was based on a probabilistic firing-rate model (Dayan and Abbott, 2001), which was shown to exhibit power-law scaling of avalanche sizes (Abbott and Rohrkemper, 2007) in networks of excitatory neurons. Poil et al. (2012) adapted the model in Abbott and Rohrkemper (2007) and added inhibitory neurons (see Materials and Methods), to study the effect of E/I balance on critical dynamics. The local E/I balance was determined by the percentage of local neighbors that each excitatory and inhibitory neuron connects to. Depending on this E/I balance, the networks exhibit oscillations in the α band (8-16 Hz range; Fig. 1*B*,*C*), along with multi-level critical behavior in the form of LRTCs in the amplitude modulation of α oscillations and power-law scaling of neuronal avalanche sizes (Poil et al., 2012).

The presence of LRTC in the amplitude of oscillations (8–16 Hz) was measured by the detrended fluctuation exponent (DFA; see Materials and Methods). In brief, the signal profile was computed as the cumulative sum of the demeaned amplitude envelope. After splitting the signal profile into overlapping windows (with 50% overlap) of a length *n*, the fluctuation function *F*(*t*) for each window was calculated as the standard deviation of the linearly detrended signal profile. Figure 1*D* shows the average fluctuation function across all windows <*F*(*t*)> for a range of window sizes. The DFA exponent is the slope of <*F*(*t*)>, and can be interpreted as the strength of the autocorrelations in signals: a value of DFA close to 0.5 indicates an uncorrelated signal, and close to 1 indicates strong LRTC.

Avalanches were defined as periods of activity crossing half a median of activity (see Materials and Methods). The scale-free nature of avalanche sizes was characterized using the κ index (Shew et al., 2009; Poil et al., 2012), which compares the probability distribution of avalanche sizes with a reference power-law scale with exponent −1.5 (Zapperi et al., 1995) as follows: a value of κ = 1 corresponds to critical networks where avalanche size distributions follow the power-law closely, κ < 1 corresponds to subcritical networks which lack large-scale avalanches when compared with the reference power-law, and κ> 1 corresponds to supercritical network with an excess in large-scale avalanches (Fig. 1*E*). When the networks are not connected (*N* = 0), we see that LRTC (as measured by the DFA exponent) and scale-free avalanches occur for the same ratio of excitatory and inhibitory synapses (Fig. 1*D–G*).

After coupling networks with the same local E/I ratio, we compared the DFA and κ index of avalanches of the two networks. We found that DFA is more similar between the coupled compared with uncoupled networks [Fig. 1*H*, Pearson correlation, *p*(*r _{uncoupled}* =

*r*)

_{coupled}*I*, Pearson correlation,

*p*(

*r*=

_{uncoupled}*r*)

_{coupled}*J*) and power-law scaling of avalanche sizes (Fig. 1

*K*). Importantly, the relationship between LRTC and avalanches is preserved after coupling the networks, and depends on a specific balance between excitation and inhibition (Fig. 1

*J*,

*K*).

### Networks exhibit robust amplitude coupling only in the critical regime

To understand the role of criticality for the ability of mutually connected networks to entrain each other, we first investigated amplitude fluctuations of the coupled networks. Co-modulation of CROS was conspicuous at long time scales (Fig. 2*A*,*B*). This was confirmed by low-pass filtering the amplitude envelope of the oscillations causing cross-correlations to increase (Fig. 2*C*). This observation is in line with empirical reports on oscillations in homologous areas (Nikouline et al., 2001; Hipp et al., 2012; Siems and Siegel, 2020), suggesting that long time-scale co-modulation may emerge from mutual entrainment and does not require a slow drive such as neuromodulation. Interestingly, the capacity of the networks to couple their amplitude fluctuations breaks down both in sub and supercritical regimes (Fig. 2*C–F*), even when considering the amplitude envelope at timescales within the range of α oscillations (Fig. 2*C*). To verify that amplitude coupling emerges from a causal interaction between the two networks, we swapped the two halves of one of the signals, which entirely abolished the amplitude correlations for all networks (Fig. 2*E*,*F*). The relationship between criticality and amplitude coupling is remarkably robust to changes in the number of connections between the two networks (Fig. 2*G*) or in the conduction delay (Fig. 2*H*). This indicates that amplitude coupling is a generic feature of criticality, i.e., it is not tied to a particular parameter choice. In fact, the robustness of amplitude entrainment across parameter values suggests that local criticality might enable coupling between networks with varying connectivity strengths and situated across a wide range of physical distances. When all excitatory neurons in the two networks are coupled to each other, some supercritical networks can also establish amplitude coupling (Fig. 2*G*); however, this may not be a physiologically realistic scenario and with such a strong integration of the two networks, they may be perceived as one. To assess whether networks with different levels of criticality can also establish amplitude coupling, we allowed for the E/I ratio to vary between the two networks in the ensemble. We found that amplitude coupling emerged only when both networks were close to the critical regime (Fig. 2*I*). Thus, local criticality of brain areas appears fundamentally important for establishing amplitude coupling.

### Phase coupling emerges at criticality and extends into (slightly) supercritical regimes

The ability of mutually connected networks to phase-couple is considered important for information exchange (Fries, 2005; Palmigiano et al., 2017). We assessed the extent of phase coupling, by computing the intersite phase clustering between oscillations in the two networks (Fig. 3*A*,*B*). We found that phase coupling emerges at criticality and further increases when moving into the supercritical regime (Fig. 3*B–E*). The level of phase coupling between critical networks in the model is similar to that between fronto-central electrode in right and left hemispheres (Nikouline et al., 2001), where the intersite phase clustering also takes values close to 0.1. To verify that phase coupling emerges from a causal interaction between the two networks, we swapped the two halves of one of the signals, which entirely abolished the phase coupling between the two signals (Fig. 3*D*,*E*). Next, we varied the number of long-range connections and observed that phase coupling is stronger when there are more long-range connections (Fig. 3*F*). Importantly, the relationship between the level of criticality and phase coupling is robust to increases in the number of long-range connections (Fig. 3*F*). Varying the conductance delays, however, revealed a clear preference to phase couple in certain ranges, consistent with a coupling in-phase or out-of-phase with α oscillations (Fig. 3*G*), as previously shown in (Woodman and Canavier, 2011; Dumont and Gutkin, 2019). Last, we assessed whether networks with different levels of criticality can also establish phase coupling. Varying the E/I ratio in the two networks independently of each other, we found that phase coupling is particularly strong when a supercritical network is connected to a critical one (Fig. 3*H*).

### Information transfer is maximized close to criticality, when networks show maximal amplitude and phase coupling

To ensure that amplitude and phase coupling are not merely epiphenomena of the critical state, but actually enable networks to transfer information over long-range connections, we applied an external stimulus to the first network (Fig. 4*A*), and measured how the stimulus was encoded in the oscillatory activity of the networks. The input stimulus was Gaussian white noise convolved with an exponential with a time constant τ = 50 ms (Mainen and Sejnowski, 1995). This convolution slows down the stimulus, allowing the networks to keep up with changes in the signal (Fig. 4*B*). Since the input was connected to the first network unidirectionally, and we expected a monotonous relationship between the input and the oscillatory network activity, we used Spearman correlation to determine information transfer from the stimulus, which improves the interpretability of the statistics compared with other information transfer methods (Conrad and Jolivet, 2020). As such, we correlated the input signal with the amplitude envelope of the two networks. We found that the input is transferred to both networks (Fig. 4*C*). As expected, the signal takes longer to reach the second network, and some information is lost across networks. As a measure of information transfer, we used the maximal possible correlation (across time lags) between the input signal and the amplitude envelope of oscillations in the second network. Although the external input makes the first network more supercritical in its activity dynamics (Fig. 4*D*), information transfer is maximized when the second, coupled network, is in a regime close to criticality (Fig. 4*E*), where amplitude coupling (Fig. 4*F*) and phase coupling (Fig. 4*G*) between the two networks are maximized. To verify that information transfer emerges from a causal interaction between the input and the two networks, we swapped the two halves of the amplitude envelopes, which entirely abolished information transfer (Fig. 4*E–G*). We found that our networks can track inputs better when the input dynamics are slower (Fig. 4*H*), and when there are more long-range connections between the networks (Fig. 4*I*). Longer conduction delays over the long-range connections increase the amount of time that it takes the input to reach the second network (Fig. 4*J*); however, the peak of the correlation of the input with the second network is unaffected (Fig. 4*K*). Altogether, regardless of the parameter manipulation, information transfer is maximized when the second network is in a critical regime, where phase and amplitude coupling are strongest.

### Co-regulation of criticality and long-range communication through modulation of synaptic gain

Since neuromodulation can alter critical dynamics (Pfeffer et al., 2018), we tested whether it also affects the coupling and information transfer between networks. Neuromodulation was modeled as a multiplicative modulation of the synaptic gain of excitatory-to-excitatory connections (Fig. 5*A*) in the second network (N2), during the time intervals of 200–400 and 600–800 s (Fig. 5*B*,*C*). A modulation by a factor of 0.75 decreases overall activity (Fig. 5*B*) and brings networks toward a more subcritical regime for the duration of neuromodulation (Fig. 5*D*), while the opposite is true for a modulation by a factor of 1.25 (Fig. 5*C*,*E*). These changes in criticality are accompanied by alterations in amplitude, phase coupling and information transfer, for the duration of neuromodulation: both subcritical and supercritical networks can transfer information once they are brought toward a critical regime through appropriate changes in synaptic strength (Fig. 5*F–K*). On the other hand, network coupling and information transfer revert back to baseline once neuromodulation is turned off. These results indicate that neuromodulation can regulate the ability of networks to establish long-range communication while maintaining a fixed structure of excitatory and inhibitory connections.

## Discussion

We have investigated how the local E/I ratio of neuronal networks shapes their spatiotemporal dynamics and capacity to communicate with each other through long-range connections. We observed that the long-range coupling of neuronal oscillations, as well as information transfer through oscillatory networks is greatly facilitated when networks operate close to the critical state, where excitation and inhibition are in balance. This relationship between criticality and network coupling is robust to changes in the number of long-range connections and is present for a wide range of conduction delays. In addition to the importance of the E/I ratio of synaptic connections, we show that modulation of synaptic gain co-regulates critical dynamics and how networks couple to each other and transfer information through their oscillatory activity. Thus, local criticality and E/I balance, as determined by both local connectivity and synaptic-gain modulation, are generic mechanisms enabling long-range communication.

Our results confirmed our hypothesis that the critical state is conducive of better long-range coupling. The intuition behind this hypothesis lies in that the critical regime is characterized by transient synchrony, lying at the transition between a subcritical quiescent regime and a regime with strong, persistent oscillations (Dalla Porta and Copelli, 2019). In a subcritical regime, networks are not likely to exhibit long-range coupling: any input from an external network quickly dies out because of the strong inhibition. In the supercritical regime, on the contrary, networks are so active that any input coming from another network drowns in the relentless activity. It is only in the critical state, where localized activations have equal chance to spread or to die out, that the networks are able to phase and amplitude couple.

In our study, coupling is restricted to only two networks. However, in the brain, the constellation of areas that are functionally coupled is continuously changing (Deco et al., 2017; Tewarie et al., 2019). During rest, an intermediate interareal connection strength, which could correspond to E/I balance with critical dynamics, enables the functional connectome to spontaneously reach the highest number of possible states (Deco et al., 2017). Our study suggests that, in addition to the strength of anatomic connections, local E/I balance and the associated critical dynamics are also important for flexible interareal coupling.

The high dimensionality of the state space in the critical regime is thought to allow functional networks to form and dissolve quickly in a task-dependent manner (Deco and Kringelbach, 2016). One possible mechanism for task dependent changes in coupling is demonstrated by Palmigiano et al. (2017), where a constant input attached to one of two networks, selectively increases phase coupling and biases information transfer. In our computational model, an external input altered the criticality of the receiving network making it more supercritical. Importantly, information transfer of the external input was maximal when the second, coupled network was critical. Such a pairing of a supercritical to a critical network leads to stronger phase coupling than if both networks are critical, possibly because the supercritical network exhibits a strong oscillatory regime with the power to entrain the oscillatory phase of the susceptible critical network.

The mechanisms by which the brain regulates the criticality of oscillations in the context of long-range coupling remain understudied. Neuromodulation, e.g., by locus coeruleus, could alter the coupling gain, with effects on criticality (Pfeffer et al., 2018), and as shown in our study, also on network coupling and information transfer. Another way in which the brain can influence long-range coupling is by shifting the interareal connectivity strengths, for example under the direction of the thalamus (Nakajima and Halassa, 2017). Three additional networks have been identified, with different modulatory effects on α oscillations (Sadaghiani and Kleinschmidt, 2016): the cingulo-opercular, fronto-parietal, and dorsal attention networks. While cingulo-opercular and dorsal attention networks set the background level of α oscillations, it is the fronto-parietal network that establishes long-range phase coupling, and which is most involved in flexible cognitive control, e.g., in shaping the functional connectome during task switching (Gold et al., 2010; Marek and Dosenbach, 2018). However, coupling within the fronto-parietal network is altered in disorders like Alzheimer's disease and schizophrenia (Babiloni et al., 2004; Sheffield et al., 2015). In the light of our findings, these alterations in long-range coupling could be explained by a departure from criticality, since evidence for locally more subcritical α/β oscillations have been reported for these disorders compared with healthy controls (Montez et al., 2009; Nikulin et al., 2012). Additionally, structural connectivity is damaged in these disorders (Bozzali et al., 2002; Davis et al., 2003), which could also impact network coupling. This is consistent with our model, where a reduction in the number of long-range connections leads to impaired amplitude coupling and phase coupling, regardless of local criticality.

Phase and amplitude coupling represent distinct mechanisms (Siems and Siegel, 2020), which operate on different timescales (Daffertshofer et al., 2018): amplitude coupling is more persistent, whereas phase coupling is fast changing. However, there is also some degree of overlap between the underlying networks (Zhigalov et al., 2017; Siems and Siegel, 2020); in source-reconstructed MEG, it was found that the amplitude, phase connectomes are co-localized with the criticality connectomes (Zhigalov et al., 2017). While this provides correlational evidence for a link between the two forms of long-range coupling and criticality, our study shows that local criticality is actually a prerequisite for long-range coupling, as well as the ensuing information transfer.

Criticality in nonoscillatory systems has been associated with optimized information transmission in empirical studies, albeit only at small spatial scales in the millimeter range (multi-electrode arrays). For example, Shew and colleagues found that the critical state maximizes information transmission from an external stimulus to a local network (Shew et al., 2011). Our research complements this result, showing that criticality is also important for communication between distant networks exhibiting critical fluctuations in oscillatory activity.

When coupling two networks with the same level of local E/I ratio, we found that coupling supercritical to supercritical networks exhibit stronger phase coupling than critical-critical networks. Yang et al., reported similar results empirically, where local synchrony was found to be stronger when the local dynamics were in a supercritical regime (Yang et al., 2012). However, in the same study, the variability of synchronization was found to be reduced in the supercritical state, compared with criticality. Thus, the strong, stable phase coupling seen in the supercritical state might prevent networks from flexibly establishing coupling when needed.

We found that the criticality of local brain networks can be affected by changing the structural E/I ratio, modulating the conductance of synaptic gain of local connections or by attaching an external, excitatory input. Regardless of the way in which the population dynamics is modulated, local criticality remains a central requirement for the emergence of network coupling in our model. We predict that local criticality will also be important for network coupling of oscillatory activity in frequency ranges other than the α band considered here. Future studies should validate whether this holds empirically, considering the differences in generative mechanisms across frequency bands.

Sleep is considered important for homeostatic regulation of synaptic connections in the brain (Tononi and Cirelli, 2003). At a network level, sleep has also been related to the maintenance of E/I balance, critical dynamics, and long-range synchrony. Prolonged sleep deprivation leads to supercritical dynamics (Meisel et al., 2013), reduced synchronization of slow oscillations (Bu et al., 2017), alterations in network connectivity (Ben Simon et al., 2017), and increases in E/I balance (Chellappa et al., 2016; Bridi et al., 2020). The cognitive deficits (Alhola and Polo-Kantola, 2007) that follow sleep deprivation might be explained by altered functional integration with impaired long-range coupling. E/I balance, under the homeostatic influence of sleep, might serve as the mechanism that bridges between regulation at the synapse level and the emergent behavior.

Previous research has focused on the role of connectome topography and long-range coupling strength for brain-wide information integration. Here, we show that local E/I balance, as reflected in criticality, is a prerequisite for long-range coupling. Pharmaceutical interventions that target local E/I balance (Bruining et al., 2020) may thus help to also restore coupling and its functional consequences.

Source code required to run all simulations, as well as datasets and scripts required to generate all figures presented here, are available at https://doi.org/10.6084/m9.figshare.16645978.v1.

## Footnotes

This work was supported by Netherlands Organization for Scientific Research (NWO) Social Sciences Grant 406.15.256 (to K.L.-H. and A.-E.A.) and by the program committee Systems and Network Neuroscience of Amsterdam Neuroscience (K.L.-H.).

The authors declare no competing financial interests.

- Correspondence should be addressed to Klaus Linkenkaer-Hansen at klaus.linkenkaer{at}cncr.vu.nl