# Synchronization, stochasticity and phase waves in neuronal networks with spatially-structured connectivity

Anirudh Kulkarni, Jonas Ranft, Vincent Hakim

► To cite this version:

Anirudh Kulkarni, Jonas Ranft, Vincent Hakim. Synchronization, stochasticity and phase waves in neuronal networks with spatially-structured connectivity. Frontiers in Computational Neuroscience, 2020, 14, pp.87. 10.3389/fncom.2020.569644. hal-03028996

HAL Id: hal-03028996

https://hal.science/hal-03028996

Submitted on 4 Jan 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Synchronization, stochasticity and phase waves in neuronal networks with spatially-structured connectivity

Anirudh Kulkarni¹,²,³✉, Jonas Ranft¹,²✉✉, Vincent Hakim¹✉.

1 Laboratoire de Physique de l'Ecole Normale Supérieure, CNRS, Ecole Normale

Supérieure, PSL University, Sorbonne Université, Université Paris-Diderot, Paris, France

2 IBENS, Ecole Normale Supérieure, PSL University, CNRS, INSERM, Paris, France

3 Present address: Department of Bioengineering, Imperial College London, London, United Kingdom

✉ These authors contributed equally to this work.

✉ jonas.ranft@ens.fr or vincent.hakim@ens.fr

Version: June 4, 2020

## Abstract

Oscillations in the beta/low gamma range (10-45 Hz) are recorded in diverse neural structures. They have successfully been modeled as sparsely synchronized oscillations arising from reciprocal interactions between randomly connected excitatory (E) pyramidal cells and local interneurons (I). The synchronization of spatially distant oscillatory spiking E-I modules has been well studied in the rate model framework but less so for modules of spiking neurons. Here, we first show that previously proposed modifications of rate models provide a quantitative description of spiking E-I modules of Exponential Integrate-and-Fire (EIF) neurons. This allows us to analyze the dynamical regimes of sparsely synchronized oscillatory E-I modules connected by long-range excitatory interactions, for two modules, as well as for a chain of such modules. For modules with a large number of neurons (> 10⁵), we obtain results similar to previous obtained ones based on the classic deterministic Wilson-Cowan rate model, with the added bonus that the results quantitatively describe simulations of spiking EIF neurons. However, for modules with a moderate (~ 10⁴) number of neurons, stochastic variations in the spike emission of neurons are important and need to be taken into account. On the one hand, they modify the oscillations in a way that tends to promote synchronization between different modules. On the other hand, independent fluctuations on different modules tend to disrupt synchronization. The correlations between distant oscillatory modules can be described by stochastic equations for the oscillator phases that have been intensely studied in other contexts. On shorter distances, we develop a description that also takes into account amplitude modes and that quantitatively accounts for our simulation data. Stochastic dephasing of neighboring modules produces transient phase gradients and the transient appearance of phase waves. We propose that these stochastically-induced phase waves provide an explanatory framework for the observations of traveling waves in the cortex during beta oscillations.

## Introduction

Rhythms and collective oscillations at different frequencies are ubiquitous in neural structures [1]. Numerous works have been devoted to understanding their origins and
characteristics [2] which depend both on the neural area and on the activity of the animal. Gamma band oscillations (30-100 Hz) are for instance recorded in the visual cortex as well as several other structures and have been hypothesized to support various functional roles [3]. Beta oscillations (10-30 Hz) are prominent in the motor cortex during movement planning before movement initiation [4] and have traditionally been assigned a role in movement control while more general roles have also been proposed [5].

Several experimental results point out the need to model and analyze the spatial organization of oscillatory activity [6]. An early study [7], using widefield imaging and voltage-sensitive dyes, reported that stimulus-induced oscillatory activity around 10 Hz and 20 Hz was organized in plane waves and spiral waves in the turtle cortex. This spiral-like organization was also reported for pharmacologically induced 10 Hz oscillations in the rat visual cortex [8]. The underlying mechanisms and specific cells involved in the synchronization of two distant regions have more recently started to be investigated using optogenetic manipulations in mice [9].

A motivating illustrative example for the present theoretical study is the observation that beta oscillatory activity during movement preparation exhibits transient episodes of propagating waves in the motor cortex of monkeys [10–13] and humans [14]. These waves appear to propagate along particular anatomical directions with a typical wavelength of 1 cm and a velocity of about 20 cm·s⁻¹ [10]. On the structural side, the characterization of the long-range intracortical connectivity in motor cortex have given rise to several quantification endeavors in both monkeys [15, 16] and cats [17] with the strength of excitatory responses to microstimulation decaying on a 2 mm length scale. Can one connect these different observations in some ways and how? What is the proper mechanistic interpretation of these observed waves?

The study of waves in oscillatory media has been largely based on the analysis of oscillator synchronisation which itself has a long history [18]. Classic mathematical methods for studying synchronization of weakly coupled oscillators have been extended to study oscillations and travelling waves in spatially extended media in physics and chemistry [19]. In simple descriptions of neural network dynamics based on rate models introduced by Wilson and Cowan [20], also called neural-mass models, the dynamics of a whole set of neurons is reduced to a small set of differential equations approximately describing the temporal evolution of the firing rate of a “typical ”’neuron in the set. This allows one to directly apply the techniques developed to analyze synchronization of oscillators to study the synchronization properties of sets of oscillating neurons in the rate-model framework. This approach has been followed in a number of works to study the synchronization of spatially-coupled neural networks in the oscillatory regime [21–25]. Rate models with spatially-structured connectivity have also been extensively studied to analyze pattern formation in a neural network context since Amari’s work [26] (see [27–29] for reviews). A limitation is that rate models are generally difficult to quantitatively relate to models of single neurons.

Networks of model spiking neurons provide a more detailed description of neuronal network dynamics than rate models. Studies of networks with random unstructured connectivities have shown the existence of a “sparsely synchronized” oscillatory (SSO) regime [30] in which a collective oscillation exists at the whole population level while spike emission by single neurons is quite stochastic with no significant periodic component. Experimental recordings suggest that neural rhythms in the beta, gamma and higher frequency ranges operate in this regime [2]. Specifically, rhythms in the beta/low gamma range are thought to arise from sparse synchronization between excitatory (E) and inhibitory (I) neuronal populations with reciprocal interactions [31, 32], the so-called Pyramidal-Interneuron Gamma (PING) mechanism [33].

A few studies have considered the impact of structured connectivity on sparsely synchronized oscillations in the spiking-network framework. The influence of delays
in synaptic transmission was studied in [34] for the synchronization of two neuronal E-I modules oscillating in the high-gamma frequency range, using a classic rate-model formalism. A variety of dynamical regimes was found and qualitative agreement with the bifurcations observed in simulations of a spiking network was reported when the neurons in each E-I module oscillated in a well-synchronized manner. Visually induced gamma oscillations and their dependence on visual stimulus contrast [35] have been investigated by simulations of a spiking network with a two-layer ring-model architecture. Recently, simulations of 2 or 3 coupled E-I modules have been performed [36] to assess how information flow between different modules is correlated with bursts of transient synchrony at gamma frequency.

Spiking networks model more closely biological reality but are more difficult to analyze theoretically. The mathematical analysis of their oscillatory regime is essentially confined to the neighborhood of the oscillatory threshold in parameter space (parameters being the strengths of synaptic connections). General linear stability analyses of E-I spiking networks with spatially structured connectivity have been performed in the absence of transmission delays [37] or taking one into account [38]. Going beyond linear stability is feasible but already somewhat heavy for a single module with unstructured connectivity [30, 39]. An exception is the soluble case of deterministic quadratic integrate-and-fire neurons with a wide distribution of frequencies [40] which has been used to study oscillations of an E-I module [41] as well as synchronization between two weakly-coupled modules [42]. Aside from oscillations, the conditions necessary for the existence of a balanced state [43] in a spatially structured network have been studied [44] and it was found in particular that the spatial spread of excitation should be broader than that of inhibitory inputs. Firing correlations have also been studied in such a network state [45]. However, the spatial organization of sparsely synchronized oscillatory activity in spiking networks with spatially structured connectivity remains to be generally studied.

In the present study, our aim is to analyze the synchronization between different local E-I modules, induced by distance-dependent long-range excitation. We focus on the case where each module oscillates in a sparsely synchronized way [30] to model the *in vivo* situation. We combine spiking network simulations in the SSO regime with a mathematical analysis of an “improved” rate model to develop a quantitative picture of the dynamics in such a system. Relating rate models to spiking networks is a classical endeavor. This can be achieved when the collective dynamics is stationary or slowly varying [46], for instance in the presence of slow synapses [47], by an appropriate choice of the f-I curve in the rate model. More recently, different works have shown that mild modifications allow rate models to overcome this limitation of a slow-varying collective dynamics by introducing a timescale that depends of the firing rate. These “adaptive” rate models successfully produced a quantitative description of an uncoupled neuronal population submitted to time-varying inputs [48], as well as of spike synchronization [49] or oscillations driven by spike-frequency-adaptation (SFA) [50] in a recurrently coupled excitatory neural population. Building on this progress, we first develop a rate model with an adaptive timescale [48, 50] and show that it accurately describes population-level oscillations of an E-I spiking neuron module in the SSO regime. This allows us to make use of the large body of work that has been developed to study synchronization in deterministic rate-model equations [23–25]. We show that in spite of the introduction of the adaptive timescale, the model that we use behaves very similarly [21–24] to classic rate models [20]. For two modules, we find that the synchronization between the two-module oscillations depends on the specific pattern of long-range excitation connectivity. Complex dynamical regimes are produced when long-range excitation is weak and targets only excitatory neurons [21, 22], whereas synchronization of the two oscillations of the two modules is otherwise observed. For a chain of oscillatory E-I
modules with long-range excitatory coupling decreasing with distance, we similarly find that the connectivity properties of long-range excitation play an important role. In particular, when long-range excitation only targets excitatory neurons, distant modules oscillate with different phases, namely oscillation phase gradients and phase waves spontaneously appear. The bonus as compared to the classical rate-model results, is that, as we show, the obtained results quantitatively agree with simulations of spiking networks of large size.

Simulations of spiking modules of biologically relevant size, of about $10^4 - 10^5$ neurons each, display significant stochastic variations in module activities that require to go beyond the use of deterministic rate models. Building on previous work [39], we find that finite-size fluctuations can be quantitatively accounted in the adaptive-timescale rate-model framework by adding to it a stochastic component. Even when different modules would fully synchronize in a deterministic context, stochasticity in the module activities produces differences in the phases of oscillation of neighboring modules. [18, 19, 25]. Classical [19] and more recent results [25, 51, 52] show how to derive stochastic equations for the oscillatory phases of the different modules that describe stochastic dephasing. In our case, we find that this usual weak-noise reduction to a phase dynamics is accurate only for modules with a very large number of neurons. We obtain a quantitative description of stochastic effects for modules of moderate size, by computing the noise influence on both the phase and amplitude of the module oscillation at the linear level. Stochastic dephasing creates transient bouts of traveling waves in a chain of modules that would be fully synchronized in a deterministic context. We end by obtaining in a chain of E-I modules with long-range excitatory coupling, the probability of dephasing between close E-I modules and the spectrum of phase velocities for the corresponding bouts of traveling waves. We propose that this phenomenon of phase waves induced by stochasticity is at the root of the observation of waves during oscillatory episodes in cortex.

# Results

## Oscillations in an E-I module: rate-model description vs spiking network simulations

We first consider oscillations of neural activity in a local module comprising excitatory (E) and inhibitory (I) neuronal populations with a spatially unstructured connectivity. Oscillations are observed both in rate-model descriptions [20] and in simulations of spiking networks [31–33]. We choose the neurons in our spiking network simulations to be of the Exponential-Integrate-and-Fire type (EIF) [53] (see Eq. (25) for their mathematical definition), which have been shown to describe well the dynamics of cortical neurons [54].

We wish to use a rate model description that quantitatively describes oscillations of the spiking E-I module in the SSO regime [30]. Mild modifications of classic rate models [20] can give quantitatively accurate accounts of simulations of spiking neurons when spike emission has a strong stochastic component, as it is the case in the SSO regime. This was previously demonstrated both for an assembly of independent spiking neurons receiving identical inputs [48] and for spiking neurons coupled by recurrent excitation [50]. Building on these advances, we use the following adaptive timescale rate-model formulation [48, 50] to represent the activity of a group of EIF neurons:

$$ \tau(I) \frac{dI}{dt} = -I + I_0 + s(t), \qquad (1) $$

$$ r(t) = \Phi_{\sigma}[I(t)], \qquad (2) $$
where $\Phi_{\sigma}(I)$ is the f-I curve of the EIF model (25) for a noisy input current $I_{ext} + I_{syn}$ with mean $I$ and total noise strength $\sigma$. It is shown in Fig. S1A for the parameters used in the present study. This specific choice of f-I curve allows the rate-model to to agree with the rate of a spiking network of spiking EIF neuron in a stationary regime. The term $s(t)$ represents the time-dependent inputs to the neurons of the population and includes both the external and the network-averaged instantaneous internal inputs. The response time $\tau(I)$ is said to be “adaptive” [48] because it depends on the current $I$. It allows to overcome the limitation of a slowly varying dynamics [46, 47] when the population rate remains close to the stationary f-I curve. The chosen function $\tau(I)$ is displayed in Fig. S1B. It is obtained by a fitting procedure to best reproduce the dynamics of a population of uncoupled EIF spiking neurons (see Methods: Fit of the adaptive timescale in the FAT rate model for a precise description of our fitting procedure and Fig. S1C for examples of the timescale fit). In agreement with previous works [48, 50], Fig. S1D-F shows that the rate model with such a fitted adaptive timescale (FAT model) reproduces quite accurately oscillations of the population activity of a set of identical uncoupled neurons driven by a sinusoidal input current.

The connectivity of our elementary module is depicted in Fig. 1A (inset). It is composed of an excitatory (E) and an inhibitory (I) population with recurrent connections. The FAT rate-model description is obtained by describing the dynamics of each of these populations by Eq. (2),

$$ \tau_E(I_E) \frac{dI_E}{dt} = -I_E + I_E^{ext} + w_{EE} r_E - w_{EI} r_I, \quad (3) $$

$$ \tau_I(I_I) \frac{dI_I}{dt} = -I_I + I_I^{ext} + w_{IE} r_E, \quad (4) $$

with $r_E = \Phi_E(I_E)$ and $r_I = \Phi_I(I_I)$. (Here and in the following we have dropped the explicit noise strength index on $\Phi(I)$ for notational simplicity). In all numerical computations, we have chosen $\Phi_E(I) = \Phi_I(I) = \Phi_{\sigma}(I)$, $\tau_E(I) = \tau_I(I) = \tau^{(FAT)}(I)$, with $\Phi_{\sigma}(I)$ and $\tau^{(FAT)}(I)$ displayed in Fig S1. Note that recurrent inhibition has not been included in Eqs. (3,4) both to simplify the analysis and to focus on oscillations in the beta/low gamma range that are mediated by E-I reciprocal interactions. Recurrent inhibitory connections allow for high frequency oscillations in suitable parameters regimes [32, 39] and can also be important to prevent a too high firing rate of interneurons, e.g. in the balanced state [43], a role that is assigned here to the external input $I_I^{ext}$.

The stability diagram for the steady non-oscillatory regime of the FAT rate model is easily obtained [20, 24] (see Methods: Oscillatory instability for the E-I module), see Fig. 1A. As for the classic Wilson-Cowan model [20], it displays three regimes. Steady neural activity is stable when recurrent excitation, measured by the total synaptic strength $w_{EE}$, is weak enough. When recurrent excitation grows, two possibilities arise. They depend on the strength of inhibitory feedback on the excitatory population, measured by the product $w_{IE}w_{EI}$, where $w_{IE}$ is the total excitatory synaptic strength on inhibitory neurons and $w_{EI}$ the total inhibitory synaptic strength on excitatory neurons. When inhibitory feedback is weak, the steady state is subject to a non-oscillatory instability: recurrent excitation leads to steady firing at a very high rate, basically limited by the neuron refractory period in our simple model (other mechanisms, not considered here, such as SFA or pair-pulse synaptic depression can moderate this regime). When inhibitory feedback is sufficiently strong, the steady state is destabilized by an oscillatory instability. This instability can lead to finite amplitude oscillations but also to steady high frequency discharge when a steady high-rate fixed point exists (see [55] for a detailed analysis). Oscillations with high discharge and synchronous spiking are also possible outcomes. Here, we limit ourselves to considering oscillations of moderate amplitude that remain sparsely synchronized [30], a dynamical regime which appears most appropriate to describe beta/low gamma oscillations recorded in vivo. The
<figure><img src="image_2.png" /><figcaption>**Fig 1.** Dynamical regimes of an E-I module. (A) Stability and instability of the module steady discharge as a function of synaptic coupling. Stable regions with complex eigenvalues (green) and real eigenvalues (light blue) and unstable regions with complex eigenvalues (purple) and real eigenvalues (orange) are shown for the stationary state $r_E^{(s)} = 5$ Hz, $r_I^{(s)} = 10$ Hz. Inset: sketch of the excitatory and inhibitory neuronal populations (the E-I module) and their synaptic interactions. The parameters for the three cases shown are: A - $w_{EE} = 1.6$ mVs, $w_{IE}w_{EI} = 0.64$ mV$^2$s$^2$; B - $w_{EE} = 1.6$ mVs, $w_{IE}w_{EI} = 1.28$ mV$^2$s$^2$; C - $w_{EE} = 1.76$ mVs, $w_{IE}w_{EI} = 1.28$ mV$^2$s$^2$. We use the parameters of case A (solid circle) throughout the remainder of this article; for cases B and C, see Fig. S2. (B,C) Oscillations of the discharge rates $r_E, r_I$ of the E-I populations in the rate-model description, (B) as $r_I$ vs. $r_E$ or (C) as a function of time (blue, excitatory population; red, inhibitory population) for parameters corresponding to case A in panel (A) with $w_{IE} = 2$ mVs, $w_{EI} = 0.32$ mVs (the corresponding effective constants are $\alpha = 2.33, \beta = 2.15$). Oscillations in a spiking network with a large number of neurons ($N = 10^6$, resampled with time bin $\Delta t = 0.1$ ms) are also shown in (B) (blue line) and (C) (black and gray lines). (D) Time traces of the excitatory activity in a spiking network with smaller number of neurons ($N = 10^4$) and the corresponding stochastic rate model Eqs. (5,6) (sampled with time bin $\Delta t = 0.1$ ms). (E) Autocorrelation of the excitatory activity for the spiking EIF module ($N = 10^4$, solid blue) and the adaptive rate model (solid orange). The autocorrelation for the rate model with the adaptive timescale of ref. [48] is also shown (solid green). The fit of the analytical expression Eq. (56) for the autocorrelation to the adaptive rate model (shown in dashed black) allows to obtain the module's autocorrelation decay time $\tau_D$. (F) The decorrelation time $\tau_D$ shown for spiking modules (solid blue dots), for FAT modules (solid orange line) and as predicted by Eq. (10) (dashed black line).</figcaption></figure>
parameters corresponding to one such point are indicated in Fig. 1A (solid black circle) and used as reference for the figures of the present paper. The corresponding limit-cycle oscillations are shown in Fig. 1B,C. Two other representative sets of parameter values are also indicated in Fig. 1A (solid losange and square) and the corresponding dynamical traces are shown in Fig. S2.

The rate-model description is compared to simulations of spiking networks (see Methods: Simulations of spiking networks) in Fig 1. The rate-model deterministic activity trace accounts well for the network activity when the number N of spiking neurons is very large ($N \sim 10^6$) and stochastic effects at the level of the population are negligible, as shown in Fig 1B,C.

Our E-I module with unstructured connectivity is, however, intended to represent local interactions at a scale comparable to that of a cortical column [56,57], which is estimated to comprise a smaller number of neurons ($N \sim 10^4 - 10^5$). In this case, simulations show that the network population activity has a significant stochastic component, as seen in Fig. 1D. Auto-correlograms of the E population activity display decreasing oscillatory tails reflecting the corresponding dephasing of oscillations (Fig. 1E).

We take into account these stochastic effects in the rate-model description by remembering that, in a sparsely synchronized spiking network [30], the global network activity retains a stochastic component due to the finite numbers $N_E, N_I$ of excitatory and inhibitory neurons in the network, even for constant external inputs. We follow ref. [39] and assume that, in the SSO regime, the probability of an excitatory (resp. inhibitory) spike between the times $t$ and $t + \Delta t$ is given by $\Phi_E[I_E(t)]\Delta t$ (resp. $\Phi_I[I_I(t)]\Delta t$), and spikes being independently drawn for each neuron in the network. This leads to replacing the deterministic relations $r_E(t) = \Phi_E[I_E(t)]$ and $r_I(t) = \Phi_I[I_I(t)]$ by stochastic versions coming from Poissonian sampling, an approximation that has previously been shown to be quite accurate when neurons are sparsely synchronized [39],

$$r_E(t) = \Phi_E[I_E(t)] + \sqrt{\Phi_E[I_E(t)]/N_E} \xi_E(t), \quad (5)$$

$$r_I(t) = \Phi_E[I_I(t)] + \sqrt{\Phi_I[I_I(t)]/N_I} \xi_I(t), \quad (6)$$

where $\xi_E(t), \xi_I(t)$ are independent white noises ($\langle \xi_E(t)\xi_E(t') \rangle = \langle \xi_I(t)\xi_I(t') \rangle = \delta(t - t')$, $\langle \xi_E(t)\xi_I(t') \rangle = 0$) and Ito's prescription is used [58]. Eq. (5,6) transform the deterministic rate equations into stochastic ones with a noise amplitude that is inversely proportional to the size of the population

Accounting in this way for finite-size fluctuations allows the FAT rate-model description to reproduce quite well the sparsely synchronized oscillations in a moderately-sized spiking network (Fig. 1D). The similarity between the spiking network and the stochastic FAT model can be quantitatively assessed by computing the autocorrelations of the E-I module excitatory population activity, $C_{EE}(t - t') = \langle r_E(t)r_E(t') \rangle - \langle r_E \rangle^2$. As shown in Fig. 1E, the stochastic FAT model accurately describes the network autocorrelation. It should be noted that the added stochastic terms in the rate equation description have no free parameters, they are entirely determined by the assumption (admittedly approximate) that the instantaneous network spike rate is the product of an underlying Poissonian process.

Stochastic dephasing of oscillations leads to an exponential decrease of the autocorrelation amplitude with a characteristic time $\tau_D$,

$$|C_{EE}(t)| \sim \exp(-t/\tau_D). \quad (7)$$

The time $\tau_D$ increases with the number of neurons in the E-I module as shown in Fig. 1F.

Well-known results are available to analytically describe dephasing due to weak noise [19, 25, 51, 52], namely in the present case when the E-I module comprises large numbers of neurons and finite-size fluctuations around the deterministic limit cycle
are small (see Methods: Diffusion of the oscillation phase of an E-I module with a finite number of neurons). For very large numbers of neurons, the activities of the two populations ($r_E(t + \phi)$, $r_I(t + \phi)$) are oscillatory with $r_E(t)$ and $r_I(t)$ periodic functions of time, and the phase $\phi$ being arbitrary but constant. Note that in the present work, we define the phase $\phi$ to be a variable with the same dimension as time and with a period $T$ equal to the oscillation period. Another usual convention is to consider the phase as a dimensionless variable with a period of $2\pi$. (Therefore $r_E$ is here a periodic function of the phase $\phi$ with period $T$, instead of being a periodic function of a phase with period $2\pi$.) The phases in these two different conventions are simply related by a multiplicative factor of $2\pi/T$.

Stochastic fluctuations in E-I modules with large number of neurons dominantly produce a random drift of the phase $\phi$ with a diffusive behavior in time,

$$ \langle [\phi(t) - \phi(0)]^2 \rangle = D_N t. \qquad (8) $$

Similarly to the amplitude of the finite-size fluctuations of activity, the diffusion constant $D_N$ vanishes as the numbers $N_E$ of excitatory neurons and $N_I$ of inhibitory neurons in the module grow,

$$ D_N = \frac{D_E}{N_E} + \frac{D_I}{N_I}. \qquad (9) $$

The constants $D_E$ and $D_I$ can be expressed and computed in terms of integrals along the deterministic limit cycle [19] (Eqs. (53,54) in Methods). This provides an explicit expression of the decorrelation time as a function of $D_N$ and the period of the limit cycle,

$$ \tau_D = \frac{1}{2\pi^2} \frac{T^2}{D_N}. \qquad (10) $$

Eq. (10) predicts that the dephasing time increases linearly with the size $N$ of the E-I module. Expression (10) is displayed in Fig. 1F. It agrees well with the simulation results for large $N$. Quantitative agreement deteriorates as $N$ diminishes and fluctuations become stronger. Fluctuations then strongly affect the shape of limit cycle itself and the description by a pure dephasing becomes less accurate.

## Dynamical regimes of two oscillatory E-I modules coupled by long-range excitation

We start by considering the simplest case of coupling between distant modules, namely two identical E-I modules coupled by long-range excitation. In spite of the addition of the adaptive timescale in the rate model we use, the results we obtain are very similar to results [21, 22] obtained for the classic Wilson-Cowan model [20]. It matters for synchronization whether long-range excitation targets only excitatory neurons or both excitatory and inhibitory neurons. We thus distinguish the two cases.

### Long-range excitation targeting excitatory neurons

We analyze the synchronization properties between two E-I modules with oscillatory activity when long-range excitation only targets excitatory neurons, the “$E \to E$” connectivity case depicted in Fig. 2A.

The dynamics of two coupled E-I modules are described in the rate-model framework (3,4) by

$$ \tau_E(I_{E,1}) \frac{dI_{E,1}}{dt} = -I_{E,1} + I_E^{ext} + w_{EE} [(1 - f_{lr})r_{E,1} + f_{lr}r_{E,2}] - w_{EI} r_{I,1} \quad (11) $$

$$ \tau_I(I_{I,1}) \frac{dI_{I,1}}{dt} = -I_{I,1} + I_I^{ext} + w_{IE} r_{E,1} \qquad (12) $$
<figure><img src="image_3.png" /><figcaption>Fig 2. (see next page)</figcaption></figure>

with the same two equations with permuted indices 1 and 2 describing the dynamics of module 2.

We first consider the dynamics of this two-module network for two coupled deterministic FAT rate models, when the firing rates, $r_{E,n}$ and $r_{I,n}$ for $n = 1, 2$, are given in terms of the respective currents by the f-I curve (Eq. (2)). Mathematical analysis and simulations show the existence of a number of different dynamical regimes. In order to obtain a precise view of the various cases as a function of the different synaptic couplings, we choose the parameters of the E-I modules so that they remain at a fixed oscillatory location in the parameter diagram of Fig. 1A. This fixes in each module the total strength of recurrent excitatory connections $w_{EE}$ and the total inhibitory feedback $w_{IE}w_{EI}$, and leaves as variable parameters the strength of excitation on inhibitory interneurons $w_{IE}$ (or equivalently the strength of inhibition on excitatory neurons $w_{EI}$)
**Fig 2 (previous page). Synchronization and dynamical regimes for two E-I modules coupled by long-range excitation.** (A) Sketch of the coupled E-I modules and their synaptic interactions. Top: In the $E \to E$ connectivity case, long-range excitation only targets excitatory neurons in other modules. Bottom: In the $E \to E, I$ connectivity case, long-range excitation targets excitatory and inhibitory neurons in distant modules in equal proportion. (B) Different dynamical behaviors for the coupled modules as a function of the coupling strength $f_{lr}$ and the synaptic excitatory strength on inhibitory neurons $w_{IE}$. The labeled points (solid black circles) mark the parameters of the dynamical regimes shown in (B1-6), the cross (light grey) corresponds to the parameter used in Fig. 3A. (B1) Fully synchronized state for sufficiently strong long-range excitation; $f_{lr} = 0.03$. (B2) Finite phase difference: the activity of one module is greater than the other, the amplitudes of their oscillatory activities are constant in time but the phases of their oscillations are different; $f_{lr} = 0.025$. (B3) Modulated dominance: one of the two modules is more active than the other but the amplitudes of the oscillatory activities of both modules vary themselves periodically in time; $f_{lr} = 0.02$. (B4) Alternating dominance: the activities of the two E-I modules successively dominate; $f_{lr} = 0.015$. (B5) Antiphase regime at very low coupling and strong local excitation on inhibitory neurons; $f_{lr} = 0.001$. (B6) Finite phase difference regime at very low coupling and weak local excitation on inhibitory neurons; $f_{lr} = 0.001$ and $w_{IE} = 1.44$ mVs. (C) The synchronization function governing the evolution of the relative of two weakly coupled modules in the $E \to E$ connectivity case. Modules with weak local excitation of inhibitory neurons $w_{IE}$ (blue solid line), strong local excitation of inhibitory neurons (orange dashed-dotted line) and $w_{IE}$ close to the threshold strength of local excitation separating the two regimes (red dashed line). (D) The synchronization function governing the evolution of the relative phase of two weakly coupled modules in the $E \to E, I$ connectivity case. The fully synchronized state is the only stable state (zero-crossing with negative slope) for the three shown $w_{IE}$ values. (E) The difference of excitatory activity between the two modules decreases exponentially ($w_{IE} = 2$ mVs, $f_{lr} = 0.001$). (F) The measured exponential restabilisation in two-modules FAT simulation with $E \to E, I$ connectivity (crosses, $f_{lr} = 0.001$) matches well the prediction of Eq. (15) and Eq. (66) (solid black line).

and the fraction $f_{lr}$ of excitatory connections on excitatory neurons that corresponds to long-range excitation. The found dynamical regimes of the two E-I networks are displayed in Fig. 2B as a function of $w_{IE}$ and $f_{lr}$ for one of our reference points (case A marked by the black solid circle in Fig. 1A; the analogous diagrams for the two other points are shown in Fig. S3). The individual panels Fig. 2B1-6 show representative traces of the activities of the two excitatory populations in the different regimes. We describe them in turn.

When the E-I modules are not coupled, each of them oscillates with a period $T$ and an arbitrary phase with respect to the other module. For very weak long-range excitation, the effect on each E-I module of the excitatory inputs coming from the other one simply leads to a slow change of the phase of its oscillation [19, 21, 23, 24]. The synchronization dynamics itself can be characterized by the evolution of the relative phase $\Delta\phi$ between the oscillations of the two modules. For $\Delta\phi = 0$, the two modules oscillate in phase whereas for $\Delta\phi = T/2$ they oscillate in antiphase. The evolution of $\Delta\phi$ is found to be governed by the following dynamics [19, 21, 23, 24] (see Methods: Synchronization function for two weakly coupled E-I modules):

$$ \frac{d\Delta\phi}{dt} = f_{lr} S_E(\Delta\phi). \qquad (13) $$

The function $S_E(\Delta\phi)$ is shown in Fig. 2C for different synaptic parameter values $w_{IE}$. For two identical modules, symmetry implies that the phase differences $\Delta\phi = 0$ and $\Delta\phi = T/2$ are always zeros of $S_E$. Therefore, in-phase as well as antiphase oscillations are always possible oscillating states. In the present case, the positive slope of the zero-crossing at $\Delta\phi = 0$, $S'_E(\Delta\phi) > 0$ shows that in-phase oscillations are always unstable. When local excitation targets sufficiently strongly interneurons, i.e. for $w_{IE} > w_{IE}^*$,
$S'_E(T/2) < 0$ and antiphase oscillations are the only possible stable phase difference. An example of such oscillations is displayed in Fig. 2B5. On the contrary for $w_{IE} < w_{IE}^*$, both in-phase and antiphase are unstable. The function $S_E(\Delta\phi)$ displays a new zero at a non-trivial $\Delta\phi$ with a negative slope (Fig. 2C) and the two modules stably oscillate with this non-trivial phase difference (Fig. 2B6). Such a synchronization regime was also obtained [21] for the classic Wilson-Cowan rate model [20].

For a larger strength of long-range excitation, the interaction between the two modules does not reduce to changing their oscillation phases. Nonetheless, the stability of the in-phase oscillations can be computed (see Methods: Stability analysis of full synchronization for two coupled E-I modules). The analysis shows that fully synchronous oscillations (Fig. 2B1) are stable when long-range excitation is sufficiently strong. This actually happens for a relatively weak long-range excitation strength, at about $f_{lr} \sim 2-7\%$, depending on the strength $w_{IE}$ of the local excitation on the inhibitory population in each module (Fig. 2B).

More complex dynamical regimes hold for intermediate strengths of the long-range excitation, below that required for full synchrony and above the excitation strength leading to finite-phase difference phase at very low $f_{lr}$. A first decrease of $f_{lr}$ below that necessary for full synchrony leads to a phase where the two-modules oscillate at the same frequency and different phases. In contrast to the finite-phase difference at very low $f_{lr}$, oscillation amplitude in one module strongly dominates the other one, as illustrated in Fig. 2B2. This finite phase difference regime itself looses stability with a further decrease of $f_{lr}$ and is replaced by a “modulated dominance” regime: beta/low gamma oscillations in one module continue to be of larger amplitude in one module than in the other one but these two unequal amplitudes are themselves modulated at a very low frequency in the 1 Hz range, as shown in Fig. 2B3. This regime stems from a Hopf bifurcation of the finite-difference regime, with the unusual low frequency coming from the small value of $f_{lr}$ (see Methods: Stability analysis of full synchronization for two coupled E-I modules). A further decrease of $f_{lr}$ transforms the modulated-dominance regime into an “alternating dominance” regime: oscillations in one module dominate the oscillations in the other one, but the dominant module switches periodically, again with a slow frequency in the 1 Hz range, as illustrated in Fig. 2B4. Similarly to modulated dominance, alternating dominance stems from a Hopf bifurcation of the antiphase regime at very low $f_{lr}$, when $f_{lr}$ is increased. A sequence of bifurcations and dynamical states similar to the one reported in Fig. 2 has been described long ago [22] for the Wilson-Cowan rate model upon increase of the excitatory coupling between the two modules. The regimes that we call here “alternating dominance” and “modulated dominance” are referred to as “antisymmetric torus” and “nonsymmetric torus” in ref. [22] and a short discussion of the transition between them is provided.

## Long-range excitation targeting excitatory and inhibitory neurons

We now consider the case in which long-range excitation targets both excitatory and inhibitory neurons with strengths proportional to that of local excitation. In this case, Eq. (12) is replaced by

$$ \tau_I(I_{I,1}) \frac{dI_{I,1}}{dt} = -I_{I,1} + I_I^{ext} + w_{IE}[(1-f_r)r_{E,1} + f_r r_{E,2}] \quad (14) $$

with the same equation with permuted index 1 and 2 to describe the dynamics of the inhibitory population of module 2.

In this “$E \to E, I$” connectivity case, sketched in Fig. 2A, long-range excitation is found to always have a synchronizing effect.

For very weak coupling, the synchronization dynamics can be reduced as above to the evolution of the relative phase of the two modules. It is described by Eq. (13) with the
<table><caption>**Table 1.** Computed constants. Values of the constants defined in the text for the E-I module limit cycle with the reference parameters (solid circle in Fig. 1A) $w_{EE} = 1.6$ mVs, $w_{EI} = 0.32$ mVs, $w_{IE} = 2.0$ mVs.</caption><thead><tr><th>Constant</th><th>Value</th><th>Definition</th><th></th></tr></thead><tbody><tr><td><i>T</i></td><td>63.7 ms</td><td>limit cycle period</td><td></td></tr><tr><td><i>D</i><sub>E</sub></td><td>1.2 &middot; 10<sup>4</sup> ms</td><td>Excitatory noise amplitude</td><td>(Eq. (53))</td></tr><tr><td><i>D</i><sub>I</sub></td><td>2.0 &middot; 10<sup>3</sup> ms</td><td>Inhibitory noise amplitude</td><td>(Eq. (54))</td></tr><tr><td><i>D</i><sub>N</sub></td><td>2.5 ms (for <i>N</i> = 10<sup>4</sup>)</td><td>Finite-size noise strength</td><td>(Eq. (9))</td></tr><tr><td><i>D</i><sub>&phi;</sub><sup>E</sup></td><td>-0.31 ms<sup>-1</sup></td><td>Phase diffusion for <i>E</i> &rarr; <i>E</i> coupling</td><td>(Eq. (63))</td></tr><tr><td><i>D</i><sub>&phi;</sub><sup>EI</sup></td><td>9.4 &middot; 10<sup>-2</sup> ms<sup>-1</sup></td><td>Phase diffusion for <i>E</i> &rarr; <i>E</i>, <i>I</i> coupling</td><td>(Eq. (66))</td></tr></tbody></table>

function $S_{EI}(\Delta\phi)$ plotted in Fig. 2D. For a small phase difference, $S_{EI}(\Delta\phi)$ behaves as

$$S_{EI}(\Delta\phi) \approx -2\Delta\phi D_{\phi}^{EI}, \qquad (15)$$

with the constant $D_{\phi}^{EI}$ depending on the characteristics of the oscillation cycle (see Eq. (66) in Methods). The constant $D_{\phi}^{EI}$ is found to be positive (see Table 1 for its value at the reference point $w_{IE} = 2$ mVs). Eq. (15) shows that a zero phase lag between the two modules is a dynamically stable configuration, in contrast to our previous results for $E \to E$ connectivity (compare with Fig. 2C). Eq. (15) predicts that an initial activity difference between two oscillating E-I modules should vanish exponentially. As shown in Fig. 2E, this is indeed observed in simulations with a rate of decrease that agrees well with the one predicted by Eq. (15) (see Fig. 2F).

We have checked that this synchronizing effect of long-range excitation persists for stronger long-range excitation strengths. We have found the two-module fully synchronized regime to be linearly stable for all couplings for which we computed its stability matrix (see Methods: Stability analysis of full synchronization for two coupled E-I modules). Moreover, numerical simulations with different initial conditions did not show any other stable pattern. Therefore, we conclude that in the context of the present model with instantaneous synapses and no propagation delays, long-range excitation with $E \to E, I$ connectivity between two E-I modules always tends to fully synchronize their activities.

## Comparison with spiking networks and influence of stochastic activity fluctuations

For a single E-I module, we already found a good match between the oscillatory dynamics in the rate-model framework and simulations of a spiking E-I network, as described in the first section of the results. Finite-size fluctuations were found to play an important role for modules of moderate sizes and were found to be well accounted for by our stochastic rate model. We now pursue this comparison in the two-module case by taking the firing rates $r_{E,n}$ and $r_{I,n}$ in Eqs. (11), (12) or (14) expressed in terms of the respective currents by the stochastic relations Eq. (5) and Eq. (6) which account for the finite number of neurons in the networks.

Fig. 3 shows the close correspondence between the rate-model results and networks with a large number of neurons ($N \sim 10^6$) for $E \to E$ connectivity. The different dynamical regimes for two E-I modules coupled by long-range $E \to E$ connectivity are all observed in simulations of networks spiking neurons. Both antiphase (Fig. 3A), alternating dominance (Fig. 3B), modulated dominance (Fig. 3C) and full synchrony (Fig. 3D) are observed for parameters that quantitatively correspond to those predicted by the rate-model analysis. Similarly, in the case of $E \to E, I$ connectivity, two large E-I modules oscillate in full synchrony as predicted by the rate description (not shown).
<figure><img src="image_4.png" /><figcaption>Fig 3. Two large E-I modules coupled by long-range excitation. (1,2) E-I modules with a large number of neurons, $N = 1.6 \cdot 10^6$, and $E \to E$ connectivity. The different dynamical regimes depicted in Fig. 2 for different coupling strengths are clearly visible both in (1) spiking network simulations and in (2) the corresponding rate-model simulations. The respective auto- and crosscorrelations of the excitatory activities are shown in panels (3). The coupling strengths are as follows: (A) $f_{lr} = 0.01$ (light grey cross in Fig. 2B), (B) $f_{lr} = 0.015$, (C) $f_{lr} = 0.02$, (D) $f_{lr} = 0.05$.</figcaption></figure>

For smaller number of neurons, stochastic effects compete with deterministic effects both in simulations of spiking networks and in the rate-model description, when stochasticity arising from finite-size effects is accounted for. Interestingly, both descriptions continue to agree well when fluctuations play a major role in the dephasing of the two E-I modules, as shown in Fig. 4.

For $E \to E$ connectivity, the correlation between the activities of the two networks clearly displays the signature of the antiphase regime at weak coupling for $N = 2 \cdot 10^5$ neurons, whereas it is not apparent for $N = 10^4$ when stochasticity is stronger (Fig. 4A,B). At larger coupling, even small fluctuations lead to modules exchanging their roles and make more complex regimes difficult to identify even for $N = 2 \cdot 10^5$ neurons. This is illustrated by simulations with a coupling parameter in the modulated dominance regime ($f_{lr} = 0.02$, Fig. 4C). The ratio of cross- and auto-correlation functions at zero time lag shows that the transition between different regimes as a function of coupling is blurred by stochasticity as $N$ decreases (Fig. 4E,F).

The corresponding results for $E \to E, I$ connectivity are shown in Fig. 4G-L. In
<figure><img src="image_5.png" /><figcaption>**Fig 4.** Stochastic effects due to network size for two E-I modules coupled by long-range excitation. Auto- and crosscorrelations $C_{11}$ and $C_{12}$ of the activities of E-I modules with a finite number of neurons $N = 10^4$ and $2 \cdot 10^5$ for different coupling strengths $f_{lr} = 0.01$ and $0.02$ with long-range excitation targeting only excitatory populations (A-D) or targeting both excitatory and inhibitory neurons (G-J) (parameter values are stated in the panels' lower left corner). Summary plots with the ratio of equal-time crosscorrelations over autocorrelations for varying $N$ and $f_{lr}$ are shown in (E,F) and (K,L), respectively. In the $E \to E$ case, oscillatory synchrony between the two modules disappears when the coupling strength $f_{lr}$ diminishes ((C) vs. (D)), but the dynamical regimes of Fig. 2 are blurred by stochastic effects when the network size $N$ decreases ((A,B) vs. (C,D)). In the $E \to E, I$ case, no other dynamical regimes than oscillatory synchrony are expected, but synchrony increases both with coupling strength $f_{lr}$ and, more strongly, with network size $N$.</figcaption></figure>

this case, two effects can be observed. On the one hand, synchrony decreases at fixed coupling strength when $N$ is decreased due to larger noise (Fig. 4G,H). On the other hand, synchronization of the two networks increases when coupling strength increases at fixed $N$ (fixed noise) as shown in Fig. 4G,I and Fig. 4H,J. Both effects are quantified in Fig. 4K,L.

One should note that whereas the precise targeting of the long-range excitation plays a crucial role in the instability of in-phase oscillations for large networks ($N \sim 10^6$), it is much less important for smaller ($N \sim 10^4$) networks. For these smaller networks, the dominant role of fluctuations in the dynamics of the two individual E-I modules modifies the connectivity synchronization properties. It leads to similar synchronization
effects of long-range excitation, whether one module only targets the other E population ($E \to E$ connectivity) or targets both E and I populations in the other module ($E \to E, I$ connectivity) (compare Fig. 4B to Fig. 4H).

## Analysis of the competition between synchronization and noise

We thought it interesting to try and analyze more precisely the competition between synchronization and stochastic fluctuations, given that it arises for E-I modules of biologically realistic sizes. We chose to precisely study the case of $E \to E, I$ connectivity when stochastic fluctuations simply compete with full synchronization promoted by long-range connectivity. We developed two approaches which both assume that fluctuations are of an amplitude moderate enough that they can be treated as perturbations of the deterministic limit cycle.

A classical approach is available for weakly coupled modules and small noise. It was previously used to study noisy coupled oscillators [19] as well as neural systems [25]. In the present case, the outcome of this approach is that finite-size fluctuations in the rate-model description add a stochastic component to the previous Eq. (13) for the phase difference between the two modules (see Eqs. (78) and (84) in Methods). (It also provides an estimate to the small shift of oscillation frequency due to small noise [51, 52] which we do not consider here). This leads to a very explicit quantification of the difference of oscillation phases between the two modules induced by fluctuations,

$$ \langle (\phi_1 - \phi_2)^2 \rangle = \frac{D_N}{2 f_{lr} D_{\phi}^{EI}}. \qquad (16) $$

In this expression, the constant $D_N$ is simply the amplitude of stochastic fluctuations in a single module (Eq. (9)), whereas $f_{lr}D_{EI}$ quantifies the synchronization strength due to long-range connectivity (Eqs. (66) and (15)). As intuitively expected, the mean square phase difference between the two oscillators increases with the noise strength and decreases when the synchronization strength increases.

Eq. (16) provides a simple estimate of the competition between stochastic dephasing and synchronization. In order to compare it to numerical simulations, one can compute the respective phase of oscillations of two oscillatory E-I modules. Alternatively, Eq. (16) can be transformed into a prediction of the average spike-rate difference between the two modules, taking into account that the firing rate is a periodic function of the oscillatory phase (Eq. (82) in Methods). The latter is chosen to compare with results of numerical simulations of both FAT rate-model and spiking network implementations of two coupled E-I modules in Fig. 5. The mean square difference of excitatory rates between the two modules are quite comparable in the noisy FAT (Fig. 5A) and spiking networks (Fig. 5B), again showing the quantitative accuracy of the FAT rate-model description. However, Fig. 5 shows that the predicted rate difference fluctuations obtained from Eq. (16) (“phase approximation” in Fig. 5) are accurate only for very small coupling between the two modules.

The underlying assumption in Eq. (16) is that the phase difference plays a dominant role in the difference of activities between the two E-I modules because the phase difference mode has a much smaller restoring strength than the other (amplitude) mode between the two modules. The respective eigenvalues of the two modes are shown in Fig. 6G as a function of the fraction of long-range connection $f_{lr}$. Indeed, for weak coupling ($f_{lr} \ll 1$), one eigenvalue $\mu_1^A$ is close to 1, so that the restoring strength $1 - \mu_1^A$ of the corresponding (phase) mode is much smaller than $1 - \mu_2^A$, the restoring strength of the other mode. However, an explicit computation shows that $\mu_1^A$ quickly decreases as $f_{lr}$ increases (see Eq. (76) in Methods). Therefore, the two modes soon play comparable roles as $f_{lr}$ increases.
<figure><img src="image_6.png" /><figcaption>Fig 5. Competition between synchronization and noise for two E-I modules with $E \to E$, *I connectivity: theory vs. simulations*. The variance of excitatory activity between the two modules $\langle \delta r_E^2 \rangle = \langle (r_{1,E} - r_{2,E})^2 \rangle / 2$ is shown as a function of coupling strength $f_{lr}$. (A) Simulations of the stochastic rate model (FAT, duration 100 s, averaged over 3 repetitions) for different noise strengths (color symbols), as specified in the insert in term of the number of neurons *N* in the individual E-I module. These simulation results are plotted against theoretical predictions with different approximations. For very small coupling strengths, the phase approximation with the full limit cycle dynamics (Eqs. (16, 82), dotted lines) gives a good account of the observed synchronization, while for larger couplings both phase and amplitude modulations have to be considered (two-mode calculation, Eq. (97), solid black line). Note that for very small couplings, the two-mode calculation coincides with the phase approximation with linearized limit cycle dynamics (Eqs. (16,80), dashed black line), which does not take into account the periodicity of the firing rate as a function of the phase and which predicts an unphysical diverging activity difference for vanishing coupling. (B) Same as in (A) with the results for the network simulations of Fig. 4K,L (colored solid lines corresponding to different E-I module sizes, as specified in the inset), shown as $\langle \delta r_E^2 \rangle = C_{11}(0) - C_{12}(0)$.</figcaption></figure>

We developed another approach, valid for small noise but for arbitrary coupling $f_{lr}$, in order to check this origin of the inaccuracy of Eq. (16) as well as to obtain a quantitatively precise description of the competition between synchronization and stochasticity. The approach is simply based on accounting for the effect of stochastic perturbations at a linear level around the fully synchronized state where all modules are in the exact same states (see Methods: Competition between synchronization and stochasticity for two coupled E-I modules of finite size. Weak noise at arbitrary coupling). As shown in Fig. 5, this two-mode computation produces precise estimates of the mean square excitatory rate difference of the two modules, $\langle (r_{1,E} - r_{2,E})^2 \rangle$, even for modules with only $10^4$ neurons when stochasticity is quite high. This second approach is thus quite successful in quantitatively capturing the observed competition between synchronization and stochasticity.

## Phase gradients in a chain of oscillatory E-I modules coupled by long-range excitation

Having analyzed the synchronization dynamics and effects of finite-size fluctuations for two E-I modules, we now turn to a fuller description of spatially-structured connectivity. We consider a chain of E-I modules, with each module's long-range excitation targeting other modules in a distance-dependent way. Interestingly, the effects that we uncovered in the two-module setting, such as desynchronization at weak coupling for $E \to E$ connectivity and competition of synchronization with stochastic finite-size fluctuations, reappear in a more elaborate form. Moreover, the computational techniques that we developed for the two-module case turn out to be directly applicable to the chain of modules.
<figure><img src="image_7.png" /><figcaption>**Fig 6.** Spontaneous appearance of phase gradients in a chain of E-I modules. (A) Sketch of the chain of E-I modules coupled by long-range excitation. Excitatory neurons in one module target excitatory neurons ($E \to E$) or excitatory and inhibitory neurons ($E \to E, I$) in the other modules with a distance-dependent excitation profile. (B) Fourier transform of the coupling function $\tilde{C}_{\lambda}(q)$ (Eq. (106)) for perturbations of wavenumber $q$ when the long-range excitation has an exponential profile (Eq. (20)), for space decay constants $\lambda = 1$ (turquoise) and $\lambda = 0.33$ (dark blue). For $\lambda \ll 1$, the chain of modules can be approximated by a continuous system. (C) Effective long-range coupling $f_{lr}(q; \lambda) = (1 - \tilde{C}_{\lambda}(q))/2$ for perturbations of wavenumber $q$. (D) Stability spectrum for perturbations of wavelength $q$ as function of the effective coupling $f_{lr}(q; \lambda)$ for $E \to E$ connectivity. For small coupling strengths, i.e. long wavelengths, the largest eigenvalue (thick solid line) is larger than 1, and the perturbations are unstable ($f_{lr} < f_{lr}^*$, hatched region). While both eigenvalues are real and positive for small coupling strengths, they become complex conjugated for larger couplings (real part, thick dotted black line; imaginary part, thick dotted grey line). (E) Excitatory activity $r_E$ for a chain of deterministic E-I modules over time ($L = 512$, with 0.1 mV noise on initial condition). (F) Growth of the variance of excitatory activity $\sum_{n=1}^{L}(r_{n,E} - \sum_{m=1}^{L}r_{m,E}/L)^2$ over time, with the predicted rate for the fastest-growing mode shown on top (dashed black line). (G,H) Same as (D,E) but for $E \to E, I$ connectivity. Note that all modes are stable (G) and no instability develops over time (H). (I) Rate of decay of longest wavelength perturbation $\sigma_{2\pi/L}$ for different chain lengths $L$ (orange crosses) against the theoretical prediction (Eq. (23), black line). Inset: Decay of the spatial modulation of excitatory activity over time ($L = 512$). For all simulations shown, $\lambda = 1/3$.</figcaption></figure>

**Spontaneous appearance of phase gradients due to long-range excitation specifically targeting excitatory neurons**

The setting of our analysis is sketched in Fig. 6A. E-I modules with locally unstructured connectivity are coupled by long-range excitation. For $E \to E$ connectivity, the
dynamics of L coupled E-I modules are described in the rate-model framework (3,4) by

$$ \tau_E(I_{E,n}) \frac{dI_{E,n}}{dt} = -I_{E,n} + I_E^{\text{ext}} + w_{EE} \sum_{m=1, \dots, L} C(n,m) r_{E,m} - w_{EI} r_{I,n} \quad (17) $$

$$ \tau_I(I_{I,n}) \frac{dI_{I,n}}{dt} = -I_{I,n} + I_I^{\text{ext}} + w_{IE} r_{E,n} \quad (18) $$

with $C(n, m) = C(|n - m|)$.

For $E \to E, I$ connectivity, Eq. (18) is replaced by

$$ \tau_I(I_{I,n}) \frac{dI_{I,n}}{dt} = -I_{I,n} + I_I^{\text{ext}} + w_{IE} \sum_{m=1, \dots, L} C(n,m) r_{E,m} \quad (19) $$

The connection strength decreases monotonically with the distance $l$ between modules [17], as described by the function $C(l)$. Our analytical expressions are valid for a general functional form. For the numerical simulations, we have chosen, for illustration, an exponential decrease which appears compatible with experimental data in the primate motor cortex [16]. We specifically take

$$ C(l) = A_{\lambda} [\exp(-\lambda|l|) + \exp(-\lambda(L - |l|))] , \quad A_{\lambda} = \frac{1 - \exp(-\lambda)}{[1 - \exp(-\lambda L)][1 + \exp(-\lambda)]} \quad (20) $$

where $A_{\lambda}$ implements the normalization $\sum_{l=0, \dots, L-1} C(l) = 1$. The second exponential corresponds to choosing periodic boundary conditions for the chain. While this is not realistic for the cortex, it minimizes boundary effects in the simulations and should not significantly modify the results in the regime when $L$ is large compared to $1/\lambda$ (i.e., when $\exp(-\lambda L) \ll 1$).

We first analyze the dynamics in the framework of deterministic rate equations, before considering the effects of fluctuations and comparing with spiking networks in the next section.

A perfectly synchronized oscillating chain, with zero phase difference between the oscillations of different modules, is always a possible network state at the level of the rate-model description. It needs however to be tested whether full synchronization is resistant to perturbations or whether some perturbations are amplified by the dynamics and full synchronization is unstable. Interestingly, the linear stability analysis of a periodic lattice such as the considered 1D chain, exactly reduces to that of a two E-I module network with an effective coupling $f_{lr}$. More precisely, any small perturbation of the fully synchronized state is the linear sum of periodic perturbations with different wavenumbers $q$ which evolve independently of each other. Perturbations of wavenumber $q$ evolve in the same way as perturbations in the two E-I module network considered in the previous section, provided the coupling in the two-module network is chosen as

$$ f_{lr}(q) = [1 - \tilde{C}_{\lambda}(q)]/2, \quad (21) $$

where $\tilde{C}_{\lambda}(q)$ is the Fourier transform of the long-range excitation profile $C(l)$ (Eq. (106) in Methods). $\tilde{C}_{\lambda}(q)$ and the corresponding “effective” coupling in the two-module network, $f_{lr}(q)$, are displayed in Fig. 6B,C for an exponentially decreasing long-range excitation profile $C(l) \sim \exp(-\lambda|l|)$. It is seen in Fig. 6C that small wavenumbers $q$ correspond to weak effective coupling $f_{lr}(q)$, with $f_{lr}(q) \to 0$ when $q \to 0$. While deriving the exact relation (21) requires mathematical analysis (see Methods: Stability analysis of full synchronization for a chain of oscillatory E-I modules), it can be intuitively understood. For perturbations of small wavenumbers ($q \to 0$), only distant parts of the chain are in distinct oscillating states. The relative dynamics of these distant parts is
thus effectively only weakly coupled by excitation since long-range excitation decreases with distance. The relation (21) allows one to directly transcribe for the E-I module chain the results previously obtained for the two E-I module case. As before, we compare two connectivities for long-range excitation.

We first consider the case of $E \to E$ connectivity, when long-range excitation only targets excitatory neurons. The effective coupling $f_{lr}(q)$ decreases monotonically with the wavenumber $q$ and vanishes as $q \to 0$. The stability diagram displayed in Fig. 2B then shows that full in-phase synchronization is unstable for chains sufficiently long to accommodate modulations of small enough wavenumber $q$. The associated stability eigenvalues are shown in Fig. 6D. In-phase synchronization is unstable when $f_{lr}$ becomes smaller than a critical $f_{lr}^*$ that depends on the single E-I module synaptic parameters. For the chain of modules, this implies that spatially-periodic perturbations grow when their wavenumber $q$ lies below a critical $q^*$, defined by

$$ \cos(q^*) = \frac{1 - 2f_{lr}^* \cosh(\lambda)}{1 - 2f_{lr}^*}, \quad \text{or} \quad q^* \approx \sqrt{\frac{2f_{lr}^*}{1 - 2f_{lr}^*}} \lambda, \qquad (22) $$

where $\lambda$ controls the exponential decrease with distance of the long-range excitation. The second equality in Eq. (22) is valid for small $\lambda$; the used approximate relation between $f_{lr}$ and $q$ is shown in Fig. 6C, it can be seen that it is already very accurate for $\lambda = 1/3$.

Eq. (22) predicts that for a network of sufficient spatial extent $L > 2\pi/q^*$, spontaneous phase gradients should appear along the chain. This is confirmed by direct numerical simulations of the rate-model equations, shown in Fig. 6E. The exponential growth of the modulation of excitatory activity along the chain is shown in Fig. 6F, starting from a weakly perturbed fully synchronized state. This directly confirms the linear instability of full synchronization. One further notable feature in Fig. 6E is that the phase gradients develop on a spatial scale that is quite long (i.e. $\sim 30$ modules) compared to the spatial scale of the long-range excitation profile $C(l)$ ($\sim 1$ in Fig. 6D). This stems from the fact that in-phase synchronization is unstable only for small $f_{lr}(q)$, that is, for small $q$ or long wavelengths compared to the characteristic module coupling length. For our reference parameters and $w_{IE} = 2$ mVs, Eq. (22) gives for the critical wavelength $q^* \simeq 0.24\lambda$. Similarly, one obtains $q_m \simeq 0.15\lambda$ for the wavenumber $q_m$ of the fastest-growing mode, which corresponds to the largest eigenvalue in Fig. 6D.

For $E \to E, I$ connectivity, when long-range excitation targets both excitatory and inhibitory neurons, the previous two E-I module synchronization analysis shows that full synchronization is stable up to vanishingly low coupling $f_{lr}$. The stability eigenvalues are shown in Fig. 6G. For all wavelengths, their real part is less than 1. This linear analysis thus predicts that any modulation of activity along the chain should disappear exponentially. This is indeed what is observed in direct simulations of a chain of modules, as illustrated in Fig. 6H. Moreover, our stability analysis shows that the longest wavelength modes are the slowest to vanish. Their relaxation rate is obtained (Eq. (107) in Methods) as

$$ \sigma_q = \frac{D_{\phi}^{EI}}{2[\cosh(\lambda) - 1]} q^2. \qquad (23) $$

The relaxation rate $\sigma_q$ vanishes when $q \to 0$ and decreases with increasing modulation wavelength. The smallest wavenumber that can be accommodated in a chain of $L$ modules is $2\pi/L$. Fig. 6I shows that $\sigma_{2\pi/L}$ as given by Eq. (23) precisely describes the approach to full synchronization in our simulations of a E-I module chains of different lengths with long-range $E \to E, I$ connectivity.
<figure><img src="image_8.png" /><figcaption>Fig 7. (see next page)</figcaption></figure>

Phase gradients emerging from stochastic fluctuations in E-I modules of finite size

Spontaneous stochastic fluctuations of activity are present in each E-I module with a finite number of neurons, as already emphasized. In the two E-I module case, these were seen to compete with synchronization and promote phase differences between different E-I modules. Finite-size fluctuations similarly act in a chain of E-I modules. This is shown in Fig. 7, where the results of stochastic rate-model simulations are displayed for both $E \to E$ and $E \to E, I$ connectivity. Both the phases of oscillation of different modules (Fig. 7A,D) and the activity of excitatory neurons (7B,E) show that in both connectivity cases differences between neighboring modules are present. The histograms of phase differences between nearest modules are displayed in Fig. 7C,F for $E \to E$ and for $E \to E, I$ connectivity, respectively. In both cases, differences in phases decrease when the number of neurons in the modules increases, i.e., when stochastic fluctuations decrease. They also decrease when the range of long-range excitation increases and neighboring modules become more strongly coupled. The similarity between the two connectivity cases shows that stochastic fluctuations play the dominant role in creating phase differences between neighboring modules, for modules of size lower than $N = 10^5$.

The stochastic fluctuations along the two chains of E-I modules are further quantified in Fig. 7G,H by computing the power in different spatial modes. These clearly show that for both connectivities, power is highest in the longest wavelengths. This qualitatively
**Fig 7 (previous page).** Appearance of phase gradients and stochastic fluctuations in a chain of E-I modules. (A,B) Simulation of the stochastic rate model (FAT) for a chain with $E \rightarrow E$ connectivity, a space constant $\lambda = 0.33$ and network size $N = 10^4$. (A) Instantaneous phase as obtained from the Hilbert transform of bandpass-filtered excitatory activity, where the (unfiltered) rates are shown in (B). (C) Distributions of local phase gradients $\nabla\phi = \phi_{x+1} - \phi_x$ for different noise strengths and space constants for $E \rightarrow E$ connectivity. Inset: distributions of gradients normalized to the respective mean value. (D,E) Phase and excitatory activity for the same parameters as in (A,B) but with $E \rightarrow E, I$ connectivity. (F) Same as in (C) but for $E \rightarrow E, I$ connectivity. (G,H) Power spectra of the spatial activity profiles (averaged over 5 s) for different module sizes and space constants $\lambda$, for $E \rightarrow E$ (G) and $E \rightarrow E, I$ (H) activity, respectively. The spectra obtained from spiking network simulations ($N = 5000: L = 128; N = 10^4: L = 32$) are shown in light thick lines, spectra obtained from rate-model simulations are shown in dark thin lines ($L = 256$). (I) Comparison of the spectra obtained from simulations (FAT, solid lines) and theoretical predictions for $E \rightarrow E, I$ connectivity. Similar to the two-module case, the modes with larger effective couplings are well described by the two-mode calculation (dashed lines, Eq. (136)), while for smaller couplings, or vanishing $q$, non-linear effects due to the periodicity of the limit cycle, which are not captured by the linearized two-mode dynamics, become dominant. This limit is better captured by the phase approximation and taking the periodicity of the limit cycle dynamics into account (dash-dotted lines, Eq. (123)). Network size $N$ and space constant $\lambda$ as in (G,H). (J) Example of a spatially extended phase gradient (upper panel) that gives rise to a the apparent propagation from left to right of a (negative) bump of activity in space (lower panel). Phase and excitatory firing rate are shown over a length of 15 modules for different instances in time, see color code legend. (K) Correlation of spatially extended phase gradients $\overline{\nabla\phi}$ (averaged over 5 modules and retained when the absolute average gradient is larger than the standard deviation of local phase differences) with phase propagation velocities $v_\phi$ determined from frame-to-frame correlations between spatial phase profiles (see Methods). (L) Summary plot of average phase gradients for different network sizes $N$ and connection decay space constants $\lambda$ for $E \rightarrow E$ connectivity (circles) and $E \rightarrow E, I$ connectivity (diamonds).

agrees with a generalization of the two-module result Eq. (16) for the chain. This result, given in Eq. (127), is obtained by assuming that oscillatory phase differences dominate the fluctuations between modules. While this is true for very long wavelengths, it does not quantitatively agree with the measured spectra for the smaller wavelengths of interest in the present context (Fig. 7I). A generalization of the two-module computation that we developed to take into account both phase and amplitude modes (see Methods: Characterization of the E-I module chain stochastic activity profile) is required to quantitatively describe the simulation results, which is also shown in Fig. 7I.

The above results were obtained by simulating chains of E-I modules using the FAT rate-model description. Equivalent simulations for spiking networks are too demanding in computational power for us to perform. We have however checked that for shorter chains of modules ($L = 128$ for $N = 5000$, $L = 32$ for $N = 10^4$), the stochastic FAT rate model reproduces well the results of spiking network simulations (Fig. 7G,H). This renders us confident that this would be the case for longer chains as well.

## Phase waves and propagation of activity

We have shown that stochastic fluctuations in a chain of E-I modules create phase differences in the oscillations of E-I modules along the chain. In the presence of phase differences, activity in different modules peaks at different times and may produce “phase waves” (also described as traveling waves, e.g. in [2, p1219]), as shown in Fig. 7J. The propagation velocity $v_\phi$ of phase waves is kinematically related to the phase gradient
and the distance between E-I modules (which is also the linear size of an E-I module),

$$v_{\phi} = \frac{d}{\nabla\phi}, \qquad (24)$$

where $\nabla\phi$ denotes the phase variation between neighboring modules (recall that the phase has the unit of time with our convention) and $d$ is the physical distance between these module centers.

These transient phase waves bear some resemblance to the observation of traveling waves in the motor cortex during episodes of beta oscillations [10–14]. Thus, it appears of interest to better characterize them in the present setting and to compare them to propagation of activity. Phase waves covering a notable spatial range correspond to phase gradients that extend over several modules. Although phase differences between modules are mostly correlated on short length scales, extended phase gradients nevertheless exist. For our purpose, we determined gradients of activity that extend over five modules (see Methods: Simulations and numerical analysis. Detection of propagation events in an E-I module chain), and computed the corresponding phase wave velocity. Furthermore, we independently assessed the spatio-temporal propagation of activity for these occurrences by correlating successive frames of simulations (see Methods: Simulations and numerical analysis. Detection of propagation events in an E-I module chain). These two independently measured velocities are strongly correlated as shown in Fig. 7K. In our simulations, phase gradients and the associated phase waves thus explain the detected propagation of neural activities.

The probability distributions of spatially extended phase gradients are shown in Fig. S4, for the two different considered connectivities and two space decay constants $\lambda$ of long-range excitation, as well as for different numbers of neurons in the E-I modules. The corresponding mean values of these gradients are summarized in Fig. 7L. In order to compare these values to neural recordings, values for the size of an E-I module and the range of long-range horizontal connections should be determined. Data are available for different species and for different areas [16, 17, 59–61]. For instance, for mouse somatosensory barrel cortex, identifying an E-I module with a barrel column would give $d \simeq 175 \mu m$ for the distance between neighboring modules (which is also the linear size of an E-I module) and $\lambda \simeq 0.7$ for layer 2/3 horizontal connections since synaptic excitatory charge is measured to decay to 26% two columns away from an excited column [60].

For the motor cortex of adult monkeys, data are available from intra-cortical micro stimulation and recordings with single and multielectrodes [16]. Response amplitudes vary with the protocol and exact location of the electrodes, but a space decay constant of 1.5 mm seems representative. An E-I module size of $d \simeq 500 \mu m$ and $N \simeq 2 \cdot 10^4$ neurons then gives $\lambda \simeq 0.33$ for the long-range connectivity. For this last case, our numerical results give $\nabla\phi \simeq 0.5$ ms ($E \to E, I$) or 1 ms ($E \to E$) depending on the precise connectivity (Fig. 7K). With the above estimate of $d$, Eq. (24) then gives a phase velocity $v_{\phi} \sim 1$ m/s or $0.5$ m/s, respectively. In both cases, this is significantly larger than the beta wave velocities recorded in the motor cortex of monkeys. Increasing the local noise on each module increases the dephasing between modules, e.g. $\nabla\phi \simeq 1.4$ ms for $N = 5 \cdot 10^3$ and ($E \to E, I$) connectivity (Fig. 7K), and correspondingly lowers $v_{\phi}$. Thus, in this framework, stochastic dephasing between modules could account for the experimental results only if other sources of “noise” are present beyond the one coming from the finite size of local modules. This would suggest that different E-I modules receive different inputs, shared at the level of individual modules but varying independently from module to module, as further discussed below.

We assessed one other possibility of increasing local noise with respect to the long-range synchronization of neighboring modules. It would consist of diminishing the importance of long-range coupling relative to purely local excitation, without modifying
the characteristic length scale of long-range connectivity as such. The associated coupling function then takes the form $C(l) \sim (1 - \alpha)\delta_{0l} + \alpha \exp(-\lambda|l|)$, where $\alpha$ measures the fraction of truly long-range connectivity. For $\alpha = 0.5$, the resulting average phase gradients are about two times steeper for all $\lambda$ and $N$, see Fig. S5, implying phase waves that are about two times slower than for $\alpha = 1$ (compare Fig. S5E to Fig. 7L).

While noise increase brings the phase wave speed in simulations closer to the range of recorded ones, other features may be responsible for the quantitative differences between the present simplified model and biological recordings, as discussed below.

## Discussion and conclusion

In the present work, we have studied the synchronization of spatially distant E-I modules connected by long-range excitation.

At a technical level, our analysis was greatly facilitated by the formulation of a stochastic rate model with an adaptive timescale that accounts quite precisely for the activity of spiking networks. This extends the finding of previous works [48, 50] to E-I oscillations. The use of a quantitatively precise rate model offers two main advantages. First, it eases numerical simulations that are computationally very demanding for spiking networks with spatially-structured connectivity, since they require both a significant amount of locally connected neurons and a finite number of such modules for spatial extension. Second, the rate-model formulation lends itself much more easily to mathematical analysis [21, 24, 25, 27] than its spiking network counterpart. Therefore, the approach may certainly prove useful in other contexts.

Our application of this methodology provides a description of spatially-structured networks with a local oscillatory dynamics in the SSO regime, based on E-I recurrent interactions. Our first general finding is that, as previously obtained in the Wilson-Cowan rate model framework [21, 22, 27], the detail of long-range connectivity matters. Based on observations in several neural areas [15–17, 61, 62], we have confined ourselves to considering long-range excitatory connections. Long-range excitatory connections only targeting distant excitatory neurons result in more complex dynamical properties than long-range excitatory connections that have the same connection specificity as local ones. Experimental investigations of the targeting specificity of long-range excitatory connections appear rather scarce at present (but see [60, 63]). Results such as ours and previous ones [21, 22, 27] will hopefully provide an incentive to gather further data on this question for which experimental tools are now available.

Complex dynamical regimes for two spiking oscillatory modules connected by long-range excitations were previously found [34] for fast oscillations driven by recurrent connections between inhibitory interneurons with synaptic delays, with delays also being present in the excitatory connections. Our results, in agreement with previous rate model calculations [22] extend these previous findings to a different oscillatory regime, for slower beta oscillations, in the absence of recurrent inhibitory connections and significant synaptic delays. It will be interesting to see if these complex regimes are observed in experiments like the one performed in ref. [9] which record and perturb distant neural ensembles oscillating in the beta/low gamma range. They may be also be relevant for oscillations in large ensembles of neurons coupled by long-distance connections, e.g. inter-areal dynamics [64].

We have however observed that these complex oscillatory regimes for two coupled E-I modules are very sensitive to stochastic fluctuations of activity. This has led us to develop a quantitative description of the competition between finite-size fluctuations and synchronization in order to describe the dynamics at the scale of a cortical column. Other studies have proposed more refined and complex ways to take into account finite-size fluctuations (see e.g. [65–67]) than the one we used. But, at least for our purposes, the
simple procedure of ref. [39] proved sufficient to account for spiking network simulations.

We have also found that the instability of full synchronization for the two-module network translates in a long-distance instability of synchronization in the extended spatial setting of a chain of modules when long-distance excitation primarily targets excitatory pyramidal cells ($E \rightarrow E$ connectivity). This instability is analogous to those occurring in oscillatory media such as the one created by an oscillatory chemical reaction taking place in a layer of liquid (see [68] for a review). It can similarly be captured by considering the spatio-temporal dynamics of the local oscillation phases as originally proposed in a general coupled-oscillator context [19] and extended in previous works to neural systems [21–25].

This phase instability results in the spontaneous creation of phase differences between neighboring spatial locations, namely oscillatory phase gradients. Phase gradients kinematically lead to the observation of traveling waves of activity (see [2, 23] for previous discussions) and seem at the origin of traveling waves associated with oscillatory activity observed in different contexts. The instability wavelength that we determined appears too long however to account for observations in motor cortex [10–14] taking a cortical column typical size of 500 µm and available data on the range of excitatory connections [15–17]. Moreover, this would require that $E \rightarrow E$ connectivity prevails in motor cortex, which is certainly not clear at present.

Stochasticity arising from the finite number of neurons in a cortical column plays a role as important in a spatial setting as for two modules. It competes with the synchronizing effect of connectivity in the $E \rightarrow E, I$ case or in the $E \rightarrow E$ case at length scales smaller than the phase instability critical wavelength. This competition itself induces phase differences between spatially neighboring neural oscillatory ensembles. These produce the stochastic and transient appearance of phase waves, as we have shown (Fig. 7J,K). This has led us to search and obtain a quantitatively precise understanding of stochastic dephasing between different modules. Phase differences are most simply described by stochastic equations for the spatially varying oscillation phase that we derived, and that have been intensely studied in other contexts [69]. We have found however that the stochastic phase description accurately describes oscillatory phase differences only for sufficiently distant spatial locations which appear beyond those relevant for waves within a neural area. On shorter distances, relevant for intra-areal dynamics, modes beyond the phase modes have to be included to obtain a quantitatively accurate description of neural activity and its fluctuations. Dephasing of oscillators in space and the allied phase waves account for the propagating activity events in our simulations (Fig. 7J, K). However, the measured phase wave velocities appear too high to account for the measured velocities [10, 12, 13] of the transiently appearing traveling waves in the motor cortex during beta episodes. This points to the existence of further effects beyond the ones we have taken into account in the simple model that we have analyzed here. First, there may exist sources of dephasing between modules beyond finite-size noise. It seems plausible that this may arise from inputs from other neural areas which vary from module to module. This seems consistent with the topographic connectivity between the motor cortex and body muscles [15]. Whether these inputs can be treated as white noise or are better approximated by noise with a significant correlation time remains to be seen. Another hypothetic origin of steeper phase gradients suggested by our model would be an overrepresentation of local connections relative to the long-range excitatory connections studied here. This question will have to be addressed quantitatively in light of precise data that will eventually be available. A third non-exclusive possibility is that the slow conduction velocity along non-myelinated axons [59, 62] and the associated propagation delays amplify the phase gradient between different modules. The present analysis provides a solid foundation to incorporate these effects that we plan to examine in forthcoming work together with a detailed comparison
with neural recordings.

An interesting further aspect of our results is that the correlation of neural activity at different spatial points is related to the profile of long-range connectivity via the effective coupling. This relation could certainly be tested when the appropriate data becomes available. Interestingly, it could also serve to deduce profiles of long-range connectivity from measurements of neural activity, i.e., measured correlations between spatially distant LFP recordings.

Several other directions will be interesting to pursue to extend the present study. We have here restricted ourselves to consider a one-dimensional chain of E-I modules. Recordings not only display planar waves but also rotating waves and more complex wave patterns. Extension to two dimensions is required for a precise comparison with experimental data [12, 13]. Recorded planar waves also display preferred propagation axes [10] that it will be interesting to link to anisotropies in connectivity or possibly to geometrical properties of the motor and premotor cortices. Finally, we have analyzed a very simple E-I module. It will be important to decipher the roles of the different interneuron types (cf. [9]) of the different cortical layers and of their connectivity [60] to better understand cortical dynamics.

# Methods

## Simulations and numerical analysis

### Simulations of spiking networks

For our spiking network simulations, we used exponential integrate-and-fire (EIF) neurons that exhibit the following membrane potential dynamics:

$$ \tau_m \frac{d}{dt} V_i = E_L - V_i + \Delta_T e^{\frac{V_i - V_T}{\Delta_T}} + I_{ext,i} + I_{syn,i}, \quad (25) $$

where $\tau_m = 10$ ms is the membrane time constant, $E_L = -65$ mV is the resting potential set by a leak current, and $\Delta_T = 3.5$ mV and $V_T = -59.9$ mV are the sharpness of the spike onset and the spike threshold, respectively. The external current $I_{ext}$ represents inputs from neurons outside the network and comprises a mean current shared among all neurons of the excitatory or inhibitory population and a noise term specific to each neuron,

$$ I_{ext,i} = I_X^{ext} + \sigma_X \sqrt{\tau_m} \xi_i, \quad (26) $$

where $\xi_i$ is a Gaussian white noise with zero mean and unit variance, $\langle \xi_i(t) \xi_j(t') \rangle = \delta_{ij}\delta(t - t')$, and $X = E/I$ depending on whether neuron $i$ is excitatory or inhibitory. Synaptic currents $I_{syn,i}$ are given by

$$ I_{syn,i} = \tau_m \sum_j J_{ij} S_j(t), \quad (27) $$

where the $J_{ij}$ are the synaptic weights between post- and presynaptic neurons $i$ and $j$, respectively, and $S_j(t) = \sum_k \delta(t - t_k^{(j)})$ is the spike train of neuron $j$ with spike times $t_k^{(j)}$.

If not specified otherwise, we considered an all-to-all connectivity at the level of an individual module, such that $J_{ij} \in \{J_{EE}, J_{IE}, J_{EI}, J_{II}\}$ depending on the type (excitatory vs. inhibitory) of neurons $i$ and $j$. The synaptic weights $J_{XY}$ relate to population weights $w_{XY}$ according to $J_{XY} = w_{XY}/(N_Y \tau_m)$. Here, $N_E$ and $N_I$ are the numbers of excitatory and inhibitory neurons in a module; for a given network size $N$, $N_E = 0.8N$ and $N_I = 0.2N$.
We chose the external currents $I_E^{ext}$, $I_I^{ext}$ such that at steady state, excitatory and inhibitory neurons fire at an average rate of $r_E^{(s)} = 5$ Hz and $r_I^{(s)} = 10$ Hz, respectively. We furthermore fixed the noise strengths $\sigma_E, \sigma_I$ by imposing a total noise (external plus synaptic) of $\sigma_{tot} = 10$ mV, using $\sigma_{tot}^2 = \sigma_X^2 + J_{XE}^2 N_E \tau_m r_E^{(s)} + J_{XI}^2 N_I \tau_m r_I^{(s)}$, $X \in \{E, I\}$. Further parameters were: the numerical threshold potential after which a spike was registered and the membrane potential reset, $V_{thresh} = -30$ mV; the reset potential, $V_{reset} = -68$ mV; and a refractory period of $\tau_{ref} = 1.7$ ms. We used a time step $dt = 0.01$ ms throughout.

For reasons of computational efficiency, we used the following connectivity scheme for the two-module simulations. Because all-to-all connectivity can be implemented without actually processing spikes over all $O(N^2)$ synapses, we maintained all-to-all connectivity with lowered synaptic weights $\tilde{J}_{XE} = (1 - f_{lr}) J_{XE}$ for those synaptic connections within a given module that are supplemented by long-range excitation ($X = E$ in the $E \to E$, $X = E, I$ in the $E \to E, I$ connectivity case, respectively). The long-range connections were implemented as a finite number $C_{lr} = f_{lr} N_E$ (rounded to the nearest integer) of excitatory synaptic connections across modules with unchanged synaptic weights $J_{XE}$. Only for the very large networks ($N = 1.6 \cdot 10^6$) used in Fig. 3, we also replaced the connections across modules with all-to-all connectivity, where we reduced the synaptic weights accordingly to $\tilde{J}_{XE} = f_{lr} J_{XE}$.

The long-range excitatory connections between modules in the chain simulations were set up as follows: We first determined the distance-dependent number of excitatory connections as a function of $N_E$ and $\lambda$ such that the $N_E$ input synapses are distributed according to their relative contribution $p(\Delta x) = e^{-\lambda|\Delta x|}/(1+2/(e^{\lambda}-1))$ across modules, using periodic boundary conditions so that each neuron could in principle receive input from $2L-1$ modules. Input connections from different modules were then randomly drawn in corresponding numbers for each neuron receiving long-range excitatory connections.

All network simulations were performed using brian2, a convenient Python package for implementing network simulations [70].

Population firing rates were calculated as the total number of spikes within a time bin $\Delta t$ divided by the size of the population; we typically used $\Delta t = 0.1$ ms if not explicitly stated otherwise.

## Computation of correlation functions

The auto- and cross-correlation functions $C_{1j}(\tau)$ with $j = 1, 2$, respectively, as shown in Figures 1-4, were calculated as

$$C_{1j}(\tau) = \frac{1}{M} \sum_{t_i} \delta r_{1,E}(t_i) \delta r_{j,E}(t_i + \tau), \qquad (28)$$

with $M$ the number of discrete time points $t_i$ considered, and $\delta r(t_i) = r(t_i) - \bar{r}$ the deviatory part of the firing rate relative to its temporal mean. Note that $M$ was limited by the duration of the simulation and we generally discarded an initial transient of 250 ms. We furthermore subtracted the expected zero-time Poisson contribution due to stochastic firing, $\langle \delta r_E(t)^2 \rangle_{Poisson} = \bar{r}_E/(N_E \Delta t)$, from the autocorrelation at $\tau = 0$ to focus on the correlations induced by the network dynamics.

## Detection of propagation events in an E-I module chain

We first resampled the firing rates obtained by simulations with a new time step $\Delta t = 1$ ms by averaging the activity in 1 ms bins to reduce the sampling noise $\propto 1/dt$. We then determined the Hilbert phase from the band-pass filtered version of the deviatory firing rate $\delta r = r - \bar{r}$, using a fifth-order Butterworth filter with lower and upper corner
frequencies of 15 and 40 Hz, respectively. The Hilbert transform was obtained using the Python scipy package.

In a next step, we determined locally averaged phase gradients by averaging phase differences $\nabla\phi_x(t_i) = \phi_{x+1}(t_i) - \phi_x(t_i)$ over five neighboring segments, i.e., $\overline{\nabla\phi_x(t)} = \frac{1}{5}\sum_{k=0}^4 \nabla\phi_{x+k}(t_i) = (\phi_{x+5} - \phi_x)/5$ (note that all phase differences are taken modulo a full oscillation period). For the detection of propagation events, we retained gradients where the absolute average phase gradient $|\overline{\nabla\phi_x(t_i)}|$ is larger than the standard deviation of the five contributing local phase differences. To determine the propagation of the local phase profile $\phi_x(t_i) = (\phi_x(t_i), \phi_{x+1}(t_i), \dots, \phi_{x+4}(t_i))$, we calculated the Euler distance of laterally shifted versions at later time bins $t_{i+j}$,

$$d(\Delta, j; x, t_i) = \|\phi_x(t_i) - \phi_{x+\Delta}(t_{i+j})\|^2.$$

We subsequently inferred the phase propagation speed $v_\phi(x, t_i)$ from the optimal shift $\Delta_{opt}(j; x, t_i)$ that minimizes the distance between phase profiles after a given time $j\Delta t$ according to:

$$v_{\phi}(x, t_i) = \frac{1}{5} \sum_{j=1}^{5} \frac{\Delta_{opt}(j; x, t_i)}{j \Delta t}.$$

For reasons of computational efficiency, we restricted the search for an optimal shift $\Delta$ over shifts of $\pm 15$ modules.

## Fit of the adaptive timescale in the FAT rate model

In the adaptive rate model (2), the response time depends on the current $I$. In a previous work [48], Ostojic and Brunel (OB) proposed the analytically-motivated choice,

$$\tau^{(OB)}(I) = \tau_m \Delta_T \frac{\Phi'_{\sigma}(I)}{\Phi_{\sigma}(I)}.$$

In the present paper, the adaptive response time was chosen following [50], by taking the exponential kernel that best fitted the firing-rate response to an oscillating current at different frequencies. Precisely, we numerically evaluated the analytical expression for the linear firing-rate response $\hat{R}_1(f)$ of EIF neurons [71] for frequencies from 1 Hz to 1 kHz in steps of 1 Hz for all values of the input current on the grid of input currents $I$ with step size $dI = 0.1$ mV, and fitted the modulus of the firing-rate response by $A/\sqrt{1 + (2\pi f \tau_{fit})^2}$, which is the Fourier transformation of an exponential decay in the time domain (see Fig. S1C). We used these values of $\tau_{fit}$ as a look-up table and interpolated in between. The resulting “fitted adaptive timescale” (FAT) curve $\tau^{(FAT)}(I)$ is shown in Fig. S1B together with $\tau^{(OB)}(I)$. We also used look-up tables and interpolation for $\Phi_\sigma(I)$ and the first and second derivatives of $\Phi_\sigma(I)$ to speed up the numerical calculations. The tabulated functions $\Phi_\sigma(I)$ and $\tau^{(FAT)}(I)$ used for the present paper are provided as supporting material S1 File.

## Simulations of the FAT rate-model for an E-I module

For a network with synaptic weights $w_{EE}$, $w_{IE}$ and $w_{EI}$, the currents $I_E^{ext}$, $I_I^{ext}$ are chosen so as to impose the steady firing rates $r_E^{(s)}$, $r_I^{(s)}$ for the excitatory and inhibitory populations. Namely,

$$I_E^{ext} = I_E^{(s)} - w_{EE} r_E^{(s)} + w_{EI} r_I^{(s)},$$

$$I_I^{ext} = I_I^{(s)} - w_{IE} r_E^{(s)},$$
where $I_E^{(s)}$ and $I_I^{(s)}$ are the currents corresponding to the chosen rates according to
$r_E^{(s)} = \Phi_\sigma(I_E^{(s)}), r_I^{(s)} = \Phi_\sigma(I_I^{(s)})$. The same external currents are used in the spiking
network simulations.

We implemented the stochastic rate-model equations in Python using an Euler scheme with a time step $dt = 0.01$ ms. At each time $t_i$, we calculated the noisy firing rate $r_X(t_i)$ by drawing a number of emitted spikes $n_X(t_i)$ from a Poisson distribution with mean $N_X \Phi[I_X(t_i)]dt$, and setting $r_X(t_i) = n_X(t_i)/(N_X dt)$, $X = E, I$, in order to avoid negative firing rates.

## Oscillatory instability for an E-I module

The E-I module (3,4) either has its E and I populations stably firing at the chosen rates $r_E^{(s)}, r_I^{(s)}$ (i.e., the chosen fixed point is stable) or departures from $r_E^{(s)}, r_I^{(s)}$ are amplified and the rates cannot be maintained (i.e., the fixed point is unstable). These different cases are shown in parameter space in Fig. 1A. The regions of Fig. 1A are obtained in a standard way [20, 24] by linearizing Eqs. (3,4) around the fixed point currents $\mathbf{I}^{(s)} = (I_E^{(s)}, I_I^{(s)})$. Using a vector notation for the currents, $\mathbf{I} = (I_E, I_I)$, the perturbation $\delta\mathbf{I} = \mathbf{I} - \mathbf{I}^{(s)}$ evolves according to

$$ \frac{d}{dt} \delta \mathbf{I} = L_{EI} \cdot \delta \mathbf{I}, \qquad (29) $$

where matrix multiplication is denoted by a dot and $L_{EI}$ is the 2 × 2 linear stability matrix

$$ L_{EI} = \begin{pmatrix} \frac{[-1 + w_{EE} \Phi'_{E}(I_E^{(s)})]}{\tau_E(I_E^{(s)})} & -\frac{w_{EI} \Phi'_{I}(I_I^{(s)})}{\tau_E(I_E^{(s)})} \\ w_{IE} \Phi'_{E}(I_E^{(s)}) & -\frac{1}{\tau_I(I_I^{(s)})} \end{pmatrix}. \quad (30) $$

One can note that the derivatives $\tau'_E(I_E^{(s)})$, $\tau'_I(I_I^{(s)})$ do not appear at the linear level since these timescales multiply time derivatives which are themselves small quantities at the linear level. Thus, the adaptive character of $\tau_E$ and $\tau_I$ in the FAT model does not modify the computation.

The solutions of the linear system (29) are linear combinations of two exponentials in time, $\exp(\kappa t)$, with the two arguments $\kappa$ being solutions of

$$ \tau_E \tau_I \kappa^2 + [\tau_E + \tau_I(1-\alpha)]\kappa + 1 - \alpha + \beta = 0. \qquad (31) $$

Here, we have defined the two positive constants

$$ \alpha = \Phi'_{E}(I_{E}^{(s)}) w_{EE} \quad \text{and} \quad \beta = \Phi'_{I}(I_{I}^{(s)}) \Phi'_{E}(I_{E}^{(s)}) w_{EI} w_{IE}. \qquad (32) $$

The conditions for the stability of the steady state are that the real parts of the two roots of Eq. (31) are negative, namely that the sum of the two roots is negative and their product is positive, respectively

$$ \alpha < 1 + \frac{\tau_E}{\tau_I}, \qquad (33) $$

$$ \beta > \alpha - 1. \qquad (34) $$

If the latter condition (34) on the product of the two roots is violated, the two solutions are real and of different signs, the instability is purely exponential. When condition (34) holds but not the one on the sum of the two roots (33), the instability is either purely
exponential or oscillatory depending on how strongly the second condition is violated. It is oscillatory when $\alpha - 1$ is not too large, namely in the domain

$$ \alpha > 1 + \frac{\tau_E}{\tau_I}, \qquad (35) $$

$$ \beta > \frac{\tau_I}{4\tau_E}(\alpha - 1)^2 + \frac{1}{2}(\alpha - 1) + \frac{\tau_E}{4\tau_I} = \frac{\tau_I}{4\tau_E}\left(\alpha - 1 + \frac{\tau_E}{\tau_I}\right)^2. \qquad (36) $$

## Phase dynamics of a forced oscillatory E-I module

We consider an E-I module in the oscillatory regime forced by a small time-dependent current, as described by Eqs. (3,4) with

$$ I_E^{\text{ext}}(t) = I_{E,0}^{\text{ext}} + I_E^{\text{forc}}(t), \quad I_I^{\text{ext}}(t) = I_{I,0}^{\text{ext}} + I_I^{\text{forc}}(t). \qquad (37) $$

Different choices of the forcing currents, $I_E^{\text{forc}}(t)$ and $I_I^{\text{forc}}(t)$, allow us to describe the influence of finite-size stochastic effects as well as the synchronization dynamics of coupled E-I modules. The computation of the drift and entrainment of oscillators is a classical topic [18]. Methods to compute the effect of a weak forcing are well-known [19] and have been applied to neural rate models in numerous works (see [23–25] for reviews). We provide a derivation of these classical formulas for the FAT rate model, for the convenience of the reader.

As above (e.g. Eq. (29)), we adopt a vector notation for the currents. We suppose that the synaptic parameters are such that, with the current $I_0^{\text{ext}}$ alone, the E-I module follows the limit cycle $\mathbf{I}^{(0)}(t+\phi) = (I_E^{(0)}(t+\phi), I_I^{(0)}(t+\phi))$ of period $T$ (i.e., $\mathbf{I}^{(0)}(t+T) = \mathbf{I}^{(0)}(t)$). Here, the phase $\phi$ is an arbitrary constant number reflecting the invariance of the dynamics under time translation when $\mathbf{I}^{\text{ext}}$ is constant and there is no time-dependent forcing. When perturbed by a small time-dependent current $\mathbf{I}^{\text{forc}}(t)$, the E-I module currents $\mathbf{I}(t)$ closely follow the unperturbed limit cycle $\mathbf{I}^{(0)}(t+\phi)$. However, the phase $\phi$ of the followed limit cycle slowly drifts in time, since it is not submitted to any restoring force.

Note that, for simplicity, in the present work, we define the phase to be a periodic variable between 0 and $T$, instead of normalizing it to $2\pi$ as it is sometimes done.

In order to determine how $\mathbf{I}^{\text{forc}}(t)$ entrains the phase $\phi$, we compute perturbatively the correction $\delta\mathbf{I}(t+\phi)$ to the zeroth-order evolution $\mathbf{I}^{(0)}(t+\phi)$. From Eqs. (3) and (4), we obtain for $\delta\mathbf{I}$

$$ \frac{d\delta\mathbf{I}(\zeta)}{d\zeta} = \mathcal{L}[\mathbf{I}^{(0)}(\zeta)] \cdot \delta\mathbf{I}(\zeta) + \mathbf{F}(\zeta) \qquad (38) $$

with $\zeta = t + \phi$ and

$$ \mathbf{F}(t+\phi) = \mathbf{F}^{\text{forc}}(t+\phi) - \frac{d\mathbf{I}^{(0)}(\zeta)}{d\zeta} \frac{d\phi}{dt}, \qquad (39) $$

$$ \mathbf{F}^{\text{forc}}(t) = \left( \frac{I_E^{\text{forc}}}{\tau_E[I_E^{(0)}(t)]}, \frac{I_I^{\text{forc}}}{\tau_I[I_I^{(0)}(t)]} \right). \qquad (40) $$

The matrix $\mathcal{L}$ describes the linearized dynamics around the limit cycle,

$$ \mathcal{L}[\mathbf{I}^{(0)}(t)] = \begin{pmatrix} \frac{-1+w_{EE} \Phi'_{E}[I_E^{(0)}(t)] - \tau_E[I_E^{(0)}(t)] I_E^{(0)}(t)}{\tau_E[I_E^{(0)}(t)]} & -\frac{w_{EI} \Phi'_{I}[I_I^{(0)}(t)]}{\tau_E[I_E^{(0)}(t)]} \\ \frac{w_{IE} \Phi'_{E}[I_E^{(0)}(t)]}{\tau_I[I_I^{(0)}(t)]} & -\frac{1+\tau'_I[I_I^{(0)}(t)] I_I^{(0)}(t)}{\tau_I[I_I^{(0)}(t)]} \end{pmatrix}. \qquad (41) $$

The last term in Eq. (39), proportional to the time derivative of the phase $\phi$, comes from the slow phase evolution $\phi = \phi(t)$ that will be required for preventing the appearance of growing secular terms [72], as explained below.
Eq. (38) is a linear equation with coefficients which are periodic functions of time. For this Floquet problem [73], it is helpful to consider the evolution of $\delta\mathbf{I}(t)$ from one period to the next. For this purpose, we introduce the so-called monodromy matrix $M(t)$ such that

$$ \frac{d}{dt} M = \mathcal{L}[\mathbf{I}^{(0)}(t)] \cdot M \qquad (42) $$

with $M(0) = \mathbf{I}_d$, the identity matrix. The matrix $M(t)$ describes the linear evolution along the limit cycle and $M(T)$ the evolution after one full turn around the limit cycle. The vector tangent to the limit cycle $\mathbf{e}_1 = (I_E^{0'}(0), I_I^{0'}(0))$ is an eigenvector of $M(T)$ with eigenvalue $\mu_1 = 1$. The other eigenvector of $M(T)$, $\mathbf{e}_2$, has an eigenvalue $\mu_2$ with $|\mu_2| < 1$ since the limit cycle is stable. Eq. (38) can be explicitly integrated by writing $\delta\mathbf{I}(t)$ and $\mathbf{F}(t)$ in the basis obtained by translating $\mathbf{e}_1, \mathbf{e}_2$ along the limit cycle,

$$ \delta\mathbf{I}(t) = x_1(t) \mathbf{e}_1(t) + x_2(t) \mathbf{e}_2(t), \qquad (43) $$

$$ \mathbf{F}(t) = f_1(t) \mathbf{e}_1(t) + f_2(t) \mathbf{e}_2(t) \qquad (44) $$

with

$$ \mathbf{e}_j(t) = M(t) \mathbf{e}_j, \quad j = 1, 2. \qquad (45) $$

Substitution in Eq. (38) gives

$$ \frac{dx_j}{dt} = f_j(t), \quad j = 1, 2. \qquad (46) $$

That is after one turn,

$$ x_j((n+1)T) = \mu_j \left[ x_j(nT) + \int_{nT}^{(n+1)T} dt f_j(t) \right], \qquad (47) $$

where we have used that $e_j(T) = \mu_j e_j$. When $|\mu_j| < 1$, the above recurrence relation implies that $x_j(nT)$ has a magnitude comparable to the integral on the r.h.s. of Eq. (47). In the simple case when $f_j(t)$ is periodic of period $T$, the integral is constant and $x_j(nT)$ tends toward the finite value $x_j^*$,

$$ x_j^* = \frac{\mu_j}{1 - \mu_j} \int_0^T dt f_j(t). \qquad (48) $$

However for $\mu_1 = 1$, Eq. (47) is an arithmetic series and $x_1$ grows with time, if the integral term in Eq. (47) does not vanish. This would lead to a breakdown of the validity of the perturbation expansion. This growth of the tangent component to the limit cycle reflects a slow phase change of the oscillations that we have anticipated in the last term in Eq. (39). We can now choose the phase evolution so that the integral term in Eq. (47) exactly vanishes and the validity of the perturbative expansion is preserved, requiring that

$$ \int_t^{T+t} du f_1(u) = 0. \qquad (49) $$

This determines the E-I module phase drift dynamics to be

$$ \frac{d\phi}{dt} = \frac{1}{T} \int_t^{T+t} du \, \mathbf{g}_1(u) \cdot \mathbf{F}_1^{\text{forc}}(u), \qquad (50) $$

where we have computed the component $f_1(t)$ of $\mathbf{F}$ (Eq. (39)) with the help of the first vector $\mathbf{g}_1(t) = (g_{1,E}(t), g_{1,I}(t))$ of the bi-orthogonal basis $(\mathbf{g}_1(t), \mathbf{g}_2(t))$, such that $\mathbf{g}_i(t) \cdot \mathbf{e}_j(t) = \delta_{ij}$. Eq. (50) can be interpreted in a usual way as the product of the phase response curves of the E-I module multiplied by the perturbations due to forcing [19,23].
Diffusion of the oscillation phase in an E-I module with a finite number of neurons

The finite number of spiking neurons in a module gives rise to stochastic effects and produces a diffusive behavior for the oscillation phase of an E-I module in the oscillatory regime. This can be well-captured in the rate-model description by including supplementary stochastic terms, as described by Eqs. (5,6).

These stochastic terms can be interpreted as forcing the deterministic E-I module oscillations. The resultant stochastic phase dynamics can be described by Eq. (50) with $F^{forc}(t)$ given by

$$F^{forc}(t) = \left( \frac{1}{\tau_E[I_E(t)]} \left[ w_{EE} \sqrt{\frac{\Phi_E[I_E(t)]}{N_E}} \xi_E(t) - w_{EI} \sqrt{\frac{\Phi_I[I_I(t)]}{N_I}} \xi_I(t) \right] \frac{1}{\tau_I[I_I(t)]} w_{IE} \sqrt{\frac{\Phi_E[I_E(t)]}{N_E}} \xi_E(t) \right). \quad (51)$$

Substitution in Eq. (50) and averaging over the noise gives a diffusive evolution for the phase $\phi$ [19],

$$\frac{d\phi}{dt} = \eta(t), \quad \langle\eta(t)\eta(t')\rangle = D_N\delta(t-t'). \quad (52)$$

The finite number of neurons in both the excitatory and inhibitory populations contributes to the diffusion constant $D_N$, as described by Eq. (9). The two constants $D_E$ and $D_I$ in Eq. (9) are expressed as integrals along the limit cycle as

$$D_E = \frac{1}{T} \int_0^T dt \left( \frac{g_{1,E}(t) w_{EE}}{\tau_E[I_E^{(0)}(t)]} + \frac{g_{1,I}(t) w_{IE}}{\tau_I[I_I^{(0)}(t)]} \right)^2 \Phi_E[I_E^{(0)}(t)], \quad (53)$$

$$D_I = \frac{1}{T} \int_0^T dt \left( \frac{g_{1,E}(t) w_{EI}}{\tau_E[I_E^{(0)}(t)]} \right)^2 \Phi_I[I_I^{(0)}(t)]. \quad (54)$$

For the synaptic weights $w_{EE} = 1.6$ mVs, $w_{EI} = 0.32$ mVs, $w_{IE} = 2.0$ mVs, Eqs. (53,54) give $D_E = 1.2 \cdot 10^4$ ms, $D_I = 2.0 \cdot 10^3$ ms (Table 1).

Stochastic effect produces a diffusion of the oscillatory phase, the mean square amplitude of which grows linearly in time, as described by Eq. (8). This gives rise to a loss of oscillation coherence with time and to the decay of the activity correlation functions. Stochastic effects in the two neuronal populations also produce small changes in the oscillation frequency [19, 51, 52] of magnitudes $1/N_E$ and $1/N_I$, changes that we do not consider here.

The autocorrelation decay time $\tau_D$ can be directly related to the phase diffusion constant $D_N$ as follows. We consider, for definiteness, the activity of the excitatory population $r_E(t)$. For a very large number of neurons, noise is negligible and $r_E(t) = r_E^{(0)}(t)$ is a periodic function of period $T$. It can be written in a Fourier series as

$$r_E^{(0)}(t) = \sum_{n=-\infty}^{+\infty} \tilde{r}_n \exp \left( i \frac{2\pi n t}{T} \right). \quad (55)$$

For a finite but large number of neurons, the activity of the excitatory population can approximately be described by taking the dominant effect of phase diffusion into account, namely approximating $r_E(t)$ by $r_E^{(0)}(t + \phi(t))$. With the help of Eq. (55), it is then easy
to compute the auto-correlation of the excitatory activity $C_{EE}(t)$. That is,

$$ 
\begin{aligned}
C_{EE}(t) & \approx \langle r_E^{(0)}(s+\phi(s))r_E^{(0)}(s+t+\phi(s+t)) \rangle - \langle r_E^{(0)}(s) \rangle^2 \\
&= \sum_{n=-\infty}^{+\infty} \tilde{r}_n \tilde{r}_{-n} \exp \left( i \frac{2\pi nt}{T} \right) \exp \left( -\frac{2\pi^2 n^2}{T^2} \langle [\phi(s) - \phi(s+t)]^2 \rangle \right) - \tilde{r}_0^2 \\
&= 2 \sum_{n=1}^{+\infty} |\tilde{r}_n|^2 \cos \left( \frac{2\pi nt}{T} \right) \exp \left( -\frac{2\pi^2 n^2}{T^2} D_N t \right),
\end{aligned} 
\quad (56) $$

where the angular brackets indicate averages over different realizations of the noise, i.e., the fluctuating values of the phase $\phi(t)$. The auto-correlation decay time $\tau_D$ is thus obtained as given by Eq. (10). Eqs. (9,10) show that the oscillation coherence time $\tau_D$ grows linearly with the number of neurons in the E-I module. For our reference parameters (solid circle in Fig. 1A) where $T = 64$ ms, Eq. (10) gives $\tau_D \simeq 83$ ms for an E-I module with $N = 10^4$ neurons ($D_N = 2.5$ ms). For a module with ten times more neurons ($N = 10^5$), the decay time is ten times longer.

A fit of the activity auto-correlation decay with the above expression allows one to directly measure $\tau_D$ in simulations. This procedure is followed in Fig. 1F to compare the theoretical prediction to the simulation results.

## Synchronization function for two weakly coupled oscillatory E-I modules

The phase dynamics of a single forced oscillatory module (Eq. (50)) allows us to derive Eq. (13) governing the evolution of the respective phases of two weakly coupled oscillatory modules [19, 21, 24, 25].

We suppose that the parameters are such that for a vanishing coupling $f_{lr}$ (Eqs. (11,14)), the two E-I modules follow the limit cycle $\mathbf{I}^{(0)}(t) = (I_E^{(0)}(t), I_I^{(0)}(t))$ of period $T$ with different phases $\phi_1$ and $\phi_2$,

$$ \mathbf{I}_j(t) = \mathbf{I}^{(0)}(t + \phi_j), \quad j = 1, 2. \quad (57) $$

In order to determine how a small coupling $f_{lr}$ makes the phases $\phi_1$ and $\phi_2$ evolve, we compute perturbatively in $f_{lr}$ the corrections $\delta\mathbf{I}_j(t + \phi_j)$, $j = 1, 2$, to the zeroth-order evolution (57). We obtain that $\delta\mathbf{I}_1$ follows Eq. (38) with the following forcing current (Eq. (37)) that arises from the coupling to the second E-I module,

$$ \mathbf{F}^{\text{forc}}(t + \phi_1) = f_{lr} \mathbf{K}[\mathbf{I}^{(0)}(t + \phi_1), \mathbf{I}^{(0)}(t + \phi_2)]. \quad (58) $$

The perturbation comprises both long-range excitation coming from the other module and missing local excitation within the module, which is instead used as long-range excitation onto the other module.

More precisely, for $E \to E$ connectivity (Eqs. (11,12)), the forcing current is given by

$$ \mathbf{K}[\mathbf{I}^{(0)}(t), \mathbf{I}^{(0)}(t + \phi_2 - \phi_1)] = \left( \begin{array}{c} \frac{w_{EE}}{\tau_E[I_E^{(0)}(t)]} \left( \Phi_E[I_E^{(0)}(t + \phi_2 - \phi_1)] - \Phi_E[I_E^{(0)}(t)] \right) \\ 0 \end{array} \right). \quad (59) $$

Eq. (50) then determines the phase drift of the first E-I module to be

$$ 
\begin{aligned}
\frac{d\phi_1}{dt} &= \frac{f_{lr}}{T} \int_0^T dt \, \mathbf{g}_1(t) \cdot \mathbf{K}[\mathbf{I}^{(0)}(t), \mathbf{I}^{(0)}(t + \phi_2 - \phi_1)] \\
&= \frac{f_{lr} w_{EE}}{T} \int_0^T dt \frac{g_{1,E}(t)}{\tau_E[I_E^{(0)}(t)]} \left( \Phi_E[I_E^{(0)}(t + \phi_2 - \phi_1)] - \Phi_E[I_E^{(0)}(t)] \right).
\end{aligned} 
\quad (60) $$
The evolution of $\phi_2$ is obtained by permutation of $\phi_1$ and $\phi_2$ in Eq. (60). Subtraction of these two equations provides the evolution of the relative phase $\Delta\phi = \phi_1 - \phi_2$ of the two modules as given by Eq. (13), with the explicit expression of the coupling function $S_E$ being

$$S_E(\Delta\phi) = \frac{w_{EE}}{T} \int_0^T dt \frac{g_{1,E}(t)}{\tau_E[I_E^{(0)}(t)]} \left( \Phi_E[I_E^{(0)}(t - \Delta\phi)] - \Phi_E[I_E^{(0)}(t + \Delta\phi)] \right). \quad (61)$$

The subscript E of $S_E$ denotes that the expression holds in the case when long-range excitation only targets excitatory neurons.

For a small phase difference, one readily obtains to first order in $\Delta\phi$

$$S_E(\Delta\phi) \approx -2\Delta\phi D_{\phi}^{E}, \quad (62)$$

with

$$D_{\phi}^{E} = \frac{w_{EE}}{T} \int_{0}^{T} dt \frac{g_{1,E}(t)}{\tau_{E}[I_{E}^{(0)}(t)]} \Phi'_{E}[I_{E}^{(0)}(t)] \frac{d}{dt} I_{E}^{(0)}(t). \quad (63)$$

The constant $D_{\phi}^E$ is actually negative (see Table 1 for our reference parameters) so that full oscillatory synchrony is unstable for small long-range excitatory coupling $f_{flr}$ between the two modules.

For $E \to E, I$ connectivity (Eqs. (11,14)), the coupling function $K$ of Eq. (59) is replaced by

$$K[I^{(0)}(t), I^{(0)}(t + \phi_2 - \phi_1)] = \left( \begin{array}{c} \frac{w_{EE}}{\tau_E[I_E^{(0)}(t)]} \left( \Phi_E[I_E^{(0)}(t + \phi_2 - \phi_1)] - \Phi_E[I_E^{(0)}(t)] \right) \\ \frac{w_{IE}}{\tau_I[I_I^{(0)}(t)]} \left( \Phi_E[I_E^{(0)}(t + \phi_2 - \phi_1)] - \Phi_E[I_E^{(0)}(t)] \right) \end{array} \right) \quad (64)$$

This gives in turn the coupling function for the phase difference $\Delta\phi$ between the two modules,

$$S_{EI}(\Delta\phi) = \int_{0}^{T} \frac{dt}{T} \left( \frac{g_{1,E}(t) w_{EE}}{\tau_E[I_E^{(0)}(t)]} + \frac{g_{1,I}(t) w_{IE}}{\tau_I[I_I^{(0)}(t)]} \right) \left( \Phi_E[I_E^{(0)}(t - \Delta\phi)] - \Phi_E[I_E^{(0)}(t + \Delta\phi)] \right). \quad (65)$$

For a small phase difference between the two modules, analogously to Eq. (61), Eq. (65) reduces to Eq. (15), the analog of Eq. (62) for $E \to E, I$ connectivity, with

$$D_{\phi}^{EI} = \frac{1}{T} \int_{0}^{T} dt \left( \frac{g_{1,E}(t) w_{EE}}{\tau_{E}[I_{E}^{(0)}(t)]} + \frac{g_{1,I}(t) w_{IE}}{\tau_{I}[I_{I}^{(0)}(t)]} \right) \Phi'_{E}[I_{E}^{(0)}(t)] \frac{d}{dt} I_{E}^{(0)}(t). \quad (66)$$

In this case, the constant $D_{\phi}^{EI}$ is positive (see Table 1 for our reference parameters) and full oscillatory synchrony is stable for small long-range excitatory coupling between the two modules.

## Stability analysis of full synchronization for two coupled E-I modules

The fully synchronized oscillatory state is always a possible dynamical state for two coupled identical E-I modules. Namely, $I_1(t) = I_2(t) = I^{(0)}(t)$ is an exact solution of the coupled Eqs. (17) and (18/19) since the missing local excitation is supplied by the distant modules. The stability of this oscillatory state can be determined irrespective of the strength of coupling between the two modules.

We consider slightly perturbed evolutions for the two modules $I_j(t) = I^{(0)}(t) + \delta I_j(t)$, $j = 1, 2$. The dynamical behavior of the perturbations $\delta I_j(t)$ is obtained by linearizing the dynamics around the fully synchronized oscillatory state,

$$\frac{d\delta\mathbf{I}_1}{dt} = \mathcal{L}_{11}[\mathbf{I}^{(0)}(t)]\delta\mathbf{I}_1 + \mathcal{L}_{12}[\mathbf{I}^{(0)}(t)]\delta\mathbf{I}_2, \quad (67)$$
together with the analogous equation with permuted indices 1 and 2.

Since the dynamics is invariant under the interchange of module 1 and 2, the 4-dimensional linear evolution of $(\delta\mathbf{I}_1(t), \delta\mathbf{I}_2(t))$ can be reduced to the study of the 2-dimensional evolution of symmetric $(\delta\mathbf{I}_1 = \delta\mathbf{I}_2)$ and antisymmetric $(\delta\mathbf{I}_1 = -\delta\mathbf{I}_2)$ perturbations. The evolution of symmetric perturbations is governed by the $2 \times 2$ matrix $\mathcal{L}_S(t)$ which is identical to the matrix $\mathcal{L}(t)$ (Eq. (41)) for a single E-I module,

$$ \mathcal{L}_S(t) = \mathcal{L}_{11} + \mathcal{L}_{12} = \mathcal{L}. \qquad (68) $$

The coupled system of two modules thus inherits its stability with respect to symmetric perturbations from the stability of the limit cycle of a single E-I module. On the contrary, antisymmetric perturbations evolve according to the new $2 \times 2$ matrix $\mathcal{L}_A(t)$,

$$ \mathcal{L}_A[\mathbf{I}^{(0)}(t)] = \mathcal{L}_{11} - \mathcal{L}_{12}. \qquad (69) $$

The matrix $\mathcal{L}_A(t)$ explicitely depends on the coupling $f_{lr}$. The stability of the synchronized state is determined by assessing whether antisymmetric perturbations have grown or decayed after one period $T$. As above, this is governed by the eigenvalues of the monodromy matrix $M_A(T)$ associated with $\mathcal{L}_A(t)$,

$$ \frac{d}{dt} M_A = \mathcal{L}_A[\mathbf{I}^{(0)}(t)] \cdot M, \quad M_A(0) = \mathrm{Id}. \qquad (70) $$

For $E \to E$ connectivity (Eqs. (11,12)), the matrices $\mathcal{L}_{11}$ and $\mathcal{L}_{12}$ read

$$ \mathcal{L}_{11}^{E}[\mathbf{I}^{(0)}(t)] = \begin{pmatrix} \frac{-1+w_{EE}(1-f_{lr})\Phi'_{E}[I_{E}^{(0)}(t)]-\tau'_{E}[I_{E}^{(0)}(t)]I'_{E}^{(0)}(t)}{\tau_{E}[I_{E}^{(0)}(t)]} & -\frac{w_{EI}\Phi'_{I}[I_{I}^{(0)}(t)]}{\tau_{E}[I_{E}^{(0)}(t)]} \\ \frac{w_{IE}\Phi'_{E}[I_{E}^{(0)}(t)]}{\tau_{I}[I_{I}^{(0)}(t)]} & -\frac{1+\tau'_{I}[I_{I}^{(0)}(t)]I'_{I}^{(0)}(t)}{\tau_{I}[I_{I}^{(0)}(t)]} \end{pmatrix}, \qquad (71) $$

$$ \mathcal{L}_{12}^{E}[\mathbf{I}^{(0)}(t)] = \begin{pmatrix} \frac{f_{lr} w_{EE} \Phi'_{E} [I_{E}^{(0)}(t)]}{\tau_{E}[I_{E}^{(0)}(t)]} & 0 \\ 0 & 0 \end{pmatrix}. \qquad (72) $$

Integration of Eq. (70) for different $f_{lr}$ and other parameters fixed, shows that the two eigenvalues of $M_A^E(T)$ are smaller than 1 for $f_{lr}$ above a threshold coupling $f_{lr}^*$, while one eigenvalue is larger than 1 for $f_{lr} < f_{lr}^*$. The obtained $f_{lr}^*$ is shown as a solid red line in Fig. 2B.

For $E \to E, I$ connectivity (Eqs. (11,14)), the calculation is identical but the resulting matrix $\mathcal{L}_A^{EI}(t)$ reads

$$ \mathcal{L}_A^{EI}[\mathbf{I}^{(0)}(t)] = \begin{pmatrix} \frac{-1+w_{EE}(1-2f_{lr})\Phi'_{E}[I_{E}^{(0)}(t)]-\tau'_{E}[I_{E}^{(0)}(t)]I'_{E}^{(0)}(t)}{\tau_{E}[I_{E}^{(0)}(t)]} & -\frac{w_{EI}\Phi'_{I}[I_{I}^{(0)}(t)]}{\tau_{E}[I_{E}^{(0)}(t)]} \\ \frac{(1-2f_{lr})w_{IE}\Phi'_{E}[I_{E}^{(0)}(t)]}{\tau_{I}[I_{I}^{(0)}(t)]} & -\frac{1+\tau'_{I}[I_{I}^{(0)}(t)]I'_{I}^{(0)}(t)}{\tau_{I}[I_{I}^{(0)}(t)]} \end{pmatrix}. \qquad (73) $$

The eigenvalues of the associated monodromy matrix $M_A^{EI}(T)$ are found to be smaller than 1 irrespective of the value of $f_{lr}$ (see Fig. 6G), showing that the fully synchronized state is stable in this case.

In the previous section, it was shown that the stability of the synchronized state for small $f_{lr}$ was dependent on the sign of the constants $D_{\phi}^E$ (Eq. (63)) and $D_{\phi}^{EI}$ (Eq. (66)), respectively. This can also be seen directly by computing the largest eigenvalue $\mu_1^A$ of $M_A(T)$ (Eq. (70)) for small $f_{lr}$ as we shall demonstrate below.
When $f_{lr} = 0$, the two modules are not coupled and $M_A(T)$ reduces to $M(T)$, the matrix for a single module (Eq. (42)), with $\mu_1^A = \mu_1 = 1$. For small $f_{lr} \to 0$, $\mu_1^A$ can be computed perturbatively. We write the expansion of $\mathcal{L}_A(t)$ and $M_A(T)$ in powers of $f_{lr}$ to first order as

$$ \mathcal{L}_A(t) = \mathcal{L}(t) + f_{lr} \delta\mathcal{L}_A(t), \quad M_A(T) = M(T) + f_{lr} \delta M_A, \qquad (74) $$

where $\mathcal{L}(t)$ is given by Eq. (41) and the matrix $M$ by Eq. (42). The matrix $\delta M_A$ is obtained from the integration of Eq. (70) to linear order in $f_{lr}$,

$$ \delta M_A = M(T) \int_0^T du M^{-1}(u) \delta \mathcal{L}_A(u) M(u). \qquad (75) $$

First-order perturbation theory gives the development of the eigenvalue $\mu_1^A$ as a function of $\delta M_A$,

$$ \mu_1^A = 1 + f_{lr} \mathbf{g}_1 \cdot \delta M_A \mathbf{e}_1, \qquad (76) $$

where $\mathbf{g}_1$ the left eigenvector of $M(T)$ associated to the eigenvalue 1, such that $\mathbf{g}_1 \cdot \mathbf{e}_1 = 1$ (i.e., the first vector of the bi-orthogonal basis $(\mathbf{g}_1, \mathbf{g}_2)$). Explicit formulas are obtained after the replacement of $\delta M_A$ by its expression (75) together with the identity $\mathbf{g}_1(u) = \mathbf{g}_1 M^{-1}(u)$ and the expressions of $\delta \mathcal{L}_A$ for the two considered connectivities. For $E \to E, I$ connectivity, Eq. (73) gives $\delta \mathcal{L}_A^{EI}$, and subsequently the expression for the largest eigenvalue of $M_A^{EI}(T)$,

$$ \mu_1^A - 1 = -2f_{lr} D_{\phi}^{EI} T + \mathcal{O}(f_{lr}^2). \qquad (77) $$

The same formula with $D_{\phi}^{EI}$ replaced by $D_{\phi}^E$ is obtained for the largest eigenvalue of $M_A^{EI}(T)$ with the help of Eqs. (71) and (72). Since $D_{\phi}^{EI} > 0$ and $D_{\phi}^E < 0$, one recovers that $E \to E, I$ connectivity synchronizes the two modules' oscillations ($\mu_1^A < 1$) while full synchronization is unstable ($\mu_1^A > 1$) at small coupling for $E \to E$ connectivity.

## Competition between synchronization and stochasticity for two coupled E-I modules of finite size

For $E \to E, I$ connectivity, long-range excitation between two E-I modules tends to synchronize their oscillations. In contrast, finite-size stochasticity tends to make their oscillatory phases drift independently. The competition between these two effects can be analytically estimated when noise is weak. We first consider the simplest case when the coupling between the two modules is also weak, before considering the case of an arbitrary coupling.

### Weakly-coupled modules

The phase dynamics of weakly-coupled modules is simply described by the linear addition of the stochastic diffusion of the two modules' phases (Eqs. (8,52)) and synchronization due to long-range coupling (Eq. (13)). Thus, one obtains for the phase difference dynamics

$$ \frac{d(\phi_1 - \phi_2)}{dt} = f_{lr} S_{EI}(\phi_1 - \phi_2) + \eta_{12}(t) \approx -2f_{lr} D_{\phi}^{EI}(\phi_1 - \phi_2) + \eta_{12}(t), \qquad (78) $$

where in the second approximate equality we have replaced $S_{EI}(\Delta\phi)$ by its linear approximation for small phase differences (Eq. (15)). The noise term $\eta_{12}(t)$ is the difference of the independent finite-size noises for the two modules (52). Being a linear combination of white noises, $\eta_{12}(t)$ is also a white noise, with an amplitude twice as large as for a single module,

$$ \eta_{12}(t) = \eta_1(t) - \eta_2(t), \quad \langle \eta_{12}(t) \eta_{12}(t') \rangle = 2D_N \delta(t - t'). \qquad (79) $$
Eqs. (78) and (79) readily give that for a small noise amplitude ($D_N \ll f_{lr} D_{\phi}^{EI} T$), the two-module oscillation phase difference has Gaussian fluctuations with a mean square amplitude given by Eq. (16).

We can use that result to estimate the mean square amplitude of the excitatory firing rate difference $\delta r_E = (r_{1,E} - r_{2,E})/\sqrt{2}$ between the two modules. In the case of weak coupling, the firing rates can be approximated as $r_{j,E}(t) \approx r_E^{(0)}(t) + r_E^{(0)}(t) \phi_j(t)$, $j = 1, 2$, as long as $|\phi_j| \ll T$. The firing rate fluctuations are then directly obtained from the phase differences as

$$ \begin{aligned} \langle \delta r_E(t)^2 \rangle_T &= \frac{1}{2} \langle (\phi_1 - \phi_2)^2 \rangle \int_0^T \frac{dt}{T} \left[ \Phi'_E [I_E^{(0)}(t)] \frac{dI_E^{(0)}}{dt}(t) \right]^2 \\ &= \frac{D_N}{4 f_{lr} D_{\phi}^{EI}} \int_0^T \frac{dt}{T} \left[ \Phi'_E [I_E^{(0)}(t)] \frac{dI_E^{(0)}}{dt}(t) \right]^2, \end{aligned} \quad (80) $$

where we used $r_E^{(0)}(t) = \Phi'_E [I_E^{(0)}(t)] I_E^{(0)}(t)$ and furthermore averaged over one oscillation period, indicated by the subscript $T$. Note that we neglected here the additional direct, i.e. not phase-mediated, contribution due to Poissonian sampling, which would add an additional term $r_E^{(0)}(t)/(N_E \Delta t)$ with $\Delta t$ being the spike sampling interval.

The above expression diverges with vanishing $f_{lr}$. This arises from the fact that the phase difference between the two oscillating modules becomes arbitrarily large when the inter-module coupling vanishes. However a large phase difference does not imply large rate differences since for completely uncorrelated oscillating modules the rate fluctuations are given by $\langle r_E^2 \rangle_T - \langle r_E \rangle_T^2$. The unphysical divergence of expression (80) for small $f_{lr}$ comes from the neglect of the periodicity of the rate as function of the oscillation phase. We can obtain an expression for $\langle \delta r_E(t)^2 \rangle_T$ that takes this periodicity into account. Using the Fourier decomposition Eq. (55) and $r_{j,E}(t) \approx r_E^{(0)}(t + \phi_j(t))$, $j = 1, 2$, one obtains analogously to the single module autocorrelation

$$ \begin{aligned} \langle r_{1,E}(t)r_{2,E}(t) \rangle_T &\approx \int_0^T \frac{dt}{T} \sum_{n,m} \tilde{r}_n \tilde{r}_m \exp \left( i \frac{2\pi(n+m)t}{T} \right) \left\langle \exp \left( i \frac{2\pi[n\phi_1(t) + m\phi_2(t)]}{T} \right) \right\rangle \\ &= \sum_n |\tilde{r}_n|^2 \left\langle \exp \left( i \frac{2\pi n}{T} [\phi_1(t) - \phi_2(t)] \right) \right\rangle \\ &= \sum_n |\tilde{r}_n|^2 \exp \left( - \frac{2\pi^2 n^2}{T} \langle [\phi_1(t) - \phi_2(t)]^2 \rangle \right). \end{aligned} \quad (81) $$

Here, we first used that $\frac{1}{T} \int_0^T dt e^{i2\pi(n+m)t/T} = \delta_{n,-m}$ and that the relative phases $\phi_1(t), \phi_2(t)$ vary slowly with time with respect to the oscillatory dynamics, i.e., can be considered constant over one oscillation period.

Using $\langle \delta r_E^2 \rangle_T = (\langle r_{1,E}^2 \rangle_T + \langle r_{2,E}^2 \rangle_T)/2 - \langle r_{1,E} r_{2,E} \rangle_T$, it is now straightforward to obtain

$$ \langle \delta r_E(t)^2 \rangle_T = 2 \sum_{n=1}^{+\infty} |\tilde{r}_n|^2 \left[ 1 - \exp \left( - \frac{2\pi^2 n^2}{T} \langle (\phi_1 - \phi_2)^2 \rangle \right) \right]. \quad (82) $$

A comparison of this theoretical result with simulations of the stochastic FAT rate model and spiking networks is shown in Fig. 5.

### Weak noise at arbitrary coupling

Eq. (16) provides a simple estimate of the competition between stochastic dephasing and synchronization for weak coupling between the two modules. A more precise estimate
for small noise but arbitrary coupling is obtained by considering general fluctuations around the fully synchronous state that result from the modules' stochastic activities. Namely, we describe the dynamics around the fully synchronized state, at the linear level, as in Eq. (67), $\mathbf{I}_j(t) = \mathbf{I}^{(0)}(t) + \delta\mathbf{I}_j(t)$, $j = 1, 2$, and compute the dynamics of the perturbations $\delta\mathbf{I}_j(t)$ resulting from the presence of stochasticity. For the first module, this gives

$$ \frac{d\delta\mathbf{I}_1}{dt} = \mathcal{L}_{11}[\mathbf{I}^{(0)}(t)]\delta\mathbf{I}_1 + \mathcal{L}_{12}[\mathbf{I}^{(0)}(t)]\delta\mathbf{I}_2 + \mathbf{B}_1 + f_{lr}\mathbf{B}_{12}, \quad (83) $$

with the stochastic terms

$$ \mathbf{B}_1 = \left( \frac{\frac{w_{EE}}{\tau_E[I_E^{(0)}(t)]}}{\frac{w_{IE}}{\tau_I[I_I^{(0)}(t)]}} \right) (1 - f_{lr}) \sqrt{\frac{\Phi_E[I_E^{(0)}(t)]}{N_E}} \xi_{E,1}(t) - \left( \frac{\frac{w_{EI}}{\tau_E[I_E^{(0)}(t)]}}{0} \right) \sqrt{\frac{\Phi_I[I_I^{(0)}(t)]}{N_I}} \xi_{I,1}(t), $$

$$ \mathbf{B}_{12} = \left( \frac{\frac{w_{EE}}{\tau_E[I_E^{(0)}(t)]}}{\frac{w_{IE}}{\tau_I[I_I^{(0)}(t)]}} \right) \sqrt{\frac{\Phi_E[I_E^{(0)}(t)]}{N_E}} \xi_{E,2}(t) \quad (84) $$

Permuting indices 1 and 2 in Eqs. (83) and (84) gives the dynamics of the linear departure $\delta\mathbf{I}_2$ of the second module from the fully synchronized state.

As previously, the permutation symmetry between module 1 and 2 can be taken advantage of to reduce the 4-dimensional dynamics to a pair of uncoupled two-dimensional dynamics for the symmetric and antisymmetric modes. We focus on the antisymmetric mode $\delta\mathbf{I}_A(t)$ in the following,

$$ \delta\mathbf{I}_A = \frac{1}{\sqrt{2}}(\delta\mathbf{I}_1 - \delta\mathbf{I}_2). \quad (85) $$

From Eq. (83) and the analogous equation for module 2, the antisymmetric mode obeys

$$ \frac{d\delta\mathbf{I}_A}{dt} = \mathcal{L}_A^{EI}[\mathbf{I}^{(0)}(t)]\delta\mathbf{I}_A + \mathbf{B}_A^{EI}, \quad (86) $$

with the matrix $\mathcal{L}_A^{EI}$ for $E \to E, I$ connectivity given by Eq. (73). The stochastic term $\mathbf{B}_A^{EI}$ is given by

$$ \mathbf{B}_A^{EI} = \left( \frac{\frac{w_{EE}}{\tau_E[I_E^{(0)}(t)]}}{\frac{w_{IE}}{\tau_I[I_I^{(0)}(t)]}} \right) (1 - 2f_{lr}) \sqrt{\frac{\Phi_E[I_E^{(0)}(t)]}{N_E}} \xi_{E,A}(t) - \left( \frac{\frac{w_{EI}}{\tau_E[I_E^{(0)}(t)]}}{0} \right) \sqrt{\frac{\Phi_I[I_I^{(0)}(t)]}{N_I}} \xi_{I,A}(t) \quad (87) $$

with

$$ \xi_{E,A}(t) = \frac{1}{\sqrt{2}}[\xi_{E,1}(t) - \xi_{E,2}(t)], \quad \xi_{I,A}(t) = \frac{1}{\sqrt{2}}[\xi_{I,1}(t) - \xi_{I,2}(t)]. \quad (88) $$

In order to solve Eq. (86), we introduce again the matrix $M_A^{EI}(t)$, as defined by Eq. (70). The matrix $M_A^{EI}(T)$ has eigenvectors $\mathbf{e}_1^A, \mathbf{e}_2^A$ with eigenvalues $\mu_1^A$ and $\mu_2^A$. The two eigenvalues and eigenvectors are real when the coupling $f_{lr}$ between the two modules is sufficiently small ($f_{lr} \lesssim 0.13$ for our reference parameters, see Fig. 6G), but they can be complex for larger couplings.

Eq. (86) can be explicitely integrated by writing $\delta\mathbf{I}_A(t)$ and $\mathbf{B}_A(t)$ in the basis obtained by translating $\mathbf{e}_1^A, \mathbf{e}_2^A$ along the limit cycle,

$$ \mathbf{e}_j^A(t) = M_A(t) \mathbf{e}_j^A, \quad j = 1, 2. \quad (89) $$

Analogously to a single oscillatory E-I module subject to a forcing current (Eqs. (43,44)), we define

$$ \delta\mathbf{I}_A(t) = x_1(t) \mathbf{e}_1^A(t) + x_2(t) \mathbf{e}_2^A(t), \quad (90) $$

$$ \mathbf{B}_A(t) = b_1(t) \mathbf{e}_1^A(t) + b_2(t) \mathbf{e}_2^A(t). \quad (91) $$
Substitution of these expressions in the differential equation (86) simply gives

$$ \frac{dx_j}{dt} = b_j(t), \quad j = 1, 2. \qquad (92) $$

One readily obtains, for $0 \le t \le T$,

$$ x_j(t) = \int_{-\infty}^{t} du b_j(u) = \int_{0}^{t} du \mathbf{g}_j^A(u) \cdot \mathbf{B}_A(u) + \sum_{n=1}^{+\infty} \mu_j^n \int_{0}^{t} du \mathbf{g}_j^A(u) \cdot \mathbf{B}_A(u - nT). \quad (93) $$

We note that

$$ \langle \mathbf{g}_j^A(u) \cdot \mathbf{B}_A(u) \mathbf{g}_k^A(v) \cdot \mathbf{B}_A(v) \rangle = b_{jk} \delta(u - v), \qquad (94) $$

with

$$ b_{jk} = (1 - 2f_{lr})^2 \left( \frac{w_{EE} g_{j,E}^A(u)}{\tau_E[I_E^{(0)}(u)]} + \frac{w_{IE} g_{j,I}^A(u)}{\tau_I[I_I^{(0)}(u)]} \right) \left( \frac{w_{EE} g_{k,E}^A(u)}{\tau_E[I_E^{(0)}(u)]} + \frac{w_{IE} g_{k,I}^A(u)}{\tau_I[I_I^{(0)}(u)]} \right) \frac{\Phi_E[I_E^{(0)}(u)]}{N_E} \\ + \frac{w_{EI}^2 g_{j,E}^A(u) g_{k,E}^A(u)}{\tau_E^2[I_E^{(0)}(u)]} \frac{\Phi_I[I_I^{(0)}(t)]}{N_I} \quad (95) $$

The amplitude of the current fluctuations along the limit cycle is finally obtained, with the help of Eqs. (90), (93), and (94), as

$$ \langle \delta I_{A,E}(t)^2 \rangle = \sum_{j,k=1,2} e_{j,E}^A(t) e_{k,E}^A(t) \left[ \int_0^t du b_{jk}(u) + \frac{\mu_j \mu_k}{1 - \mu_j \mu_k} \int_0^T du b_{jk}(u) \right]. \quad (96) $$

Since $e_j^A(T) = \mu_j^A e_j^A$, $j = 1, 2$, the obtained expression is periodic in time, as it should.

Because we are considering fluctuations around the stationary limit cycle to linear order, the firing rates can be approximated by $r_{j,E}(t) = \Phi'_E[I_E^{(0)}(t)]\delta I_{j,E}(t)$, $j = 1, 2$, up to the stochastic contribution due to Poisson sampling. The amplitudes of the rate fluctuations $\langle \delta r_E(t)^2 \rangle = \langle [r_{1,E}(t) - r_{2,E}(t)]^2 \rangle/2$ along the limit cycle, up to the direct Poisson contribution, are thus directly obtained by multiplying the above expression with $\phi'_E[I_E^{(0)}(t)]^2$. Averaging over one oscillation period, we find

$$ \langle \delta r_E(t)^2 \rangle_T = \frac{1}{T} \int_0^T dt \sum_{j,k=1,2} \phi'_E[I_E^{(0)}(t)]^2 e_{j,E}^A(t) e_{k,E}^A(t) \\ \times \left[ \int_0^t du b_{jk}(u) + \frac{\mu_j \mu_k}{1 - \mu_j \mu_k} \int_0^T du b_{jk}(u) \right]. \quad (97) $$

For larger coupling strengths, this result is found to be much more precise when compared to numerical simulations than the expressions obtained in the weak-coupling limit, Eq. (80) or respectively (82) when taking the periodicity of the limit cycle dynamics into account, see Fig. 5. Note however that expression (97) is valid only for small perturbations and fails to account for the limited rate correlations at vanishing coupling due to the periodicity of the limit cycle dynamics, as done in Eq. (82).

In the weak-coupling regime ($f_{lr} \to 0$), Eq. (97) gives back the simple expression (16) obtained from the stochastic phase equation. This can be seen as follows. When $f_{lr} \to 0$, $\mu_1^A$ tends towards 1 (Eq. (77)), its value for a single module. Therefore, the r.h.s. of Eq. (97) is dominated by the $j = k = 1$ contribution in the sum and more precisely by its second term, the denominator of which vanishes. Namely, with the help of Eq. (77),
when $f_{lr} \to 0$ one obtains

$$ \begin{aligned} \langle \delta r_E(t)^2 \rangle_T & \approx \frac{1}{T} \int_0^T dt \, \phi'_E [I_E^{(0)}(t)]^2 \frac{(e_{1,E}^A(t))^2}{4 f_{lr} D_\phi^{EI} T} \int_0^T du \, b_{11}(u) \\ & \approx \frac{1}{T} \int_0^T dt \, \phi'_E [I_E^{(0)}(t)]^2 \frac{(e_{1,E}^A(t))^2}{4 f_{lr} D_\phi^{EI} T} \left( \frac{D_I}{N_I} + (1 - 2 f_{lr})^2 \frac{D_E}{N_E} \right), \end{aligned} \quad (98) $$

where the last equality is obtained by comparing the definition (95) to the previous expressions for $D_E$ and $D_I$ and noting that the vectors $g_j^A$ and $e_1^A$ tend toward those of the single module limit cycle, $g_j$ and $e_1$, when $f_{lr} \to 0$. Finally, in the limit of weak coupling and weak noise where the phase is well-defined, the mean square difference of excitatory firing rates between the two modules is related to their mean square phase difference,

$$ \langle \delta r_E(t)^2 \rangle_T = \frac{1}{2} \langle (\phi_1 - \phi_2)^2 \rangle \int_0^T \frac{dt}{T} \left[ \phi'_E [I_E^{(0)}(t)] \frac{dI_E^{(0)}}{dt}(t) \right]^2 . \quad (99) $$

(We used this relation before when deriving Eq. (80).) Comparison of Eq. (98) and Eq. (99) gives back the expression (16) for the mean square phase difference when one remembers that the tangent vector $e_1(t)$ is simply the velocity along the current limit cycle.

## Stability analysis of full synchronization for a chain of oscillatory E-I modules

The fully synchronized oscillatory state $\mathbf{I}_n(t) = \mathbf{I}^{(0)}(t)$, $n = 1, \dots, L$, is also always an exact solution for the dynamics of a chain of $L$ E-I modules (Eqs. (17) and (18) or (19) depending on the connectivity). Its stability can be assessed very similarly to the previous case of two coupled modules.

We consider slightly perturbed evolutions for the $L$ modules $\mathbf{I}_n(t) = \mathbf{I}^{(0)}(t) + \delta \mathbf{I}_n(t)$, $n = 1, \dots, L$. The linear evolution of the perturbations $\delta \mathbf{I}_n$ is found to be

$$ \frac{d\delta\mathbf{I}_n}{dt} = \sum_{m=1,\dots,L} \mathcal{L}_{nm}[\mathbf{I}^{(0)}(t)]\delta\mathbf{I}_m. \quad (100) $$

Since the chain is invariant by translation, one can reduce this $2L$-dimensional evolution to the evolutions of $L$ independent 2-dimensional systems. Namely, we write the perturbations as

$$ \delta \mathbf{I}_n(t) = \exp(iqn) \delta \tilde{\mathbf{I}}_q, \quad n = 1, \dots, L, \quad (101) $$

with $L$ “wavevectors” $q = 2\pi k/L$, $n = 0, \dots, L-1$. This gives, for each $q$, the evolution

$$ \frac{d\delta\tilde{\mathbf{I}}_q}{dt} = \mathcal{L}_q[\mathbf{I}^{(0)}(t)]\delta\tilde{\mathbf{I}}_q \quad (102) $$

with $2 \times 2$ matrices $\mathcal{L}_q$ that we give below.

For $E \to E$ connectivity, the matrices $\mathcal{L}_{nm}$ read, for $n = m$,

$$ \mathcal{L}_{nn}^{E}[\mathbf{I}^{(0)}(t)] = \begin{pmatrix} \frac{-1+w_{EE}C(0)\Phi'_{E}[I_{E}^{(0)}(t)]-\tau'_{E}[I_{E}^{(0)}(t)]I'_{E}^{(0)}(t)}{\tau_{E}[I_{E}^{(0)}(t)]} & -\frac{w_{EI}\Phi'_{I}[I_{I}^{(0)}(t)]}{\tau_{E}[I_{E}^{(0)}(t)]} \\ \frac{w_{IE}\Phi'_{E}[I_{E}^{(0)}(t)]}{\tau_{I}[I_{I}^{(0)}(t)]} & -\frac{1+\tau'_{I}[I_{I}^{(0)}(t)]I'_{I}^{(0)}(t)}{\tau_{I}[I_{I}^{(0)}(t)]} \end{pmatrix}, \quad (103) $$
and for $n \neq m$,

$$ \mathcal{L}_{nm}^{E}[\mathbf{I}^{(0)}(t)] = \begin{pmatrix} \frac{w_{EE}C(n-m)\Phi'_{E}[I_{E}^{(0)}(t)]}{\tau_{E}[I_{E}^{(0)}(t)]} & 0 \\ 0 & 0 \end{pmatrix}. \quad (104) $$

The corresponding matrices $\mathcal{L}_q^E$ read

$$ \mathcal{L}_q^E[\mathbf{I}^{(0)}(t)] = \begin{pmatrix} \frac{-1+w_{EE}\tilde{C}_\lambda(q)\Phi'_E[I_E^{(0)}(t)]-\tau'_E[I_E^{(0)}(t)]I'_E^{(0)}(t)}{\tau_E[I_E^{(0)}(t)]} & -\frac{w_{EI}\Phi'_I[I_I^{(0)}(t)]}{\tau_E[I_E^{(0)}(t)]} \\ \frac{w_{IE}\Phi'_E[I_E^{(0)}(t)]}{\tau_I[I_I^{(0)}(t)]} & -\frac{1+\tau'_I[I_I^{(0)}(t)]I'_I^{(0)}(t)}{\tau_I[I_I^{(0)}(t)]} \end{pmatrix} \quad (105) $$

The matrix $\mathcal{L}_q^E$ is identical to the previous matrix $\mathcal{L}_A^E$ (Eqs. (69,71,72)) governing the stability of two coupled modules with the replacement of $(1 - 2f_{lr})$ by the effective coupling $\tilde{C}_\lambda(q)$,

$$ \tilde{C}_{\lambda}(q) = \sum_{l=0, \dots, L-1} \exp(iql)C(l) = \frac{\cosh(\lambda) - 1}{\cosh(\lambda) - \cos(q)}, \quad (106) $$

where the second equality specifically holds for the exponentially decreasing excitatory interaction (20). The instability of full synchronization for two-coupled E-I modules for $f_{lr} < f_{lr}^*$ translates into an instability of full synchronization in a chain of modules for $\tilde{C}_\lambda(q)$ close enough to 1. Namely, spontaneous phase gradients appear in the chain of E-I modules for $|q| < q^*$ with the threshold $q^*$ given by Eq. (22) in the main text (for the coupling choice (20)).

For $E \to E, I$ connectivity, the matrix $\mathcal{L}_q^{EI}$ is similarly given by the previous matrix $\mathcal{L}_A^{EI}$ (Eq. (73)) with $(1 - 2f_{lr})$ replaced by $\tilde{C}_\lambda(q)$. Therefore in this case, the eigenvalues of $\mathcal{L}_q^{EI}$ are of magnitude smaller than 1, perturbations of any wavelength tend to vanish, and full synchronization of the chain modules is stable. The evolution of a long wavelength mode can be directly transcribed from our previous results in the two-module case (Eqs. (13) and (15) or Eq. (77)), using the correspondence between $f_{lr}$ and $\tilde{C}_\lambda(q)$,

$$ \frac{d\delta\tilde{\mathbf{I}}_q}{dt} = D_{\phi}^{EI} [\tilde{C}_{\lambda}(q) - 1] \delta\tilde{\mathbf{I}}_q(t) \approx -\frac{D_{\phi}^{EI} q^2}{2[\cosh(\lambda) - 1]} \delta\tilde{\mathbf{I}}_q(t), \quad (107) $$

where in the second approximate equality we used that $q \ll 1$. This provides the relaxation rate of long-wavelength modes given in Eq. (23) of the main text.

## Stochastic dynamics of the phases of oscillations along a chain of E-I modules

We analyze the competition between synchronization and stochastic fluctuation in a chain of E-I modules. We consider $E \to E, I$ connectivity (Eqs. (17,19)), for which the oscillations of all modules are stably synchronized in the absence of stochastic fluctuations of activity. We first consider the case when the variation of the phase of oscillation is small. The dynamics of these long wavelength modes can be fully described analytically and provide reference expressions. Shorter wavelength fluctuations are considered in the next section.

Taking Eq. (17) as an example, we rewrite the long-range synaptic coupling as

$$ \sum_{m=1, \dots, L} C(n, m) r_{E,m} = \Phi_E[I_{E,n}(t)] + \sum_{m=1, \dots, L} C(n, m) \left( \{\Phi_E[I_{E,m}(t)] - \Phi_E[I_{E,n}(t)]\} + \sqrt{\frac{\Phi_E[I_{E,m}(t)]}{N_E}} \xi_{E,m}(t) \right), \quad (108) $$
where we have made use of the expression (5) of the firing rate. The term between large parenthesis on the r.h.s. of Eq. (108) can be treated perturbatively when one assumes that the activities of nearby modules are close, and that the noise is not too strong. As before, we suppose that a good starting approximation for the currents $\mathbf{I}_n(t) = (I_{E,n}(t), I_{I,n}(t))$ which describe the dynamics of the nth module is provided by $\mathbf{I}^{(0)}(t + \phi_n) = (I_E^{(0)}(t + \phi_n), I_I^{(0)}(t + \phi_n))$, namely the vector of currents on the deterministic limit cycle at a particular phase $\phi_n$. Linearization of the chain dynamics around this approximation gives, similarly to Eqs. (38,39),

$$ \frac{d\delta\mathbf{I}_n(t + \phi_n)}{dt} = \mathcal{L}[\mathbf{I}^{(0)}(t + \phi_n)] \cdot \delta\mathbf{I}_n(t + \phi_n) + \mathbf{F}_n^{\text{forc}}(t + \phi_n) - \frac{d\mathbf{I}^{(0)}(t + \phi_n)}{d\phi_n} \frac{d\phi_n}{dt} \quad (109) $$

where the matrix $\mathcal{L}$ is given by Eq. (41). The forcing of the nth module $\mathbf{F}_n^{\text{forc}}$ can be decomposed as a sum of inputs coming from the mean activities, $\mathbf{I}^{(0)}(t + \phi_m)$, of nearby modules and of fluctuations of their activities,

$$ \mathbf{F}_n^{\text{forc}} = \mathbf{F}_n^{\text{f,mean}} + \mathbf{F}_n^{\text{f,stoch}}. \quad (110) $$

With the help of Eq. (108), one obtains

$$ \begin{aligned} \mathbf{F}_n^{\text{f,mean}} &= \left( \frac{\frac{w_{EE}}{\tau_E[I_E^{(0)}(t+\phi_n)]}}{\frac{w_{IE}}{\tau_I[I_I^{(0)}(t+\phi_n)]}} \right) \sum_{m=1}^{L} C(n,m) \left\{ \Phi_E[I_E^{(0)}(t+\phi_m)] - \Phi_E[I_E^{(0)}(t+\phi_n)] \right\} \\ &\approx \left( \frac{\frac{w_{EE}}{\tau_E[I_E^{(0)}(t+\phi_n)]}}{\frac{w_{IE}}{\tau_I[I_I^{(0)}(t+\phi_n)]}} \right) \sum_{m=1}^{L} C(n,m)(\phi_m - \phi_n) \Phi'_E[I_E^{(0)}(t+\phi_n)] \frac{dI_E^{(0)}(t+\phi_n)}{dt} \end{aligned} \quad (111) $$

and

$$ \begin{aligned} \mathbf{F}_n^{\text{f,stoch}} &= \left( \frac{\frac{w_{EE}}{\tau_E[I_E^{(0)}(t+\phi_n)]}}{\frac{w_{IE}}{\tau_I[I_I^{(0)}(t+\phi_n)]}} \right) \sum_{m=1}^{L} C(n,m) \sqrt{\frac{\Phi_E[I_E^{(0)}(t+\phi_m)]}{N_E}} \xi_{E,m}(t) \\ &\approx \left( \frac{\frac{w_{EE}}{\tau_E[I_E^{(0)}(t+\phi_n)]}}{\frac{w_{IE}}{\tau_I[I_I^{(0)}(t+\phi_n)]}} \right) \sum_{m=1}^{L} C(n,m) \sqrt{\frac{\Phi_E[I_E^{(0)}(t+\phi_m)]}{N_E}} \xi_{E,m}(t). \end{aligned} \quad (112) $$

In Eq. (111), we have approximated in a linear way the difference of activities between the modules $n$ and $m$. This supposes that the phase difference between the two modules stays close when their distance is small enough for $C(n, m)$ not to be negligible. (Note that this approximation simplifies the phase equation below but it is not required; without it, the linear phase difference between module $n$ and $m$ would be replaced by a nontrivial function of the phase difference, as in Eq. (60) for the two-module case.) In Eq. (112), we have neglected the phase difference between the two modules in the amplitudes of the stochastic terms since these terms are already treated perturbatively.

Substitution of the forcing term (110) along with the expressions (111,112) in Eq. (50) provides coupled equations for the oscillation phases along the chain,

$$ \frac{d\phi_n}{dt} = D_{\phi}^{EI} \sum_{m=1, \dots, L} C(n, m)[\phi_m(t) - \phi_n(t)] + \eta_n(t), \quad (113) $$

with the constant $D_{\phi}^{EI}$ given by Eq. (63). The value of $D_{\phi}^{EI}$ for the reference parameters is provided in Table 1. The stochastic terms $\eta_n$ in Eq. (113) are Gaussian and are fully characterized by their correlation functions

$$ \langle \eta_n(t) \eta_m(t') \rangle = \delta(t-t') \left[ \frac{D_I}{N_I} \delta_{n,m} + \frac{D_E}{N_E} S(n-m) \right] \quad (114) $$
with

$$ 
\begin{aligned}
S(n-m) &= \sum_{l=1, \dots, L} C(n,l)C(m,l) \\
&= A_{\lambda}^{2} \exp(-\lambda|n-m|) [\coth(\lambda) + |n-m|],
\end{aligned} 
\quad (115) $$

where the second equality in Eq. (115) holds for the exponentially decreasing coupling (20) and where, for simplicity, the expression has been written for a long chain (i.e., neglecting $\exp(-L\lambda)$ for $L \gg 1/\lambda$). Eqs. (53,54) give the expressions of the constants $D_E$ and $D_I$ that quantify the phase diffusion of an oscillatory E-I module due to stochastic fluctuations.

## Characterization of the E-I module chain stochastic activity profile

### Long-wavelength modes

As described by Eq. (113), the phases of the different oscillating E-I modules of the spatially extended network are stochastic quantities. Since Eq. (113) is linear and the stochastic terms are Gaussian, the phases also form a Gaussian field that can be fully characterized by its correlation functions. Since the chain is translation invariant, these are most easily computed by writing the phases $\phi_n$ and the stochastic terms $\eta_n$ in Fourier space,

$$ \phi_n(t) = \frac{1}{\sqrt{L}} \sum_{k=0, \dots, L-1} \tilde{\phi}_{q^{(k)}}(t) \exp(iq^{(k)}n), \quad (116) $$

$$ \eta_n(t) = \frac{1}{\sqrt{L}} \sum_{k=0, \dots, L-1} \tilde{\eta}_{q^{(k)}}(t) \exp(iq^{(k)}n), \quad (117) $$

where the L wavevectors are $\{q^{(k)} = 2\pi k/L, k = 0, \dots, L-1\}$. The stochastic terms $\tilde{\eta}_{q^{(k)}}$ are Gaussian as linear sums of the $\eta_n(t)$ and their correlations are obtained from Eqs. (114, 115) as

$$ \langle \tilde{\eta}_{q^{(k)}}(t) \tilde{\eta}_{q^{(l)}}(t') \rangle = \delta(t-t') \delta_{k+l \equiv 0 \pmod{L}} \left[ \frac{D_I}{N_I} + \frac{D_E}{N_E} \tilde{C}_{\lambda}^2(q^{(k)}), \right] \quad (118) $$

where we have made use of the identity

$$ \sum_n \exp(iqn)S(n) = \tilde{C}_{\lambda}^2(q). \quad (119) $$

The Fourier transform of the coupling function, $\tilde{C}_{\lambda}(q)$, is given in Eq. (106). One can note that its normalization implies that $\tilde{C}_{\lambda}(q) \to 1$ in the long-wavelength limit $q \to 0$. Replacement of expressions (116) and (117) into Eq. (113) gives

$$ \frac{d\tilde{\phi}_{q^{(k)}}}{dt} = D_{\phi}^{EI} [\tilde{C}_{\lambda}(q^{(k)}) - 1] \tilde{\phi}_{q^{(k)}}(t) + \tilde{\eta}_{q^{(k)}}(t). \quad (120) $$

Note that the deterministic part of Eq. (120) exactly corresponds to our previous Eq. (107) for the relaxation of the long-wavelength modes.

For $k \neq 0$, the above equation can easily be integrated to obtain $\tilde{\phi}_{q^{(k)}}(t)$,

$$ \tilde{\phi}_{q^{(k)}}(t) = \int_{-\infty}^{t} dt' \exp \left\{ D_{\phi}^{EI} \left[ \tilde{C}_{\lambda}(q^{(k)}) - 1 \right] (t - t') \right\} \tilde{\eta}_{q^{(k)}}(t'). \quad (121) $$
Therefore, the mean square amplitude for the E-I module oscillation phases for a modulation of wavevector $q^{(k)} \neq 0$ is given by

$$ \langle |\tilde{\phi}_{q^{(k)}}(t)|^2 \rangle = \frac{D_I/N_I + D_E/N_E \tilde{C}_{\lambda}^2(q^{(k)})}{2D_{\phi}^{EI}[1 - \tilde{C}_{\lambda}(q^{(k)})]}, \quad (122) $$

where we have used the expression (118) for the noise correlation. Eq. (122) is the analog for the chain of modules of Eq. (16) in the two-module case. It is relatively straightforward to transform the $\langle |\tilde{\phi}_{q^{(k)}}(t)|^2 \rangle$ into an expression for the amplitudes of the firing rate fluctuations $\langle |\tilde{r}_{q^{(k)}}(t)|^2 \rangle$ for $k \neq 0$, along the lines of the previous calculations. Namely, we find

$$ \begin{aligned} \langle |\tilde{r}_{q^{(k)}}(t)|^2 \rangle_T &= \frac{1}{L} \sum_{n,m} \exp[iq^{(k)}(m-n)] \langle r_{n,E}(t) r_{m,E}(t) \rangle_T \\ &\approx \frac{1}{L} \sum_{n,m=1}^{L} \exp[iq^{(k)}(m-n)] \sum_{l=-\infty}^{+\infty} |\tilde{r}_l|^2 \exp\left(-\frac{2\pi^2 l^2}{T} \langle [\phi_n(t) - \phi_m(t)]^2 \rangle\right) \end{aligned} \quad (123) $$

with

$$ \langle [\phi_n - \phi_m]^2 \rangle = -2 \langle \phi_n \phi_m \rangle = -2 \sum_{j=0}^{L-1} \exp[iq^{(j)}(m-n)] \langle |\tilde{\phi}_{q^{(j)}}(t)|^2 \rangle. \quad (124) $$

Expression (123) can be expected to hold at large wavelengths, or small $q$, while it will fail at shorter wavelengths where the effective coupling is stronger. It is shown (including an additional contribution $\langle r_E^{(0)} \rangle_T / (N_E \Delta t)$ due to Poisson sampling of spikes) alongside a more precise estimate valid for all couplings (see below) and FAT rate model results in Fig. 7I.

When $q \to 0$, the coupling strength vanishes. Eq. (122) shows that the mean square amplitudes of the long-wavelength modes then diverge as $q^{-2}$,

$$ \langle |\tilde{\phi}_q(t)|^2 \rangle \sim \left[ \frac{D_I}{N_I} + \frac{D_E}{N_E} \right] \frac{\cosh(\lambda) - 1}{D_{\phi}^{EI} q^2}, \quad q \to 0. \quad (125) $$

In this limit, the phase equation (121) reduces to the classical Edwards-Wilkinson equation [69],

$$ \partial_t \phi = D_{lw} \partial_{xx} \phi + \eta(x, t), \quad \langle \eta(x, t) \eta(x', t') \rangle = D_N \delta(x - x') \delta(t - t'), \quad (126) $$

where $D_N = D_I/N_I + D_E/N_E$ is the local module noise amplitude (Eq. (9)) and the long-wavelength diffusion constant $D_{lw}$ is given by

$$ D_{lw} = \frac{D_{\phi}^{EI}}{2[\cosh(\lambda) - 1]}. \quad (127) $$

## Fluctuations of arbitrary wavelengths

Eqs. (113) and (122) very explicitly quantify the dynamics of long-wavelength modes and their average amplitude. Their validity depends however on the fact that for $q \to 0$, one eigenvalue of the matrix $L_q^{EI}$ is very close to 1 and the corresponding (phase) mode dominates the fluctuation of activity at the corresponding wavelength $2\pi/q$. For modulations of smaller wavelengths, the modules are more strongly coupled. The two modes of $L_q^{EI}$ significantly contribute to the fluctuation of activity and Eq. (122) looses
its accuracy. In this case, a more precise estimate of the fluctuation amplitude is obtained, for small noise, by linearizing around the fully synchronized state. This extends to a chain of E-I modules our previous analysis for two coupled E-I modules (Eqs. (83)-(96)).

As in the above analysis of the stability of the synchronized state for the chain, we consider perturbed evolutions for the L modules around the fully synchronized states $\mathbf{I}_n(t) = \mathbf{I}^{(0)}(t) + \delta\mathbf{I}_n(t)$, $n = 1, \dots, L$. The linear evolution of the perturbations $\delta\mathbf{I}_n$ is found to be

$$ \frac{d\delta\mathbf{I}_n}{dt} = \sum_{m=1,\dots,L} \mathcal{L}_{nm}^{EI}[\mathbf{I}^{(0)}(t)]\delta\mathbf{I}_m + \mathbf{B}_n \quad (128) $$

with the stochastic terms $\mathbf{B}_n$ given by,

$$ \mathbf{B}_n = \left( \begin{array}{c} \frac{1}{\tau_E[I_E^{(0)}(t)]} \left[ w_{EE} \sqrt{\frac{\Phi_E[I_E^{(0)}(t)]}{N_E}} \sum_m C(n-m) \xi_{E,m}(t) - w_{EI} \sqrt{\frac{\Phi_I[I_I^{(0)}(t)]}{N_I}} \xi_{I,n}(t) \right] \\ \frac{1}{\tau_I[I_I^{(0)}(t)]} \left[ w_{IE} \sqrt{\frac{\Phi_E[I_E^{(0)}(t)]}{N_E}} \sum_m C(n-m) \xi_{E,m}(t) \right] \end{array} \right) \quad (129) $$

As above (Eq. (101),(117)), translation invariance allows us to simplify the coupled equations (128) by introducing the Fourier representations

$$ \delta \mathbf{I}_n(t) = \frac{1}{\sqrt{L}} \sum_{k=0, \dots, L-1} \delta \tilde{\mathbf{I}}_{q(k)}(t) \exp(iq^{(k)}n), \quad (130) $$

$$ \xi_{\nu,n}(t) = \frac{1}{\sqrt{L}} \sum_{k=0, \dots, L-1} \tilde{\xi}_{\nu,q^{(k)}}(t) \exp(iq^{(k)}n), \quad \nu = E, I. \quad (131) $$

The noise Fourier components obey,

$$ \langle \tilde{\xi}_{\nu,q^{(k)}}(t) \tilde{\xi}_{\nu',q^{(l)}}(t') \rangle = \delta(t-t') \delta_{k+l \equiv 0 \pmod{L}} \delta_{\nu,\nu'} \quad (132) $$

Replacement of Eqs. (130,131) in Eq. (129) gives decoupled dynamics for each Fourier mode of the currents $\delta\tilde{\mathbf{I}}_{q(k)}(t)$, generalizing Eq. (102) to account for stochasticity,

$$ \frac{\delta \tilde{\mathbf{I}}_{q(k)}}{dt} = \mathcal{L}_{q(k)}^{EI}[\mathbf{I}^{(0)}(t)] \delta \tilde{\mathbf{I}}_{q(k)} + \delta \tilde{\mathbf{B}}_{q(k)}(t). \quad (133) $$

The stochastic terms $\delta\tilde{\mathbf{B}}_{q(k)}$ are the Fourier components of the $\mathbf{B}_n$ (Eq. (129)),

$$ \delta \tilde{\mathbf{B}}_{q(k)} = \left( \begin{array}{c} \frac{1}{\tau_E[I_E^{(0)}(t)]} \left[ w_{EE} \sqrt{\frac{\Phi_E[I_E^{(0)}(t)]}{N_E}} \tilde{C}_\lambda(q^{(k)}) \tilde{\xi}_{E,q^{(k)}}(t) - w_{EI} \sqrt{\frac{\Phi_I[I_I^{(0)}(t)]}{N_I}} \tilde{\xi}_{I,q^{(k)}}(t) \right] \\ \frac{1}{\tau_I[I_I^{(0)}(t)]} \left[ w_{IE} \sqrt{\frac{\Phi_E[I_E^{(0)}(t)]}{N_E}} \tilde{C}_\lambda(q^{(k)}) \tilde{\xi}_{E,q^{(k)}}(t) \right] \end{array} \right). \quad (134) $$

Eq. (133) and (134) generalize to arbitrary Fourier modes the previous Eq. (86) and (87) for the antisymmetric mode in the two-module case with the with the already-noted replacement of $(1-2f_{tr})$ by the Fourier transform $\tilde{C}_\lambda(q)$ (Eq. (106)). The amplitudes of the Fourier mode modulations thus read, similarly to Eqs. (95,96),

$$ \langle |\delta \tilde{I}_{E,q^{(k)}}(t)|^2 \rangle = \sum_{j,l=1,2} e_{j,E}^{q(k)}(t) e_{l,E}^{q(k)}(t) \left[ \int_0^t du \, b_{jl}(u) + \frac{\mu_j \mu_l}{1 - \mu_j \mu_l} \int_0^T du \, b_{jl}(u) \right], \quad (135) $$

where the vectors $e_j^q(t)$ are defined for the matrices $\mathcal{L}_q^{EI}(t)$ as the vectors $e_j^A(t)$ for the matrix $\mathcal{L}_A^{EI}(t)$.
The amplitudes of the rate fluctuations are again obtained by multiplying the expression above by $\Phi'_E[I_E^{(0)}(t)]^2$, and averaged over the complete limit cycle we find

$$ \langle |\tilde{r}_{E,q^{(k)}}(t)|^2 \rangle_T = \frac{1}{T} \int_0^T dt \left\{ \sum_{j,l=1,2} \Phi'_E[I_E^{(0)}(t)]^2 e_{j,E}^{q^{(k)}}(t) e_{l,E}^{q^{(k)}}(t) \times \left[ \int_0^t du \, b_{jl}(u) + \frac{\mu_j \mu_l}{1 - \mu_j \mu_l} \int_0^T du \, b_{jl}(u) \right] \right\} + \frac{\langle r_E^{(0)} \rangle_T}{N_E \Delta T} . \quad (136) $$

Here, we explicitly added the contribution that arises due to Poissonian sampling of the spiking dynamics within bins of size $\Delta T$. The above result is compared to the stochastic FAT rate model in Fig. 7I.

## Supporting information

S1 File. Tabulated values for $\Phi_\sigma(I)$ and $\tau^{(fat)}(I)$.

S1 Fig. The FAT rate model for EIF neurons. Plots of the f-I curve and the fitted adaptive timescale (FAT) used throughout the manuscript, and an assessment of the FAT rate model's capacity to account for time-varying input current.

S2 Fig. Oscillatory E-I module dynamics for reference parameters B and C. Comparison of spiking network and FAT rate model dynamics for different synaptic strengths that those used throughout the manuscript.

S3 Fig. Dynamical regimes of two coupled E-I modules for reference parameters B and C. Phase diagrams of the dynamical regimes of two E-I modules coupled by long-range $E \to E$ excitatory connectivity, for different synaptic strengths that those used throughout the manuscript.

S4 Fig. Average phase gradients for a chain of E-I modules. The probability distribution of extended phase gradients $\nabla\phi_x = (\phi_{x+5} - \phi_x)/5$ that exceed the standard deviation of locally contributing phase differences $\nabla\phi_{x+i} = \phi_{x+i+1} - \phi_{x+i}$, $i = 0, \dots, 4$.

S5 Fig. Phase gradients in a chain of E-I modules with decreased long-range excitation. Histograms of local phase differences and extended phase gradients for different space decay constants $\lambda$ and network sizes $N$ (as well as both types of long-range connectivity studied), when the contribution of long-range excitation to the local dynamics is diminished by a factor 2.

## Acknowledgments

We would like to thank Joran Deschamps for performing preliminary rate-model simulations and computations on synchronization of E-I modules during an M2 internship. We are also grateful to Frédéric Chavane and Thomas Brochier for instructive comments as well as to Josef Ladenbauer and Nicolas Brunei for drawing our attention respectively to ref. [50] and ref. [13]. We also wish to thank Bard Ermentrout for drawing our attention to ref. [22] which had escaped our notice. AK acknowledges the support of Fondation de la Recherche Médicale (FRM) through a doctoral FRM fellowship.
References

1. Buzsaki G. Rhythms of the Brain. Oxford University Press; 2006.

2. Wang XJ. Neurophysiological and computational principles of cortical rhythms in cognition. Physiological Reviews. 2010;90(3):1195-1268.

3. Fries P. Neuronal gamma-band synchronization as a fundamental process in cortical computation. Annual Review of Neuroscience. 2009;32:209-224.

4. Sanes JN, Donoghue JP. Oscillations in local field potentials of the primate motor cortex during voluntary movement. Proceedings of the National Academy of Sciences. 1993;90(10):4470-4474.

5. Engel AK, Fries P. Beta-band oscillations—signalling the status quo? Current Opinion in Neurobiology. 2010;20(2):156-165.

6. Muller L, Chavane F, Reynolds J, Sejnowski TJ. Cortical travelling waves: mechanisms and computational principles. Nat Rev Neurosci. 2018;19(5):255-268.

7. Prechtl J, Cohen L, Pesaran B, Mitra P, Kleinfeld D. Visual stimuli induce waves of electrical activity in turtle cortex. Proceedings of the National Academy of Sciences. 1997;94(14):7621-7626.

8. Huang X, Xu W, Liang J, Takagaki K, Gao X, Wu Jy. Spiral wave dynamics in neocortex. Neuron. 2010;68(5):978-990.

9. Veit J, Hakim R, Jadi MP, Sejnowski TJ, Adesnik H. Cortical gamma band synchronization through somatostatin interneurons. Nature Neuroscience. 2017;20(7):951-959.

10. Rubino D, Robbins KA, Hatsopoulos NG. Propagating waves mediate information transfer in the motor cortex. Nature Neuroscience. 2006;9(12):1549.

11. Takahashi K, Kim S, Coleman TP, Brown KA, Suminski AJ, Best MD, et al. Large-scale spatiotemporal spike patterning consistent with wave propagation in motor cortex. Nature Communications. 2015;6:7169.

12. Denker M, Zehl L, Kilavik BE, Diesmann M, Brochier T, Riehle A, et al. LFP beta amplitude is linked to mesoscopic spatio-temporal phase patterns. Scientific Reports. 2018;8(1):5200.

13. Rule ME, Vargas-Irwin C, Donoghue JP, Truccolo W. Phase reorganization leads to transient β-LFP spatial wave patterns in motor cortex during steady-state movement preparation. Journal of Neurophysiology. 2018;119(6):2212-2228.

14. Takahashi K, Saleh M, Penn RD, Hatsopoulos NG. Propagating waves in human motor cortex. Frontiers in Human Neuroscience. 2011;5.

15. Huntley GW, Jones EG. Relationship of intrinsic connections to forelimb movement representations in monkey motor cortex: a correlative anatomic and physiological study. Journal of Neurophysiology. 1991;66(2):390-413.

16. Hao Y, Riehle A, Brochier TG. Mapping Horizontal Spread of Activity in Monkey Motor Cortex Using Single Pulse Microstimulation. Frontiers in Neural Circuits. 2016;10.
17. Capaday C, Ethier C, Brizzi L, Sik A, van Vreeswijk C, Gingras D. On the nature of the intrinsic connectivity of the cat motor cortex: evidence for a recurrent neural network topology. Journal of Neurophysiology. 2009;102(4):2131-2141.

18. Pikovsky A, Kurths J, Rosenblum M, Kurths J. Synchronization: a universal concept in nonlinear sciences. vol. 12. Cambridge University Press; 2003.

19. Kuramoto Y. Chemical oscillations, waves, and turbulence. vol. 19. Springer-Verlag, Berlin; 1984.

20. Wilson HR, Cowan JD. Excitatory and inhibitory interactions in localized populations of model neurons. Biophysical Journal. 1972;12(1):1-24.

21. Ermentrout GB, Kopell N. Multiple pulse interactions and averaging in systems of coupled neural oscillators. Journal of Mathematical Biology. 1991;29(3):195-217.

22. Borisyuk GN, Borisyuk RM, Khibnik AI, Roosevelt D. Dynamics and bifurcations of two coupled neural oscillators with different connection types. Bulletin of Mathematical Biology. 1995;57(6):809-840.

23. Ermentrout GB, Kleinfeld D. Traveling electrical waves in cortex: insights from phase dynamics and speculation on a computational role. Neuron. 2001;29(1):33-44.

24. Hoppensteadt FC, Izhikevich EM. Weakly connected neural networks. vol. 126. Springer Science & Business Media; 2012.

25. Ashwin P, Coombes S, Nicks R. Mathematical frameworks for oscillatory network dynamics in neuroscience. The Journal of Mathematical Neuroscience. 2016;6(1):2.

26. Amari Si. Dynamics of pattern formation in lateral-inhibition type neural fields. Biological Cybernetics. 1977;27(2):77-87.

27. Ermentrout B. Neural networks as spatio-temporal pattern-forming systems. Reports on Progress in Physics. 1998;61(4):353.

28. Coombes S. Waves, bumps, and patterns in neural field theories. Biological Cybernetics. 2005;93(2):91-108.

29. Bressloff PC. Spatiotemporal dynamics of continuum neural fields. Journal of Physics A: Mathematical and Theoretical. 2011;45(3):033001.

30. Brunel N, Hakim V. Sparsely synchronized neuronal oscillations. Chaos: An Interdisciplinary Journal of Nonlinear Science. 2008;18(1):015113.

31. Brunel N. Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. Journal of Computational Neuroscience. 2000;8(3):183-208.

32. Brunel N, Wang XJ. What determines the frequency of fast network oscillations with irregular neural discharges? I. Synaptic dynamics and excitation-inhibition balance. Journal of Neurophysiology. 2003;90(1):415-430.

33. Börgers C, Kopell N. Synchronization in networks of excitatory and inhibitory neurons with sparse, random connectivity. Neural Computation. 2003;15(3):509-538.

34. Battaglia D, Brunel N, Hansel D. Temporal decorrelation of collective oscillations in neural networks with local inhibition and long-range excitation. Physical Review Letters. 2007;99(23):238106.
35. Battaglia D, Hansel D. Synchronous chaos and broad band gamma rhythm in a minimal multi-layer model of primary visual cortex. PLoS Computational Biology. 2011;7(10):e1002176.

36. Palmigiano A, Geisel T, Wolf F, Battaglia D. Flexible information routing by transient synchrony. Nature Neuroscience. 2017;.

37. Kriener B, Helias M, Rotter S, Diesmann M, Einevoll GT. How pattern formation in ring networks of excitatory and inhibitory spiking neurons depends on the input current regime. Frontiers in Computational Neuroscience. 2014;7:187.

38. Senk J, Korvasová K, Schuecker J, Hagen E, Tetzlaff T, Diesmann M, et al. Conditions for traveling waves in spiking neural networks. arXiv preprint arXiv:180106046. 2018;.

39. Brunel N, Hakim V. Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. Neural Computation. 1999;11(7):1621-1671.

40. Montbrió E, Pazó D, Roxin A. Macroscopic description for networks of spiking neurons. Physical Review X. 2015;5(2):021028.

41. Dumont G, Ermentrout GB, Gutkin B. Macroscopic phase-resetting curves for spiking neural networks. Physical Review E. 2017;96(4):042311.

42. Dumont G, Gutkin B. Macroscopic phase resetting-curves determine oscillatory coherence and signal transfer in inter-coupled neural circuits. PLoS Computational Biology. 2019;15(5):e1007019.

43. Vreeswijk Cv, Sompolinsky H. Chaotic balanced state in a model of cortical circuits. Neural Computation. 1998;10(6):1321-1371.

44. Rosenbaum R, Doiron B. Balanced networks of spiking neurons with spatially dependent recurrent connections. Physical Review X. 2014;4(2):021039.

45. Rosenbaum R, Smith MA, Kohn A, Rubin JE, Doiron B. The spatial structure of correlated neuronal variability. Nature Neuroscience. 2017;20(1):107.

46. Shriki O, Hansel D, Sompolinsky H. Rate models for conductance-based cortical neuronal networks. Neural Computation. 2003;15(8):1809-1841.

47. Ermentrout B. Reduction of conductance-based models with slow synapses to neural nets. Neural Computation. 1994;6(4):679-695.

48. Ostojic S, Brunel N. From spiking neuron models to linear-nonlinear models. PLoS Computational Biology. 2011;7(1):e1001056.

49. Schaffer ES, Ostojic S, Abbott LF. A complex-valued firing-rate model that approximates the dynamics of spiking networks. PLoS Computational Biology. 2013;9(10):e1003301.

50. Augustin M, Ladenbauer J, Baumann F, Obermayer K. Low-dimensional spike rate models derived from networks of adaptive integrate-and-fire neurons: comparison and implementation. PLoS Computational Biology. 2017;13(6):e1005545.

51. Yoshimura K, Arai K. Phase reduction of stochastic limit cycle oscillators. Physical Review Letters. 2008;101(15):154101.

52. Goldobin DS, Teramae Jn, Nakao H, Ermentrout GB. Dynamics of limit-cycle oscillators subject to general noise. Physical Review Letters. 2010;105(15):154101.
53. Fourcaud-Trocmé N, Hansel D, Van Vreeswijk C, Brunel N. How spike generation mechanisms determine the neuronal response to fluctuating inputs. *Journal of Neuroscience*. 2003;23(37):11628–11640.

54. Badel L, Lefort S, Brette R, Petersen CC, Gerstner W, Richardson M.J. Dynamic IV curves are reliable predictors of naturalistic pyramidal-neuron voltage traces. *Journal of Neurophysiology*. 2008;99(2):656–666.

55. Roxin A, Compte A. Oscillations in the bistable regime of neuronal networks. *Physical Review E*. 2016;94(1):012410.

56. Mountcastle VB. The columnar organization of the neocortex. *Brain*. 1997;120 (Pt 4):701–722.

57. Defelipe J, Markram H, Rockland KS. The neocortical column. *Front Neuroanat*. 2012;6:22.

58. Gardiner C. Stochastic methods. Springer Series in Synergetics (Springer-Verlag, Berlin, 2009). 1985;.

59. Nauhaus I, Busse L, Carandini M, Ringach DL. Stimulus contrast modulates functional connectivity in visual cortex. *Nature Neuroscience*. 2009;12(1):70.

60. Adesnik H, Scanziani M. Lateral competition for cortical space by layer-specific horizontal circuits. *Nature*. 2010;464(7292):1155.

61. Huang X, Elyada YM, Bosking WH, Walker T, Fitzpatrick D. Optogenetic assessment of horizontal interactions in primary visual cortex. *Journal of Neuroscience*. 2014;34(14):4976–4990.

62. Chavane F, Sharon D, Jancke D, Marre O, Frégnac Y, Grinvald A. Lateral spread of orientation selectivity in V1 is controlled by intracortical cooperativity. *Frontiers in Systems Neuroscience*. 2011;5:4.

63. Levy RB, Reyes AD. Spatial profile of excitatory and inhibitory synaptic connectivity in mouse primary auditory cortex. *Journal of Neuroscience*. 2012;32(16):5609–5619.

64. Roberts JA, Gollo LL, Abeysuriya RG, Roberts G, Mitchell PB, Woolrich MW, et al. Metastable brain waves. *Nature Communications*. 2019;10(1):1056.

65. Buice MA, Chow CC. Beyond mean field theory: statistical field theory for neural networks. *Journal of Statistical Mechanics: Theory and Experiment*. 2013;2013(03):P03003.

66. Schwalger T, Deger M, Gerstner W. Towards a theory of cortical columns: From spiking neurons to interacting neural populations of finite size. *PLoS Computational Biology*. 2017;13(4):e1005507.

67. Dumont G, Payeur A, Longtin A. A stochastic-field description of finite-size spiking neural networks. *PLoS Computational Biology*. 2017;13(8):1–34. doi:10.1371/journal.pcbi.1005691.

68. Cross MC, Hohenberg PC. Pattern formation outside of equilibrium. *Reviews of Modern Physics*. 1993;65(3):851.

69. Barabási AL, Stanley HE. Fractal concepts in surface growth. *Cambridge university press*; 1995.
70. Goodman DF, Brette R. The brian simulator. Frontiers in Neuroscience. 2009;3(2):192.

71. Richardson M. Firing-rate response of linear and nonlinear integrate-and-fire neurons to modulated current-based and conductance-based synaptic drive. Physical Review E. 2007;76(2).

72. Bender CM, Orszag SA. Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory. Springer Science & Business Media; 2013.

73. Ince EL. Ordinary Differential Equations. Dover Publications, Inc., N.Y.; 1956.
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Fig S1. The FAT rate model for EIF neurons.** (A) The f-I curve $\Phi_{\sigma}(I)$ for the EIF neurons used in the rate-model formulation (numerically evaluated following [71]). (B) The fitted adaptive timescale $\tau^{(FAT)}(I)$ determined by fitting the linear firing-rate response and used throughout the manuscript (thick blue line). For comparison, we also show the adaptive timescale proposed by Ostojic and Brunel (OB) [48] (thin orange line). (C) Examples of the fit of the analytically determined firing-rate response with the Fourier transformation of an exponential kernel in the time domain (see Methods). (D) Dynamics of a population of (uncoupled) EIF neurons and the two adaptive rate models upon injection of a sinusoidal current $I_1 \sin(2\pi ft)$ for $f = 17$ Hz and different values of $I_1$. (E,F) Correlation between the network activity and the two adaptive rate models (FAT, thick blue; OB, thin orange) as a function of the amplitude $I_1$ of the injected current for a fixed frequency $f = 17$ Hz (E), and as a function of frequency $f$ for a fixed amplitude $I_1 = 3$ mV (F).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig S2. Oscillatory E-I module dynamics for reference parameters B and C. (A,B)**
The deterministic limit cycle dynamics for a single E-I module described by the FAT rate model with synaptic strengths corresponding to reference case B (see Fig. 1A) is shown together with spiking network simulations with $N = 10^6$. The synaptic strengths are given by $w_{EE} = 1.6$ mVs, $w_{IE} = 1.6$ mVs, and $w_{EI} = 0.8$ mVs. (C) A comparison of the stochastic FAT rate model with spiking network simulations with $N = 10^4$ for the same parameters. (D-E) Analogous plots for the reference case C, with $w_{EE} = 1.76$ mVs, $w_{IE} = 3.2$ mVs, and $w_{EI} = 0.4$ mVs.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig S3.** Dynamical regimes of two coupled E-I modules for reference parameters **B** and **C**. The analog of Fig. 2B for the two other reference parameters shown on Fig. 1A. The phase diagrams for case B ($w_{EE}$ = 1.6 mVs, $w_{IE}w_{EI}$ = 1.28 mV$^2$s$^2$) and case C ($w_{EE}$ = 1.76 mVs, $w_{IE}w_{EI}$ = 1.28 mV$^2$s$^2$) are shown on panels (A) and (B), respectively.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig S4.** Average phase gradients for a chain of E-I modules. The probability distribution of extended phase gradients $\overline{\nabla\phi}_x$ that exceed the standard deviation of locally contributing phase differences (see *Methods*), for different values of $\lambda$ and $N$, for $E \to E$ connectivity (A) and $E \to E, I$ connectivity (B).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig S5.** Phase gradients in a chain of E-I modules with decreased long-range excitation. (A,B) Simulation of the stochastic rate model (FAT) for a chain with $E \rightarrow E, I$ connectivity, a space constant $\lambda = 0.33$, and network size $N = 10^4$. In contrast to Fig. 7D,E, the contribution of the long-range connectivity to the total recurrent excitatory drive is lowered by a factor 2 at the expense of purely local recurrent excitation. (A) Instantaneous phase as obtained from the Hilbert transform of bandpass-filtered excitatory activity, where the (unfiltered) rates are shown in (B). (C,D) Histograms of the local phase differences for different values of $\lambda$ and $N$ for $E \rightarrow E$ connectivity (C) and $E \rightarrow E, I$ connectivity (D). (E,F) Histograms of the extended phase gradients $\overline{\nabla}\phi$. The averages of the distributions shown in (E,F) are plotted in (G), and are the analog of Fig. 7I for a diminished contribution (by a factor of 2) of long-range connectivity to the total recurrent excitatory drive.</figcaption></figure>