Synchronization, Stochasticity, and Phase Waves in Neuronal Networks With Spatially-Structured Connectivity

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 (> 105), we obtain results similar to previously 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 (~ 104) 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 explanative framework for the observations of traveling waves in the cortex during beta oscillations.

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 5 ), we obtain results similar to previously 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 4 ) 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 explanative 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 (Buzsaki, 2006). Numerous works have been devoted to understanding their origins and characteristics (Wang, 2010) which depend both on the neural structure 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 (Fries, 2009). Beta oscillations (10-30 Hz) are prominent in the motor cortex during movement planning before movement initiation (Sanes and Donoghue, 1993) and have traditionally been assigned a role in movement control while more general roles have also been proposed (Engel and Fries, 2010).
Several experimental results point out the need to model and analyze the spatial organization of oscillatory activity (Muller et al., 2018). An early study (Prechtl et al., 1997), 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 (Huang et al., 2010). The underlying mechanisms and specific cells involved in the synchronization of two distant regions have more recently been investigated using optogenetic manipulations in mice (Veit et al., 2017).
A motivating 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 (Rubino et al., 2006;Takahashi et al., 2015;Denker et al., 2018;Rule et al., 2018) and humans (Takahashi et al., 2011). These waves appear to propagate along particular anatomical directions with a typical wavelength of 1 cm and a velocity of about 20 cm · s −1 (Rubino et al., 2006). On the structural side, characterization of the long-range intracortical connectivity in motor cortex has given rise to several quantification endeavors in both monkeys (Huntley and Jones, 1991;Hao et al., 2016) and cats (Capaday et al., 2009) with the strength of excitatory responses to microstimulation decaying on a 2 mm length scale. We would like to connect these different observations and find a proper mechanistic framework to interpret these observed waves.
The study of waves in oscillatory media has been largely based on the analysis of oscillator synchronization which itself has a long history (Pikovsky et al., 2003). Classic mathematical methods for studying synchronization of weakly coupled oscillators have been extended to study oscillations and traveling waves in spatially extended media in physics and chemistry (Kuramoto, 1984). In simple descriptions of neural network dynamics based on rate models introduced by Wilson and Cowan (1972), 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 (Ermentrout and Kopell, 1991;Borisyuk et al., 1995;Ermentrout and Kleinfeld, 2001;Hoppensteadt and Izhikevich, 2012;Ashwin et al., 2016). Rate models with spatiallystructured connectivity have also been extensively studied to analyze pattern formation in a neural network context since Amari's work (Amari, 1977) (see Ermentrout, 1998;Coombes, 2005;Bressloff, 2011 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 (Brunel and Hakim, 2008) 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 (Wang, 2010). 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 (Brunel, 2000;Brunel and Wang, 2003), the so-called Pyramidal-Interneuron Gamma (PING) mechanism (Börgers and Kopell, 2003).
A few studies have considered the impact of structured connectivity on sparsely synchronized oscillations in the spikingnetwork framework. The influence of delays in synaptic transmission was studied in Battaglia et al. (2007) for the synchronization of two neuronal E-I modules oscillating in the high-gamma frequency range, using a classic ratemodel 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 (Battaglia and Hansel, 2011) 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 (Palmigiano et al., 2017) 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 with (Senk et al., 2018) or without (Kriener et al., 2014) transmission delays. Going beyond linear stability is feasible but already somewhat heavy for a single module with unstructured connectivity (Brunel andHakim, 1999, 2008). An exception is the soluble case of deterministic quadratic integrate-and-fire neurons with a wide (Lorentzian) distribution of firing rates (Montbrió et al., 2015) which has been used to study oscillations of an E-I module (Dumont et al., 2017a) as well as synchronization between two weakly-coupled modules (Dumont and Gutkin, 2019). Aside from oscillations, the conditions necessary for the existence of a balanced state (Vreeswijk and Sompolinsky, 1998) in a spatially structured network have been studied (Rosenbaum and Doiron, 2014) 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 (Rosenbaum et al., 2017). 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 distancedependent long-range excitation. We focus on the case where each module oscillates in a sparsely synchronized way (Brunel and Hakim, 2008) 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 has been the aim of several works since the invention of rate models (Wilson and Cowan, 1972). This can be achieved when the collective dynamics is stationary or slowly varying (Shriki et al., 2003), for instance in the presence of slow synapses (Ermentrout, 1994), by an appropriate choice of the steady state input-output transfer function (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 slowly varying collective dynamics by introducing a timescale that depends on the firing rate. These "adaptive" rate models successfully produced a quantitative description of an uncoupled neuronal population submitted to time-varying inputs (Ostojic and Brunel, 2011), as well as of spike synchronization (Schaffer et al., 2013) or oscillations driven by spike-frequency-adaptation (SFA) (Augustin et al., 2017) in a recurrently coupled excitatory neural population.
Building on this progress, we first develop a rate model with an adaptive timescale (Ostojic and Brunel, 2011;Augustin et al., 2017) 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 (Ermentrout and Terman, 2010;Hoppensteadt and Izhikevich, 2012;Ashwin et al., 2016). We show that in spite of the introduction of the adaptive timescale, the model that we use behaves very similarly (Ermentrout and Kopell, 1991;Borisyuk et al., 1995;Ermentrout and Terman, 2010;Hoppensteadt and Izhikevich, 2012) to classic rate models (Wilson and Cowan, 1972).
We then consider the synchronization of neural oscillations between two E-I modules. We find that it 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 (Ermentrout and Kopell, 1991;Borisyuk et al., 1995), whereas synchronization of the oscillations of the two modules is otherwise observed. An additional feature of our results as compared to the ones obtained with classical rate-models, is that, as we show, they 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 (Brunel and Hakim, 1999), we find that finite-size fluctuations can be quantitatively accounted in the adaptive-timescale ratemodel 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. (Pikovsky et al., 2003;Kuramoto, 1984;Ashwin et al., 2016). Classical (Kuramoto, 1984) and more recent results (Yoshimura and Arai, 2008;Goldobin et al., 2010;Ashwin et al., 2016) 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 weaknoise 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.
Finally, we study the dynamics of a chain of oscillatory E-I modules with long-range excitatory coupling that decreases with distance. For large module sizes, we find, as in the twomodule case, that the connectivity properties of long-range excitation play an important role. In particular we find that when long-range excitation only targets excitatory neurons, distant modules oscillate with different phases. Namely, spatial gradients in the module's phases of oscillation appear spontaneously, as well as phase waves. For modules of moderate sizes, the specific connectivity of long-range excitation is found to be less important. Stochastic dephasing creates transient bouts of traveling waves in a chain of modules even when it would be fully synchronized in a deterministic context. We end by obtaining, for 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.
The present work reports results of numerical simulations of spiking networks and rate models and of systematic comparisons of these simulations with analytic mathematical formulas aimed at describing these numerical results. Most of the formulas are applications of results previously derived in the literature to the specific adaptive rate model that we simulate. In the methods section, we provide references to these previous results but also give short derivations of all the formulas that we use in a single framework, for the convenience of the interested reader. We have however tried to structure the manuscript so that these mathematical derivations can be skipped without impairing the understanding of our 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 (Wilson and Cowan, 1972) and in simulations of spiking networks (Brunel, 2000;Brunel and Wang, 2003;Börgers and Kopell, 2003). We choose the neurons in our spiking network simulations to be of the Exponential-Integrate-and-Fire type (EIF) (Fourcaud-Trocmé et al., 2003) (see Equation 25 for their mathematical definition), which have been shown to describe well the dynamics of cortical neurons (Badel et al., 2008).
We wish to use a rate-model description that quantitatively describes oscillations of the spiking E-I module in the SSO regime (Brunel and Hakim, 2008). Mild modifications of classic rate models (Wilson and Cowan, 1972) can give quantitatively accurate accounts of simulations of spiking neurons when spike emission has a strong stochastic component, as is the case in the SSO regime. This was previously demonstrated both for an assembly of independent spiking neurons receiving identical inputs (Ostojic and Brunel, 2011) and for spiking neurons coupled by recurrent excitation (Augustin et al., 2017). Building on these advances, we use the following adaptive timescale ratemodel formulation (Ostojic and Brunel, 2011;Augustin et al., 2017) to represent the activity of a group of EIF neurons: where σ (I) is the f-I curve of the EIF model (Equation 25) for a noisy input with mean I and total noise strength σ . In the following, we refer to the input I as a current but note that it has the dimension of a voltage because it is the input current divided by the neuron leak conductance (see Equation  25). It is shown in Supplementary Figure 1A for the parameters used in the present study. This specific choice of f-I curve allows the rate model to agree with the rate of a network of spiking EIF neurons 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 networkaveraged instantaneous internal inputs. One way to understand the model described by Equations (1) and (2) is to view it (Ostojic and Brunel, 2011) as a generalization of a linear-non-linear model (Chichilnisky, 2001). When τ (I) is constant, the linear Equation (1) can be explicitly integrated and gives that I is related to the input s(t) through the convolution with an exponential kernel. Equation (2) then provides the neuron firing rate as a non-linear function of this linear transform. The response time τ (I) is said to be "adaptive" (Ostojic and Brunel, 2011) because it depends on the current I. It allows the rate model to describe spiking neurons when the population rate response remains close to the stationary f-I curve, without being limited to a slowly varying dynamics (Ermentrout, 1994;Shriki et al., 2003) . The chosen function τ (I) is displayed in Supplementary Figure 1B. 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 Supplementary Figure 1C for examples of the timescale fit).
In agreement with previous works (Ostojic and Brunel, 2011;Augustin et al., 2017), Supplementary Figures 1D-F show 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 Figure 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 Equation (2), with r E = E (I E ) and r I = I (I I ). (Here and in the following we have dropped the explicit noise strength index on (I) for notational simplicity). In all numerical computations, we have chosen E (I) = I (I) = σ (I), τ E (I) = τ I (I) = τ (FAT) (I), with σ (I) and τ (FAT) (I) displayed in Supplementary Figure 1. Note that recurrent inhibition has not been included in Equations (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 parameter regimes (Brunel and Hakim, 1999;Brunel and Wang, 2003) and can also be important to prevent a too high firing rate of interneurons, e.g., in the balanced state (Vreeswijk and Sompolinsky, 1998), a role that is assigned here to the external input I ext I . The stability diagram for the steady non-oscillatory dynamics of the FAT rate model is easily obtained (Wilson and Cowan, 1972;Hoppensteadt and Izhikevich, 2012) (see Methods: Oscillatory Instability for the E-I Module) (see Figure 1A). As for the classic Wilson-Cowan model (Wilson and Cowan, 1972), it displays three regimes. A constant neural activity is a stable fixed point of the dynamics 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 FIGURE 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 (s) E = 5 Hz, r (s) I = 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 Supplementary Figure 2. (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 (A) with w IE = 2 mVs, w EI = 0.32 mVs (the corresponding effective constants are α = 2.33, β = 2.15). Oscillations in a spiking network with a large number of neurons (N = 10 6 , resampled with time bin 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 (Equations 5, 6) (sampled with time bin 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 Ostojic and Brunel (2011) is also shown (solid green). The fit of the analytical expression Equation (56) for the autocorrelation to the adaptive rate model (shown in dashed black) allows to obtain the module's autocorrelation decay time τ D . (F) The decorrelation time τ D shown for spiking modules (solid blue dots), for FAT modules (solid orange line), and as predicted by Equation (10) (dashed black line).
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 Roxin and Compte, 2016 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 (Brunel and Hakim, 2008), a dynamical regime which appears most appropriate to describe beta/low gamma oscillations recorded in vivo. The parameters corresponding to one such point are indicated in Figure 1A (solid black circle) and used as reference for the figures of the present paper. The corresponding limit-cycle oscillations are shown in Figures 1B,C. Two other representative sets of parameter values are also indicated in Figure 1A (solid losange and square) and the corresponding dynamical traces are shown in Supplementary Figure 2.
The rate-model description is compared to simulations of spiking networks (see Methods: Simulations of Spiking Networks) in Figure 1. The rate-model deterministic activity trace accounts well for the network activity when the number N of spiking neurons is very large (N ∼ 10 6 ) and stochastic effects at the level of the population are negligible, as shown in Figures 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 (Mountcastle, 1997;Defelipe et al., 2012), which is estimated to comprise a smaller number of neurons (N ∼ 10 4 − 10 5 ). In this case, simulations show that the network population activity has a significant stochastic component, as seen in Figure 1D. Auto-correlograms of the E population activity display decreasing oscillatory tails reflecting the corresponding dephasing of oscillations ( Figure 1E).
We take into account these stochastic effects in the rate-model description by remembering that, in a sparsely synchronized spiking network (Brunel and Hakim, 2008), 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 Brunel and Hakim (1999) and assume that, in the SSO regime, the probability of an excitatory (respectively inhibitory) spike between the times t and t + t is given by E [I E (t)] t (respectively I [I I (t)] t), and spikes being independently drawn for each neuron in the network. This leads to replacing the deterministic relations r E (t) = E [I E (t)] and r I (t) = 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 (Brunel and Hakim, 1999), where ξ E (t), ξ I (t) are independent white noises ( = 0) and Ito's prescription is used (Gardiner, 1985). Equation (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 ( Figure 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 ′ ) = r E (t)r E (t ′ ) − r E 2 . As shown in Figure 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 τ D , The time τ D increases with the number of neurons in the E-I module as shown in Figure 1F.
Well-known results are available to analytically describe dephasing due to weak noise (Kuramoto, 1984;Yoshimura and Arai, 2008;Goldobin et al., 2010;Ashwin et al., 2016), 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+φ), r I (t+φ)) are oscillatory with r E (t) and r I (t) periodic functions of time, and the phase φ being arbitrary but constant. Note that in the present work, we define the phase φ 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π. (Therefore r E is here a periodic function of the phase φ with period T, instead of being a periodic function of a phase with period 2π.) The phases in these two different conventions are simply related by a multiplicative factor of 2π/T. Stochastic fluctuations in E-I modules with large number of neurons produce a random drift of the phase φ with a diffusive behavior in time, 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, The constants D E and D I can be expressed and computed in terms of integrals along the deterministic limit cycle (Kuramoto, 1984) (Equations 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, Equation (10) predicts that the dephasing time increases linearly with the size N of the E-I module. Expression (10) is displayed in Figure 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.

Long-Range Excitation Targeting Excitatory Neurons
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 (Ermentrout and Kopell, 1991;Borisyuk et al., 1995) obtained for the classic Wilson-Cowan model (Wilson and Cowan, 1972). It matters for synchronization whether long-range excitation targets only excitatory neurons or both excitatory and inhibitory neurons. We thus distinguish the two cases. We begin by analyzing the synchronization properties between two E-I modules with oscillatory activity when long-range excitation only targets excitatory neurons, the "E → E" connectivity case depicted in Figure 2A.
The dynamics of two coupled E-I modules are described in the rate-model framework (3,4) by τ I (I I,1 ) dI I,1 dt = −I I,1 + I ext I + w IE r E,1 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 (Equation 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 Figure 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 ) 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 Figure 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 Figure 1A; the analogous diagrams for the two other points are shown in Supplementary Figure 3). The individual panels (Figures 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 (Kuramoto, 1984;Ermentrout and Kopell, 1991;Hoppensteadt and Izhikevich, 2012;Ermentrout and Kleinfeld, 2001). The synchronization dynamics itself can be characterized by the evolution of the relative phase φ between the oscillations of the two modules. For φ = 0, the two modules oscillate in phase whereas for φ = T/2 they oscillate in antiphase. The evolution of φ is found to be governed by the following dynamics (Kuramoto, 1984;Ermentrout and Kopell, 1991;Hoppensteadt and Izhikevich, 2012;Ermentrout and Kleinfeld, 2001) (see Methods: Synchronization Function for Two Weakly Coupled E-I Modules): The function S E ( φ) is shown in Figure 2C for different synaptic parameter values w IE . For two identical modules, symmetry implies that the phase differences φ = 0 and φ = 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 φ = 0, S ′ E ( φ) > 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 Figure 2B5. On the contrary for w IE < w * IE , both in-phase and antiphase are unstable. The function S E ( φ) displays a new zero at a non-trivial φ with a negative slope ( Figure 2C) and the two modules stably oscillate with this nontrivial phase difference ( Figure 2B6). Such a synchronization regime was also obtained (Ermentrout and Kopell, 1991) for the classic Wilson-Cowan rate model (Wilson and Cowan, 1972).
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 ( Figure 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 ∼ 2 − 7%, depending on the strength w IE of the local excitation on the inhibitory population in each module ( Figure 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 finitephase 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 Figure 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 Figure 2B3. This regime stems The synchronization function governing the evolution of the relative phase of two weakly coupled modules in the E → 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 restabilization in two-modules FAT simulation with E → E, I connectivity (crosses, f lr = 0.001) matches well the prediction of Equations (15) and (66)  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 Figure 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 Figure 2 has been described long ago (Borisyuk et al., 1995) 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 "non-symmetric torus" in Borisyuk et al. (1995) 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, Equation (12) is replaced by (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 → E, I" connectivity case, sketched in Figure 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 Equation (13) with the function S EI ( φ) plotted in Figure 2D. For a small phase difference, S EI ( φ) behaves as with the constant D EI φ depending on the characteristics of the oscillation cycle (see Equation 66 in Methods). The constant D EI φ is found to be positive (see Table 1 for its value at the reference point w IE = 2 mVs). Equation (15) shows that a zero phase lag between the two modules is a dynamically stable configuration, in contrast to our previous results for E → E connectivity (compare with Figure 2C). Equation (15) predicts that an initial activity difference between two oscillating E-I modules should vanish exponentially. As shown in Figure 2E, this is indeed observed in simulations with a rate of decrease that agrees well with the one predicted by Equation (15) (see Figure 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 → 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 Equations (11), (12), or (14) expressed in terms of the respective currents by the stochastic relations Equations (5) and (6) which account for the finite number of neurons in the networks. Figure 3 shows the close correspondence between the ratemodel results and networks with a large number of neurons (N ∼ 10 6 ) for E → E connectivity. The different dynamical regimes for two E-I modules coupled by longrange E → E connectivity are all observed in simulations of networks spiking neurons. Both antiphase (Figure 3A), alternating dominance ( Figure 3B), modulated dominance ( Figure 3C) and full synchrony ( Figure 3D) are observed for parameters that quantitatively correspond to those predicted by the rate-model analysis. Similarly, in the case of E → E, I ) Values of the constants defined in the text for the E-I module limit cycle with the reference parameters (solid circle in Figure 1A) w EE = 1.6 mVs, w EI = 0.32 mVs, w IE = 2.0 mVs. connectivity, two large E-I modules oscillate in full synchrony as predicted by the rate description (not shown). 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 Figure 4.
For E → 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 · 10 5 neurons, whereas it is not apparent for N = 10 4 when stochasticity is stronger (Figures 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 · 10 5 neurons. This is illustrated by simulations with a coupling parameter in the modulated dominance regime (f lr = 0.02, Figure 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 (Figures 4E,F).
The corresponding results for E → E, I connectivity are shown in Figures 4G-L. In 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 (Figures 4G,H). On the other hand, synchronization of the two networks increases when coupling strength increases at fixed N (fixed noise) as shown in One should note that whereas the precise targeting of the long-range excitation plays a crucial role in the instability of inphase oscillations for large networks (N ∼ 10 6 ), it is much less important for smaller (N ∼ 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 → E connectivity) or targets both E and I populations in the other module (E → E, I connectivity) (compare Figure 4B to Figure 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 → 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 (Kuramoto, 1984) as well as neural systems (Ashwin et al., 2016). 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 Equation (13) for the phase difference between the two modules (see Equations 78 and 84 in Methods). (It also provides an estimate to the small shift of oscillation frequency due to small noise (Yoshimura and Arai, 2008;Goldobin et al., 2010) 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, In this expression, the constant D N is simply the amplitude of stochastic fluctuations in a single module (Equation 9), whereas f lr D EI φ quantifies the synchronization strength due to long-range connectivity (Equations 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.
Equation (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, Equation (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 (Equation 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 Figure 5. The mean square difference of excitatory rates between the two modules are quite comparable in the noisy FAT ( Figure 5A) and spiking networks (Figure 5B), again showing the quantitative accuracy of the FAT rate-model description. However, Figure 5 shows that the predicted rate difference fluctuations obtained from Equation (16) ("phase approximation" in Figure 5) are accurate only for very small coupling between the two modules.
The underlying assumption in Equation (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 Figure 6G as a function of the fraction of long-range connection f lr . Indeed, for weak coupling (f lr ≪1), one eigenvalue µ A 1 is close to 1, so that the restoring strength 1 − µ A 1 of the corresponding (phase) mode is much smaller than 1 − µ A 2 , the restoring strength of the other mode. However, an explicit computation shows that µ A 1 quickly decreases as f lr increases (see Equation 76 in Methods). Therefore, the two modes soon play comparable roles as f lr increases.
We developed another approach, valid for small noise but for arbitrary coupling f lr , in order to check this origin of the inaccuracy of Equation (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 In the E → 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.
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 Figure 5, this two-mode computation produces precise estimates of the mean square excitatory rate difference of the two modules, (r 1,E − r 2,E ) 2 , 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 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 (Equations 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, Equation 97, solid black line). Note that for very small couplings, the two-mode calculation coincides with the phase approximation with linearized limit cycle dynamics (Equations 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. a fuller description of spatially-structured connectivity. We consider a chain of E-I modules, with each module's longrange excitation targeting other modules in a distance-dependent way. Interestingly, the effects that we uncovered in the twomodule setting, such as desynchronization at weak coupling for E → 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.

Spontaneous Appearance of Phase Gradients Due to Long-Range Excitation Specifically Targeting Excitatory Neurons
The setting of our analysis is sketched in Figure 6A. E-I modules with locally unstructured connectivity are coupled by long-range excitation. For E → E connectivity, the dynamics of L coupled E-I modules are described in the rate-model framework (3,4) by τ I (I I,n ) dI I,n dt = −I I,n + I ext I + w IE r E,n with C(n, m) = C(|n − m|). For E → E, I connectivity, Equation (18) is replaced by τ I (I I,n ) dI I,n dt = −I I,n + I ext The connection strength decreases monotonically with the distance l between modules (Capaday et al., 2009), 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 (Hao et al., 2016). We specifically take where A λ implements the normalization l=0,··· ,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/λ (i.e., when exp(−λL) ≪ 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 (see e.g., Ali et al., 2016 for a similar analysis). 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  (106)) for perturbations of wavenumber q when the long-range excitation has an exponential profile (Equation 20), for space decay constants λ = 1 (turquoise) and λ = 0.33 (dark blue). For λ ≪ 1, the chain of modules can be approximated by a continuous system. (C) Effective long-range coupling f lr (q; λ) = (1 −C λ (q))/2 for perturbations of wavenumber q. (D) Stability spectrum for perturbations of wavelength q as function of the effective coupling f lr (q; λ) for E → 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 gray 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 L n=1 (r n,E − L m=1 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 → 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 σ 2π/L for different chain lengths L (orange crosses) against the theoretical prediction (Equation (23), black line). Inset: Decay of the spatial modulation of excitatory activity over time (L = 512). For all simulations shown, λ = 1/3. E-I module network considered in the previous section, provided the coupling in the two-module network is chosen as whereC λ (q) is the Fourier transform of the long-range excitation profile C(l) (Equation 106 in Methods).C λ (q) and the corresponding "effective" coupling in the two-module network, f lr (q), are displayed in Figures 6B,C for an exponentially decreasing long-range excitation profile C(l) ∼ exp(−λ|l|). It is seen in Figure 6C that small wavenumbers q correspond to weak effective coupling f lr (q), with f lr (q) → 0 when q → 0. While deriving the exact relation (21)  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 → 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 → 0. The stability diagram displayed in Figure 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 Figure 6D. Inphase 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 spatiallyperiodic perturbations grow when their wavenumber q lies below a critical q * , defined by where λ controls the exponential decrease with distance of the long-range excitation. The second equality in Equation (22) is valid for small λ; the used approximate relation between f lr and q is shown in Figure 6C, it can be seen that it is already very accurate for λ = 1/3. Equation (22) predicts that for a network of sufficient spatial extent L > 2π/q * , spontaneous phase gradients should appear along the chain. This is confirmed by direct numerical simulations of the rate-model equations, shown in Figure 6E. The exponential growth of the modulation of excitatory activity along the chain is shown in Figure 6F, starting from a weakly perturbed fully synchronized state. This directly confirms the linear instability of full synchronization. One further notable feature in Figure 6E is that the phase gradients develop on a spatial scale that is quite long (i.e., ∼30 modules) compared to the spatial scale of the long-range excitation profile C(l) (∼ 1 in Figure 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, Equation (22) gives for the critical wavelength q * ≃ 0.24λ. Similarly, one obtains q m ≃ 0.15λ for the wavenumber q m of the fastest-growing mode, which corresponds to the largest eigenvalue in Figure 6D.
For E → 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 Figure 6G. For all wavelengths, their real part is <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 Figure 6H. Moreover, our stability analysis shows that the longest wavelength modes are the slowest to vanish. Their relaxation rate is obtained (Equation 107 in Methods) as The relaxation rate σ q vanishes when q → 0 and decreases with increasing modulation wavelength. The smallest wavenumber that can be accommodated in a chain of L modules is 2π/L. Figure 6I shows that σ 2π/L as given by Equation (23) precisely describes the approach to full synchronization in our simulations of a E-I module chains of different lengths with long-range E → E, I connectivity.

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 Figure 7, where the results of stochastic rate-model simulations are displayed for both E → E and E → E, I connectivity. Both the phases of oscillation of different modules (Figures 7A,D) and the activity of excitatory neurons (Figures 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 Figures 7C,F for E → E and for E → 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 Figures 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 agrees with a generalization of the twomodule result (Equation 16) for the chain. This result, given in Equation (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 ( Figure 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 Figure 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 = 5, 000, L = 32 for N = 10 4 ), the stochastic FAT rate model reproduces well the results of spiking network simulations (Figures 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, (I) Comparison of the spectra obtained from simulations (FAT, solid lines) and theoretical predictions for E → 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, Equation (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, Equation 123). Network size N and space constant λ as in (G,H). (J) Example of a spatially extended phase gradient (upper panel) that gives rise to 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 ∇φ (averaged over five modules and retained when the absolute average gradient is larger than the standard deviation of local phase differences) with phase propagation velocities v φ 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 λ for E → E connectivity (circles) and E → E, I connectivity (diamonds). activity in different modules peaks at different times and may produce "phase waves" (also described as traveling waves, e.g., in Wang, 2010, as shown in Figure 7J. The propagation velocity v φ 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), where ∇φ 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 (Rubino et al., 2006;Takahashi et al., 2015;Denker et al., 2018;Rule et al., 2018;Takahashi et al., 2011). 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 Figure 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 Supplementary Figure 4, for the two different considered connectivities and two space decay constants λ 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 Figure 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 (Capaday et al., 2009;Nauhaus et al., 2009;Adesnik and Scanziani, 2010;Huang et al., 2014;Hao et al., 2016). For instance, for mouse somatosensory barrel cortex, identifying an E-I module with a barrel column would give d ≃ 175 µm for the distance between neighboring modules (which is also the linear size of an E-I module) and λ ≃ 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 (Adesnik and Scanziani, 2010).
For the motor cortex of adult monkeys, data are available from intra-cortical micro stimulation and recordings with single and multielectrodes (Hao et al., 2016). 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 ≃ 500 µm and N ≃ 2 · 10 4 neurons then gives λ ≃ 0.33 for the long-range connectivity. For this last case, our numerical results give ∇φ ≃ 0.5 ms (E → E, I) or 1 ms (E → E) depending on the precise connectivity ( Figure 7K). With the above estimate of d, Equation (24) then gives a phase velocity v φ ∼1 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., ∇φ ≃ 1.4 ms for N = 5·10 3 and (E → E, I) connectivity (Figure 7K), and correspondingly lowers v φ . 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) ∼ (1 − α)δ 0l + α exp(−λ|l|), where α measures the fraction of truly long-range connectivity. For α = 0.5, the resulting average phase gradients are about two times steeper for all λ and N (see Supplementary Figure 5), implying phase waves that are about two times slower than for α = 1 (compare Supplementary Figure 5E to Figure 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 (Ostojic and Brunel, 2011;Augustin et al., 2017) 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 ratemodel formulation lends itself much more easily to mathematical analysis (Ermentrout and Kopell, 1991;Ermentrout, 1998;Hoppensteadt and Izhikevich, 2012;Ashwin et al., 2016) 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 (Ermentrout and Kopell, 1991;Borisyuk et al., 1995;Ermentrout, 1998), the detail of longrange connectivity matters. Based on observations in several neural areas (Huntley and Jones, 1991;Capaday et al., 2009;Chavane et al., 2011;Huang et al., 2014;Hao et al., 2016), 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 Adesnik and Scanziani, 2010;Levy and Reyes, 2012). Results, such as ours and previous ones (Ermentrout and Kopell, 1991;Borisyuk et al., 1995;Ermentrout, 1998) 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 (Battaglia et al., 2007) 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 (Borisyuk et al., 1995) 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 Veit et al. (2017) which record and perturb distant neural ensembles oscillating in the beta/low gamma range. They may also be relevant for oscillations in large ensembles of neurons coupled by long-distance connections, e.g., inter-areal dynamics (Roberts et al., 2019).
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 finitesize 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., Buice and Chow, 2013;Schwalger et al., 2017;Dumont et al., 2017b) than the one we used. But, at least for our purposes, the simple procedure of Brunel and Hakim (1999) 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 → 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 Cross and Hohenberg, 1993 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 (Kuramoto, 1984) and extended in previous works to neural systems (Ermentrout and Kopell, 1991;Borisyuk et al., 1995;Ermentrout and Kleinfeld, 2001;Hoppensteadt and Izhikevich, 2012;Ashwin et al., 2016).
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 Ermentrout and Kleinfeld, 2001;Wang, 2010 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 (Rubino et al., 2006;Takahashi et al., 2015;Denker et al., 2018;Rule et al., 2018;Takahashi et al., 2011) taking a cortical column typical size of 500 µm and available data on the range of excitatory connections (Huntley and Jones, 1991;Hao et al., 2016;Capaday et al., 2009). Moreover, this would require that E → 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 → E, I case or in the E → 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 (Figures 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 (Barabási and Stanley, 1995). 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 (Figures 7J,K). However, the measured phase wave velocities appear too high to account for the measured velocities (Rubino et al., 2006;Denker et al., 2018;Rule et al., 2018) 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 (Huntley and Jones, 1991). 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 nonexclusive possibility is that the slow conduction velocity along non-myelinated axons (Nauhaus et al., 2009;Chavane et al., 2011) 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 (Denker et al., 2018;Rule et al., 2018). Recorded planar waves also display preferred propagation axes (Rubino et al., 2006) that it will be interesting to link to anisotropies in connectivity or possibly to geometrical properties of the motor and pre-motor 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. Veit et al., 2017) of the different cortical layers and of their connectivity (Adesnik and Scanziani, 2010) to better understand cortical dynamics.

Simulations of Spiking Networks
For our spiking network simulations, we used exponential integrate-and-fire (EIF) neurons that exhibit the following membrane potential dynamics: where τ m = 10 ms is the membrane time constant, E L = −65 mV is the resting potential set by a leak current, and T = 3.5 mV and V T = −59.9 mV are the sharpness of the spike onset and the spike threshold, respectively. These and other parameters used in our simulations are summarized in Table 2. 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, where ξ i is a Gaussian white noise with zero mean and unit variance, ξ i (t)ξ j (t ′ ) = δ ij δ(t − t ′ ), and X ∈ {E, I} depending on whether neuron i is excitatory or inhibitory. Synaptic currents I syn,i are given by where the J ij are the synaptic weights between post-and presynaptic neurons i and j, respectively, and S j (t) = k δ(t − t (j) k ) is the spike train of neuron j with spike times t (j) k . Please note that we refer to I ext,i and I syn,i as current but that we measure them in Volt in Equation (25) since we have divided both sides by the neuron leak conductance g l . In other terms what we refer to as I for notational convenience, is the input current over g l .
If not specified otherwise, we considered an all-to-all connectivity at the level of an individual module, such that J ij ∈ {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 τ 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 ext E , I ext I such that at steady state, excitatory and inhibitory neurons fire at an average rate of r (s) E = 5 Hz and r (s) I = 10 Hz, respectively. We furthermore fixed the noise strengths σ E , σ I by imposing a total noise (external plus synaptic) of σ tot = 10 mV, using σ 2 tot = σ 2 X + J 2 XE N E τ m r (s) E + J 2 XI N I τ m r (s) I , X ∈ {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 τ 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-toall connectivity with lowered synaptic weightsJ 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 → E, X = E, I in the E → 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 · 10 6 ) used in Figure 3, we also replaced the connections across modules with all-to-all connectivity, where we reduced the synaptic weights accordingly toJ XE = f lr J X E.
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 λ such that the N E input synapses are We give explicit formulas for parameters that depend on network sizes N I , N E that are systematically varied among simulations and reported together with the results. Unless specified otherwise, the connection strengths of case A were used. The inverse of the f-I curve, −1 σ (·), can be obtained by interpolation from the tabulated values given in the Supplementary Material. distributed according to their relative contribution p( x) = e −λ| x| /(1+2/(e λ −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 (Goodman and Brette, 2009). While single-module network simulations could be performed on a local machine (Intel Xeon E5 processor @ 3 GHz, 32 GB RAM) in a reasonable amount of time, the largest network simulations (large N, large L) required > 100 GB RAM and were run on the institute's (IBENS) computational cluster, taking up to several days.
Population firing rates were calculated as the total number of spikes within a time bin t divided by the size of the population; we typically used t = 0.1 ms if not explicitly stated otherwise.

Computation of Correlation Functions
The auto-and cross-correlation functions C 1j (τ ) with j = 1, 2, respectively, as shown in Figures 1-4, were calculated as with M the number of discrete time points t i considered, and δr(t i ) = r(t i ) −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, δr E (t) 2 Poisson =r E /(N E t), from the autocorrelation at τ = 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 t = 1 ms by averaging the activity in 1 ms bins to reduce the sampling noise ∝ 1/dt. We then determined the Hilbert phase from the band-pass filtered version of the deviatory firing rate δr = r −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 ∇φ x (t i ) = φ x+1 (t i ) − φ x (t i ) over five neighboring segments, i.e., ∇φ x (t) = 1 5 4 k=0 ∇φ x+k (t i ) = (φ x+5 − φ 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 |∇φ 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 φ x (t i ) = (φ x (t i ), φ x+1 (t i ), . . . , φ x+4 (t i )), we calculated the Euler distance of laterally shifted versions at later time bins t i+j , We subsequently inferred the phase propagation speed v φ (x, t i ) from the optimal shift opt (j; x, t i ) that minimizes the distance between phase profiles after a given time j t according to: For reasons of computational efficiency, we restricted the search for an optimal shift over shifts of ±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 (Ostojic and Brunel, 2011), Ostojic and Brunel (OB) proposed the analytically-motivated choice, In the present paper, the adaptive response time was chosen following Augustin et al. (2017), 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 responseR 1 (f ) of EIF neurons (Richardson, 2007) 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/ 1 + (2πf τ fit ) 2 , which is the Fourier transformation of an exponential decay in the time domain (see Supplementary Figure 1C). We used these values of τ fit as a look-up table and interpolated in between. The resulting "fitted adaptive timescale" (FAT) curve τ (FAT) (I) is shown in Supplementary Figure 1B together with τ (OB) (I). We also used look-up tables and interpolation for σ (I) and the first and second derivatives of σ (I) to speed up the numerical calculations. The tabulated functions σ (I) and τ (FAT) (I) used for the present paper are provided as Supplementary Material.

Simulations of the FAT Rate-Model for an E-I Module
For a network with synaptic weights w EE , w IE , and w EI (see Figure 1 and Table 2 for the values used in the simulations) the currents I ext E , I ext I are chosen so as to impose the steady firing rates r (s) E , r (s) I for the excitatory and inhibitory populations. Namely, where I (s) E and I (s) I are the currents corresponding to the chosen rates according to r (s) E = σ (I (s) E ), r (s) I = σ (I (s) I ). 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 [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. On an Intel Xeon E5 processor @ 3 GHz, a simulation of about 10 s evolution of a network chain runs in about 1 min wall time, which proved sufficient for our purpose.

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 (s) E , r (s) I (i.e., the chosen fixed point is stable) or departures from r (s) E , r (s) I are amplified and the rates cannot be maintained (i.e., the fixed point is unstable). These different cases are shown in parameter space in Figure 1A. The regions of Figure 1A are obtained in a standard way (Wilson and Cowan, 1972;Hoppensteadt and Izhikevich, 2012) by linearizing Equations (3, 4) around the fixed point currents I (s) = (I (s) E , I (s) I ). Using a vector notation for the currents, I = (I E , I I ), the perturbation δI = I − I (s) evolves according to where matrix multiplication is denoted by a dot and L EI is the 2 × 2 linear stability matrix (30) One can note that the derivatives τ ′ E (I (s) E ), τ ′ I (I (s) I ) 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 τ E and τ 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(κt), with the two arguments κ being solutions of Here, we have defined the two positive constants The conditions for the stability of the steady state are that the real parts of the two roots of Equation (31) are negative, namely that the sum of the two roots is negative and their product is positive, respectively 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 α − 1 is not too large, namely in the domain (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 Equations (3, 4) with Different choices of the forcing currents, I forc E (t) and I forc I (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 (Pikovsky et al., 2003). Methods to compute the effect of a weak forcing are well-known (Kuramoto, 1984) and have been applied to neural rate models in numerous works (see Ermentrout and Kleinfeld, 2001;Hoppensteadt and Izhikevich, 2012;Ashwin et al., 2016 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., Equation 29), we adopt a vector notation for the currents. We suppose that the synaptic parameters are such that, with the current I ext 0 alone, the E-I module follows the limit cycle I (0) (t + φ) = (I (0) E (t + φ), I (0) I (t + φ)) of period T [i.e., I (0) (t + T) = I (0) (t)]. Here, the phase φ is an arbitrary constant number reflecting the invariance of the dynamics under time translation when I ext is constant and there is no timedependent forcing. When perturbed by a small time-dependent current I forc (t), the E-I module currents I(t) closely follow the unperturbed limit cycle I (0) (t + φ). However, the phase φ 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π as it is sometimes done.
The matrix L describes the linearized dynamics around the limit cycle, The last term in Equation (39), proportional to the time derivative of the phase φ, comes from the slow phase evolution φ = φ(t) that will be required for preventing the appearance of growing secular terms (Bender and Orszag, 2013), as explained below.
Equation (38) is a linear equation with coefficients which are periodic functions of time. For this Floquet problem (Ince, 1956), it is helpful to consider the evolution of δI(t) from one period to the next. For this purpose, we introduce the so-called monodromy matrix M(t) such that with M(0) = Id, 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 e 1 = (I 0 with eigenvalue µ 1 = 1. The other eigenvector of M(T), e 2 , has an eigenvalue µ 2 with |µ 2 | < 1 since the limit cycle is stable. Equation (38) can be explicitly integrated by writing δI(t) and F(t) in the T-periodic basis obtained by translating e 1 , e 2 along the limit cycle, with e j (t) = M(t) e j , j = 1, 2, for 0 ≤ t < T, with e j (t) = e j (t − T) for t ≥ T. Substitution in Equation (38) gives That is after one turn, where we have used that M(T)e j = µ j e j . When |µ j | < 1, the above recurrence relation implies that x j (n T) has a magnitude comparable to the integral on the r.h.s. of Equation (47). In the simple case when f j (t) is periodic of period T, the integral is constant and x j (n T) tends toward the finite value x * j , However for µ 1 = 1, Equation (47) is an arithmetic series and x 1 grows with time, if the integral term in Equation (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 Equation (39). We can now choose the phase evolution so that the integral term in Equation (47) exactly vanishes and the validity of the perturbative expansion is preserved, requiring that T+t t du f 1 (u) = 0 .
This determines the E-I module phase drift dynamics to be where we have computed the component f 1 (t) of F (Equation 39) with the help of the first vector g 1 (t) = (g 1,E (t), g 1,I (t)) of the bi-orthogonal basis (g 1 (t), g 2 (t)), such that g i (t) · e j (t) = δ ij . Equation (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 (Kuramoto, 1984;Ermentrout and Terman, 2010).

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 Equations (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 Equation (50) with F forc (t) given by (51) Substitution in Equation (50) and averaging over the noise gives a diffusive evolution for the phase φ (Kuramoto, 1984), The finite number of neurons in both the excitatory and inhibitory populations contributes to the diffusion constant D N , as described by Equation (9). The two constants D E and D I in Equation (9) are expressed as integrals along the limit cycle as For the synaptic weights w EE = 1.6 mVs, w EI = 0.32 mVs, w IE = 2.0 mVs, Equations (53, 54) give D E = 1.2 · 10 4 ms, D I = 2.0 · 10 3 ms ( Table 1). Stochastic effects produce a diffusion of the oscillatory phase, the mean square amplitude of which grows linearly in time, as described by Equation (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 (Kuramoto, 1984;Yoshimura and Arai, 2008;Goldobin et al., 2010) of magnitudes 1/N E and 1/N I , changes that we do not consider here.
The autocorrelation decay time τ 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 (0) E (t) is a periodic function of period T. It can be written in a Fourier series as 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 (0) E (t + φ(t)). With the help of Equation (55), it is then easy to compute the auto-correlation of the excitatory activity C EE (t). That is, where the angular brackets indicate averages over different realizations of the noise, i.e., the fluctuating values of the phase φ(t). The auto-correlation decay time τ D is thus obtained as given by Equation (10). Equations (9, 10) show that the oscillation coherence time τ D grows linearly with the number of neurons in the E-I module. For our reference parameters (solid circle in Figure 1A) where T = 64 ms, Equation (10) gives τ D ≃ 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 10 times longer. A fit of the activity auto-correlation decay with the above expression allows one to directly measure τ D in simulations. This procedure is followed in Figure 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 (Equation 50) allows us to derive Equation (13) governing the evolution of the respective phases of two weakly coupled oscillatory modules (Kuramoto, 1984;Ermentrout and Kopell, 1991;Hoppensteadt and Izhikevich, 2012;Ashwin et al., 2016).
We suppose that the parameters are such that for a vanishing coupling f lr (Equations 11,14), the two E-I modules follow the limit cycle I (0) (t) = (I (0) E (t), I (0) I (t)) of period T with different phases φ 1 and φ 2 , In order to determine how a small coupling f lr makes the phases φ 1 and φ 2 evolve, we compute perturbatively in f lr the corrections δI j (t + φ j ), j = 1, 2, to the zeroth-order evolution (57). We obtain that δI 1 follows Equation (38) with the following forcing current (Equation 37) that arises from the coupling to the second E-I module, 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 → E connectivity (Equations 11, 12), the forcing current is given by Equation (50) then determines the phase drift of the first E-I module to be The evolution of φ 2 is obtained by permutation of φ 1 and φ 2 in Equation (60). Subtraction of these two equations provides the evolution of the relative phase φ = φ 1 − φ 2 of the two modules as given by Equation (13), with the explicit expression of the coupling function S E being 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 φ with The constant D 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 → E, I connectivity (Equations 11,14), the coupling function K of Equation (59) is replaced by This gives in turn the coupling function for the phase difference φ between the two modules, For a small phase difference between the two modules, analogously to Equation (61), Equation (65) reduces to Equation (15), the analog of Equation (62) for E → E, I connectivity, with (66) In this case, the constant D 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 Equations (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) + δI j (t), j = 1, 2. The dynamical behavior of the perturbations δI j (t) is obtained by linearizing the dynamics around the fully synchronized oscillatory state, 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 (δI 1 (t), δI 2 (t)) can be reduced to the study of the 2-dimensional evolution of symmetric (δI 1 = δI 2 ) and antisymmetric (δI 1 = −δI 2 ) perturbations. The evolution of symmetric perturbations is governed by the 2 × 2 matrix L S (t) which is identical to the matrix L(t) (Equation 41) for a single E-I module, 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 × 2 matrix L A (t), The matrix L A (t) explicitly 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 L A (t), For E → E connectivity (Equations 11, 12), the matrices L 11 and L 12 read Integration of Equation (70) for different f lr and other parameters fixed, shows that the two eigenvalues of M E A (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 Figure 2B.
For E → E, I connectivity (Equations 11,14), the calculation is identical but the resulting matrix L EI A (t) reads The eigenvalues of the associated monodromy matrix M EI A (T)are found to be smaller than 1 irrespective of the value of f lr (see Figure 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 E φ (Equation 63) and D EI φ (Equation 66), respectively. This can also be seen directly by computing the largest eigenvalue µ A 1 of M A (T) ( Equation 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 (Equation 42), with µ A 1 = µ 1 = 1. For small f lr → 0, µ A 1 can be computed perturbatively. We write the expansion of L A (t) and M A (T) in powers of f lr to first order as where L(t) is given by Equation (41) and the matrix M by Equation (42). The matrix δM A is obtained from the integration of Equation (70) to linear order in f lr , First-order perturbation theory gives the development of the eigenvalue µ A 1 as a function of δM A , where g 1 the left eigenvector of M(T) associated to the eigenvalue 1, such that g 1 · e 1 = 1 [i.e., the first vector of the biorthogonal basis (g 1 , g 2 )]. Explicit formulas are obtained after the replacement of δM A by its expression (75) together with the identity g 1 (u) = g 1 M −1 (u) and the expressions of δL A for the two considered connectivities. For E → E, I connectivity, Equation (73) gives δL EI A , and subsequently the expression for the largest eigenvalue of M EI A (T), The same formula with D EI φ replaced by D E φ is obtained for the largest eigenvalue of M E A (T) with the help of Equations (71) and (72). Since D EI φ > 0 and D E φ < 0, one recovers that E → E, I connectivity synchronizes the two modules' oscillations (µ A 1 < 1) while full synchronization is unstable (µ A 1 > 1) at small coupling for E → E connectivity.

Competition Between Synchronization and Stochasticity for Two Coupled E-I Modules of Finite Size
For E → 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 (Equations 8, 52) and synchronization due to long-range coupling (Equation 13). Thus, one obtains for the phase difference dynamics where in the second approximate equality we have replaced S EI ( φ) by its linear approximation for small phase differences (Equation 15). The noise term η 12 (t) is the difference of the independent finite-size noises for the two modules (52). Being a linear combination of white noises, η 12 (t) is also a white noise, with an amplitude twice as large as for a single module, η 12 (t) = η 1 (t) − η 2 (t), η 12 (t)η 12 (t ′ ) = 2D N δ(t − t ′ ) . (79) Equations (78) and (79) readily give that for a small noise amplitude (D N ≪ f lr D EI φ T), the two-module oscillation phase difference has Gaussian fluctuations with a mean square amplitude given by Equation (16).
We can use that result to estimate the mean square amplitude of the excitatory firing rate difference δr E = (r 1,E − r 2,E )/ √ 2 between the two modules. In the case of weak coupling, the firing rates can be approximated as r j, , j = 1, 2, as long as |φ j | ≪ T. The firing rate fluctuations are then directly obtained from the phase differences as where we used r ′(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 phasemediated, contribution due to Poissonian sampling, which would add an additional term r (0) E (t)/(N E t) with 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 intermodule 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 r 2 E T − r E 2 T . 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 δr E (t) 2 T that takes this periodicity into account. Using the Fourier decomposition Equation (55) and r j,E (t) ≃ r (0) E (t+φ j (t)), j = 1, 2, one obtains analogously to the single module autocorrelation Here, we first used that 1 T T 0 dte i2π (n+m)t/T = δ n,−m and that the relative phases φ 1 (t), φ 2 (t) vary slowly with time with respect to the oscillatory dynamics, i.e., can be considered constant over one oscillation period.
Using δr 2 E T = ( r 2 1,E T + r 2 2,E T )/2 − r 1,E r 2,E T , it is now straightforward to obtain A comparison of this theoretical result with simulations of the stochastic FAT rate model and spiking networks is shown in Figure 5.

Weak Noise at Arbitrary Coupling
Equation (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 Equation (67), I j (t) = I (0) (t) + δI j (t), j = 1, 2, and compute the dynamics of the perturbations δI j (t) resulting from the presence of stochasticity. For the first module, this gives dδI 1 dt = L 11 [I (0) (t)]δI 1 + L 12 [I (0) (t)]δI 2 + B 1 + f lr B 12 , (83) with the stochastic terms Permuting indices 1 and 2 in Equations (83) and (84) gives the dynamics of the linear departure δ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 δI A (t) in the following, From Equation (83) and the analogous equation for module 2, the antisymmetric mode obeys with the matrix L EI A for E → E, I connectivity given by Equation (73). The stochastic term B EI A is given by with (88) In order to solve Equation (86), we introduce again the matrix M EI A (t), as defined by Equation (70). The matrix M EI A (T) has eigenvectors e A 1 , e A 2 with eigenvalues µ A 1 and µ A 2 . The two eigenvalues and eigenvectors are real when the coupling f lr between the two modules is sufficiently small (f lr 0.13 for our reference parameters, see Figure 6G), but they can be complex for larger couplings.
Equation (86) can be explicitly integrated by writing δI A (t) and B A (t) in the T-periodic basis obtained by translating e A 1 , e A 2 along the limit cycle, for 0 ≤ t < T, with e A j (t) = e A j (t − T) for t ≥ T. Analogously to a single oscillatory E-I module subject to a forcing current (Equations 43, 44), we define Substitution of these expressions in the differential equation (86) simply gives One readily obtains, for 0 ≤ t ≤ T, We note that The amplitude of the current fluctuations along the limit cycle is finally obtained, with the help of Equations (90), (93), and (94), as Since e A j (T) = µ A j e A j , 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) = ′ E [I (0) E (t)]δI j,E (t), j = 1, 2, up to the stochastic contribution due to Poisson sampling. The amplitudes of the rate fluctuations δr E (t) 2 = [r 1,E (t) − r 2,E (t)] 2 /2 along the limit cycle, up to the direct Poisson contribution, are thus directly obtained by multiplying the above expression with φ ′ E [I (0) E (t)] 2 . Averaging over one oscillation period, we find (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, Equation (80) or respectively (82) when taking the periodicity of the limit cycle dynamics into account (see Figure 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 Equation (82).
In the weak-coupling regime (f lr → 0), Equation (97) gives back the simple expression (16) obtained from the stochastic phase equation. This can be seen as follows. When f lr → 0, µ A 1 tends toward 1 (Equation 77), its value for a single module. Therefore, the r.h.s. of Equation (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 Equation (77), when f lr → 0 one obtains 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 A j and e A 1 tend toward those of the single module limit cycle, g j and e 1 , when f lr → 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, (99) We used this relation before when deriving Equation (80). Comparison of Equations (98) and (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 I n (t) = I (0) (t), n = 1, · · · , L, is also always an exact solution for the dynamics of a chain of L E-I modules (Equations 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 I n (t) = I (0) (t) + δI n (t), n = 1, · · · , L. The linear evolution of the perturbations δI n is found to be 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 δI n (t) = exp(iqn)δĨ q , n = 1, · · · , L, with L "wavevectors" q = 2πk/L, n = 0, · · · , L − 1. This gives, for each q, the evolution with 2 × 2 matrices L q that we give below. For E → E connectivity, the matrices L nm read, for n = m, and for n = m, The corresponding matrices L E q read The matrix L E q is identical to the previous matrix L E A (Equations 69, 71, 72) governing the stability of two coupled modules with the replacement of (1 − 2f lr ) by the effective couplingC λ (q), 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 forC λ (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 Equation (22) in the main text [for the coupling choice (20)]. For E → E, I connectivity, the matrix L EI q is similarly given by the previous matrix L EI A (Equation 73) with (1 − 2f lr ) replaced byC λ (q). Therefore, in this case, the eigenvalues of L EI q 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 (Equations 13 and 15 or Equation 77), using the correspondence between f lr andC λ (q), where in the second approximate equality we used that q ≪ 1. This provides the relaxation rate of long-wavelength modes given in Equation (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 → E, I connectivity (Equations 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 Equation (17) as an example, we rewrite the long-range synaptic coupling as where we have made use of the expression (5) of the firing rate. The term between large parenthesis on the r.h.s. of Equation (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 I n (t) = (I E,n (t), I I,n (t)) which describe the dynamics of the nth module is provided by I (0) (t + φ n ) = (I (0) E (t + φ n ), I (0) I (t + φ n )), namely the vector of currents on the deterministic limit cycle at a particular phase φ n .
Linearization of the chain dynamics around this approximation gives, similarly to Equations (38, 39), dδI n (t + φ n ) dt = L[I (0) (t + φ n )] · δI n (t + φ n ) + F forc n (t + φ n ) − dI (0) (t + φ n ) dφ n dφ n dt (109) where the matrix L is given by Equation (41). The forcing of the nth module F forc n can be decomposed as a sum of inputs coming from the mean activities, I (0) (t + φ m ), of nearby modules and of fluctuations of their activities, With the help of Equation (108), one obtains In Equation (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 non-trivial function of the phase difference, as in Equation (60) for the two-module case.) In Equation (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 Equation (50) with the constant D EI φ given by Equation (63). The value of D EI φ for the reference parameters is provided in Table 1. The stochastic terms η n in Equation (113)  where the second equality in Equation (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λ) for L ≫ 1/λ]. Equations (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 4.12.1. Long-Wavelength Modes As described by Equation (113), the phases of the different oscillating E-I modules of the spatially extended network are stochastic quantities. Since Equation (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 φ n and the stochastic terms η n in Fourier space, φ n (t) = 1 √ L k=0,··· ,L−1φ q (k) (t) exp(iq (k) n) , η n (t) = 1 √ L k=0,··· ,L−1η q (k) (t) exp(iq (k) n) , where the L wavevectors are {q (k) = 2πk/L, k = 0, · · · , L − 1}. The stochastic termsη q (k) are Gaussian as linear sums of the η n (t) and their correlations are obtained from Equations (114, 115) as where we have made use of the identity n exp(iqn)S(n) =C 2 λ (q).
Note that the deterministic part of Equation (120) exactly corresponds to our previous Equation (107) for the relaxation of the long-wavelength modes. For k = 0, the above equation can easily be integrated to obtainφ q (k) (t), (121) Therefore, the mean square amplitude for the E-I module oscillation phases for a modulation of wavevector q (k) = 0 is given by where we have used the expression (118) for the noise correlation. Equation (122) is the analog for the chain of modules of Equation (16) in the two-module case. It is relatively straightforward to transform the |φ q (k) (t)| 2 into an expression for the amplitudes of the firing rate fluctuations |r q (k) (t)| 2 for k = 0, along the lines of the previous calculations. Namely, we find with [φ n −φ m ] 2 = −2 φ n φ m = −2 L−1 j=0 exp iq (j) (m − n) |φ q (j) (t)| 2 .
(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 r (0) E T /(N E t) due to Poisson sampling of spikes) alongside a more precise estimate valid for all couplings (see below) and FAT rate model results in Figure 7I.
When q → 0, the coupling strength vanishes. Equation (122) shows that the mean square amplitudes of the long-wavelength modes then diverge as q −2 , |φ q (t)| 2 ∼ D I N I + D E N E cosh(λ) − 1 D EI φ q 2 , q → 0 .
In this limit, the phase equation (121) reduces to the classical Edwards-Wilkinson equation (Barabási and Stanley, 1995), where D N = D I /N I +D E /N E is the local module noise amplitude (Equation 9) and the long-wavelength diffusion constant D lw is given by

Fluctuations of Arbitrary Wavelengths
Equations (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 → 0, one eigenvalue of the matrix L EI q is very close to 1 and the corresponding (phase) mode dominates the fluctuation of activity at the corresponding wavelength 2π/q. For modulations of smaller wavelengths, the modules are more strongly coupled. The two modes of L EI q significantly contribute to the fluctuation of activity and Equation (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 (Equations 83-96).
The stochastic terms δB q (k) are the Fourier components of the B n (Equation 129), N EC λ (q (k) )ξ E,q (k) (t) w IE (134) Equations (133) and (134) generalize to arbitrary Fourier modes the previous Equations (86) and (87) for the antisymmetric mode in the two-module case with the with the already-noted replacement of (1 − 2f lr ) by the Fourier transformC λ (q) (Equation 106). The amplitudes of the Fourier mode modulations thus read, similarly to Equations (95,96), where the vectors e q j (t) are defined for the matrices L EI q (t) as the vectors e A j (t) for the matrix L EI A (t). The amplitudes of the rate fluctuations are again obtained by multiplying the expression above by ′ E [I (0) E (t)] 2 , and averaged over the complete limit cycle we find Here, we explicitly added the contribution that arises due to Poissonian sampling of the spiking dynamics within bins of size T. The above result is compared to the stochastic FAT rate model in Figure 7I.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
AK, JR, and VH performed the research and wrote the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
AK acknowledges the support of Fondation de la Recherche Médicale (FRM) through a doctoral FRM fellowship.