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 (> 10⁵), 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⁴) 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.

1. 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 10Hz and 20Hz was organized in plane waves and spiral waves in the turtle cortex. This spiral-like organization was also reported for pharmacologically induced 10Hz 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⁻¹ (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 spatially-structured 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 spiking-network 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 rate-model formalism. A variety of dynamical regimes was found and qualitative agreement with the bifurcations observed in simulations of a spiking network was reported when the neurons in each E–I module oscillated in a well-synchronized manner. Visually induced gamma oscillations and their dependence on visual stimulus contrast (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 and Hakim, 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 distance-dependent 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⁴ − 10⁵ 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 rate-model framework by adding to it a stochastic component. Even when different modules would fully synchronize in a deterministic context, stochasticity in the module activities produces differences in the phases of oscillation of neighboring modules. (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 weak-noise reduction to a phase dynamics is accurate only for modules with a very large number of neurons. We obtain a quantitative description of stochastic effects for modules of moderate size, by computing the noise influence on both the phase and amplitude of the module oscillation at the linear level.

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 two-module 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.

2. Results

2.1. 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 rate-model formulation (Ostojic and Brunel, 2011; Augustin et al., 2017) to represent the activity of a group of EIF neurons:

$\begin{array}{ll}\tau (I)\frac{d}{dt}I=-I+I_0+s(t),& (1)\end{array}$

$\begin{array}{ll}r(t)={\Phi }_{\sigma }[I(t)],& (2)\end{array}$

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 network-averaged 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),

$\begin{array}{ll}{\tau }_E(I_E)\frac{d{I}_E}{dt}=-I_E+I_E^{ext}+w_{EE}{r}_E-w_{EI}{r}_I,& (3)\end{array}$

$\begin{array}{ll}{\tau }_I(I_I)\frac{d{I}_I}{dt}=-I_I+I_I^{ext}+w_{IE}{r}_E,& (4)\end{array}$

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), ${\tau }_E(I)={\tau }_I(I)={\tau }^{(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_I^{ext}$.

FIGURE 1

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_E^{(s)}=5\text{Hz}$, $r_I^{(s)}=10\text{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\text{m}{\text{V}}^2{\text{s}}^2$; B—w_EE = 1.6 mVs, $w_{IE}{w}_{EI}=1.28\text{m}{\text{V}}^2{\text{s}}^2$; C—w_EE = 1.76 mVs, $w_{IE}{w}_{EI}=1.28\text{m}{\text{V}}^2{\text{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⁶, 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⁴) 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⁴, 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).

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_IEw_EI, where w_IE is the total excitatory synaptic strength on inhibitory neurons and w_EI the total inhibitory synaptic strength on excitatory neurons. When inhibitory feedback is weak, the steady state is subject to a non-oscillatory instability: recurrent excitation leads to steady firing at a very high rate, basically limited by the neuron refractory period in our simple model (other mechanisms, not considered here, such as SFA or pair-pulse synaptic depression can moderate this regime). When inhibitory feedback is sufficiently strong, the steady state is destabilized by an oscillatory instability. This instability can lead to finite amplitude oscillations but also to steady high frequency discharge when a steady high-rate fixed point exists (see 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⁶) 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⁴ − 10⁵). 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),

$\begin{array}{ll}{r}_E(t)={\Phi }_E[I_E(t)]+\sqrt{{\Phi }_E[I_E(t)]/N_E}{\xi }_E(t),& (5)\end{array}$

$\begin{array}{ll}{r}_I(t)={\Phi }_E[I_I(t)]+\sqrt{{\Phi }_I[I_I(t)]/N_I}{\xi }_I(t),& (6)\end{array}$

where ξ_E(t), ξ_I(t) are independent white noises ($\langle {\xi }_E(t){\xi }_E(t^{\prime })\rangle =\langle {\xi }_I(t){\xi }_I(t^{\prime })\rangle =\delta (t-t^{\prime })$, $\langle {\xi }_E(t){\xi }_I(t^{\prime })\rangle =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^{\prime })=\langle r_E(t)r_E(t^{\prime })\rangle -{\langle r_E\rangle }^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,

$\begin{array}{ll}|C_{EE}(t)|~\text{exp}(-t/{\tau }_D).& (7)\end{array}$

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,

$\begin{array}{ll}\langle {[\varphi (t)-\varphi (0)]}^2\rangle =D_{N}t.& (8)\end{array}$

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,

$\begin{array}{ll}{D}_N=\frac{D_E}{N_E}+\frac{D_I}{N_I}.& (9)\end{array}$

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,

$\begin{array}{ll}{\tau }_D=\frac{1}{2{\pi }^2}\frac{T^2}{D_N}.& (10)\end{array}$

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.

2.2. Dynamical Regimes of Two Oscillatory E–I Modules Coupled by Long-Range Excitation

2.2.1. 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.

FIGURE 2

Figure 2. Synchronization and dynamical regimes for two E–I modules coupled by long-range excitation. (A) Sketch of the coupled E–I modules and their synaptic interactions. (Top) In the E → E connectivity case, long-range excitation only targets excitatory neurons in other modules. (Bottom) In the E → E, I connectivity case, long-range excitation targets excitatory and inhibitory neurons in distant modules in equal proportion. (B) Different dynamical behaviors for the coupled modules as a function of the coupling strength f_lr and the synaptic excitatory strength on inhibitory neurons w_IE. The labeled points (solid black circles) mark the parameters of the dynamical regimes shown in (B1–6), the cross (light gray) corresponds to the parameter used in Figure 3A. (B1) Fully synchronized state for sufficiently strong long-range excitation; f_lr = 0.03. (B2) Finite phase difference: the activity of one module is greater than the other, the amplitudes of their oscillatory activities are constant in time but the phases of their oscillations are different; f_lr = 0.025. (B3) Modulated dominance: one of the two modules is more active than the other but the amplitudes of the oscillatory activities of both modules vary themselves periodically in time; f_lr = 0.02. (B4) Alternating dominance: the activities of the two E–I modules successively dominate; f_lr = 0.015. (B5) Antiphase regime at very low coupling and strong local excitation on inhibitory neurons; f_lr = 0.001. (B6) Finite phase difference regime at very low coupling and weak local excitation on inhibitory neurons; f_lr = 0.001 and w_IE = 1.44 mVs. (C) The synchronization function governing the evolution of the relative phase of two weakly coupled modules in the E → E connectivity case. Modules with weak local excitation of inhibitory neurons w_IE (blue solid line), strong local excitation of inhibitory neurons (orange dashed-dotted line) and w_IE close to the threshold strength of local excitation separating the two regimes (red dashed line). (D) The synchronization function governing the evolution of the relative phase of two weakly coupled modules in the E → 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) (solid black line).

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

$\begin{array}{l}{\tau }_E(I_{E,1})\frac{d{I}_{E,1}}{dt}=-I_{E,1}+I_E^{ext}+w_{EE}[(1-f_{lr})r_{E,1}+f_{lr}{r}_{E,2}]\\ -w_{EI}{r}_{I,1}& (11)\end{array}$

$\begin{array}{ll}{\tau }_I(I_{I,1})\frac{d{I}_{I,1}}{dt}=-I_{I,1}+I_I^{ext}+w_{IE}{r}_{E,1}& (12)\end{array}$

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_IEw_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):

$\begin{array}{ll}\frac{d\Delta \varphi }{dt}=f_{lr}{S}_E(\Delta \varphi ).& (13)\end{array}$

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^{\prime }(\Delta \varphi )>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^{\prime }(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 non-trivial 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 finite-phase difference phase at very low f_lr. A first decrease of f_lr below that necessary for full synchrony leads to a phase where the two-modules oscillate at the same frequency and different phases. In contrast to the finite-phase difference at very low f_lr, oscillation amplitude in one module strongly dominates the other one, as illustrated in 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 1Hz range, as shown in Figure 2B3. This regime stems from a Hopf bifurcation of the finite-difference regime, with the unusual low frequency coming from the small value of f_lr (see Methods: Stability Analysis of Full Synchronization for Two Coupled E–I Modules). A further decrease of f_lr transforms the modulated-dominance regime into an “alternating dominance” regime: oscillations in one module dominate the oscillations in the other one, but the dominant module switches periodically, again with a slow frequency in the 1Hz 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.

2.2.2. 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

$\begin{array}{ll}{\tau }_I(I_{I,1})\frac{d{I}_{I,1}}{dt}=-I_{I,1}+I_I^{ext}+w_{IE}[(1-f_r)r_{E,1}+f_r{r}_{E,2}]& (14)\end{array}$

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

$\begin{array}{ll}{S}_{EI}(\Delta \varphi )\simeq -2\Delta \varphi D_{\varphi }^{EI},& (15)\end{array}$

with the constant $D_{\varphi }^{EI}$ depending on the characteristics of the oscillation cycle (see Equation 66 in Methods). The constant $D_{\varphi }^{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).

TABLE 1

Table 1. Computed constants.

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.

2.2.3. 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 rate-model results and networks with a large number of neurons (N ~ 10⁶) for E → E connectivity. The different dynamical regimes for two E–I modules coupled by long-range 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 connectivity, two large E–I modules oscillate in full synchrony as predicted by the rate description (not shown).

FIGURE 3

Figure 3. Two large E–I modules coupled by long-range excitation. (1,2) E–I modules with a large number of neurons, N = 1.6 · 10⁶, and E → E connectivity. The different dynamical regimes depicted in Figure 2 for different coupling strengths are clearly visible both in (1) spiking network simulations and in (2) the corresponding rate-model simulations. The respective auto- and crosscorrelations of the excitatory activities are shown in panels (3). The coupling strengths are as follows: (A) f_lr = 0.01 (light gray cross in Figure 2B), (B) f_lr = 0.015, (C) f_lr = 0.02, (D) f_lr = 0.05.

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.

FIGURE 4

Figure 4. Stochastic effects due to network size for two E–I modules coupled by long-range excitation. Auto- and crosscorrelations C₁₁ and C₁₂ of the activities of E–I modules with a finite number of neurons N = 10⁴ and 2 · 10⁵ for different coupling strengths f_lr = 0.01 and 0.02 with long-range excitation targeting only excitatory populations (A–D) or targeting both excitatory and inhibitory neurons (G–J) (parameter values are stated in the panels' lower left corner). Summary plots with the ratio of equal-time cross-correlations over autocorrelations for varying N and f_lr are shown in (E,F,K,L), respectively. In the E → E case, oscillatory synchrony between the two modules disappears when the coupling strength f_lr diminishes [(C) vs. (D)], but the dynamical regimes of Figure 2 are blurred by stochastic effects when the network size N decreases [(A,B) vs. (C,D)]. 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.

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⁵ neurons, whereas it is not apparent for N = 10⁴ 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⁵ 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 Figures 4G–J. Both effects are quantified in Figures 4K,L.

One should note that whereas the precise targeting of the long-range excitation plays a crucial role in the instability of in-phase oscillations for large networks (N ~ 10⁶), it is much less important for smaller (N ~ 10⁴) 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).

2.2.4. 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,

$\begin{array}{ll}\langle {({\varphi }_1-{\varphi }_2)}^2\rangle =\frac{D_N}{2f_{lr}{D}_{\varphi }^{EI}}.& (16)\end{array}$

In this expression, the constant D_N is simply the amplitude of stochastic fluctuations in a single module (Equation 9), whereas $f_{lr}{D}_{\varphi }^{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.

FIGURE 5

Figure 5. Competition between synchronization and noise for two E–I modules with E → E, I connectivity: theory vs. simulations. The variance of excitatory activity between the two modules $\langle \delta r_E^2\rangle =\langle {(r_{1,E}-r_{2,E})}^2\rangle /2$ is shown as a function of coupling strength f_lr. (A) Simulations of the stochastic rate model (FAT, duration 100 s, averaged over three repetitions) for different noise strengths (color symbols), as specified in the insert in term of the number of neurons N in the individual E–I module. These simulation results are plotted against theoretical predictions with different approximations. For very small coupling strengths, the phase approximation with the full limit cycle dynamics (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 98, 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. (B) Same as in (A) with the results for the network simulations of Figures 4K,L (colored solid lines corresponding to different E–I module sizes, as sprecified in the inset), shown as $\langle \delta r_E^2\rangle =C_{11}(0)-C_{12}(0)$.

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 ${\mu }_1^A$ is close to 1, so that the restoring strength $1-{\mu }_1^A$ of the corresponding (phase) mode is much smaller than $1-{\mu }_2^A$, the restoring strength of the other mode. However, an explicit computation shows that ${\mu }_1^A$ quickly decreases as f_lr increases (see Equation 76 in Methods). Therefore, the two modes soon play comparable roles as f_lr increases.

FIGURE 6

Figure 6. Spontaneous appearance of phase gradients in a chain of E–I modules. (A) Sketch of the chain of E–I modules coupled by long-range excitation. Excitatory neurons in one module target excitatory neurons (E → E) or excitatory and inhibitory neurons (E → E, I) in the other modules with a distance-dependent excitation profile. (B) Fourier transform of the coupling function ${\tilde{C}}_{\lambda }(q)$ (Equation (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;\lambda )=(1-{\tilde{C}}_{\lambda }(q))/2$ for perturbations of wavenumber q. (D) Stability spectrum for perturbations of wavelength q as function of the effective coupling f_lr(q; λ) 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 $\sum _{n=1}^L{(r_{n,E}-\sum _{m=1}^L{r}_{m,E}/L)}^2$ over time, with the predicted rate for the fastest-growing mode shown on top (dashed black line). (G,H) Same as (D,E) but for E → 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.

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 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, $\langle {(r_{1,E}-r_{2,E})}^2\rangle$, even for modules with only 10⁴ neurons when stochasticity is quite high. This second approach is thus quite successful in quantitatively capturing the observed competition between synchronization and stochasticity.

2.3. Phase Gradients in a Chain of Oscillatory E–I Modules Coupled by Long-Range Excitation

Having analyzed the synchronization dynamics and effects of finite-size fluctuations for two E–I modules, we now turn to a fuller description of spatially-structured connectivity. We consider a chain of E–I modules, with each module's long-range excitation targeting other modules in a distance-dependent way. Interestingly, the effects that we uncovered in the two-module setting, such as desynchronization at weak coupling for E → 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.

2.3.1. 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

$\begin{array}{l}{\tau }_E(I_{E,n})\frac{d{I}_{E,n}}{dt}=-I_{E,n}+I_E^{ext}+w_{EE}{\displaystyle \sum _{m=1,\cdots ,L}}C(n,m)r_{E,m}\\ -w_{EI}{r}_{I,n}& (17)\end{array}$

$\begin{array}{ll}{\tau }_I(I_{I,n})\frac{d{I}_{I,n}}{dt}=-I_{I,n}+I_I^{ext}+w_{IE}{r}_{E,n}& (18)\end{array}$

with C(n, m) = C(|n − m|).

For E → E, I connectivity, Equation (18) is replaced by

$\begin{array}{ll}{\tau }_I(I_{I,n})\frac{d{I}_{I,n}}{dt}=-I_{I,n}+I_I^{ext}+w_{IE}{\displaystyle \sum _{m=1,\cdots ,L}C}(n,m)r_{E,m}.& (19)\end{array}$

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

$\begin{array}{l}C(l)=A_{\lambda }\left[\text{exp}(-\lambda |l|)+\text{exp}(-\lambda (L-|l|))\right],\\ A_{\lambda }=\frac{1-\text{exp}(-\lambda )}{[1-\text{exp}(-\lambda L)][1+\text{exp}(-\lambda )]}& (20)\end{array}$

where A_λ implements the normalization $\sum _{l=0,\cdots ,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 E–I module network considered in the previous section, provided the coupling in the two-module network is chosen as

$\begin{array}{ll}{f}_{lr}(q)=[1-{\tilde{C}}_{\lambda }(q)]/2,& (21)\end{array}$

where ${\tilde{C}}_{\lambda }(q)$ is the Fourier transform of the long-range excitation profile C(l) (Equation 106 in Methods). ${\tilde{C}}_{\lambda }(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) requires mathematical analysis (see Methods: Stability Analysis of Full Synchronization for a Chain of Oscillatory E–I Modules), it can be intuitively understood. For perturbations of small wavenumbers (q → 0), only distant parts of the chain are in distinct oscillating states. The relative dynamics of these distant parts is thus effectively only weakly coupled by excitation since long-range excitation decreases with distance. The relation (21) allows one to directly transcribe for the E–I module chain the results previously obtained for the two E–I module case. As before, we compare two connectivities for long-range excitation.

We first consider the case of E → 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. In-phase synchronization is unstable when f_lr becomes smaller than a critical $f_{lr}^{*}$ that depends on the single E–I module synaptic parameters. For the chain of modules, this implies that spatially-periodic perturbations grow when their wavenumber q lies below a critical q^*, defined by

$\begin{array}{ll}\text{cos}(q^{*})=\frac{1-2f_{lr}^{*}\text{cosh}(\lambda )}{1-2f_{lr}^{*}},\text{ or }{q}^{*}\simeq \sqrt{\frac{2f_{lr}^{*}}{1-2f_{lr}^{*}}}\lambda ,& (22)\end{array}$

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

$\begin{array}{ll}{\sigma }_q=\frac{D_{\varphi }^{EI}}{2[\text{cosh}(\lambda )-1]}{q}^2.& (23)\end{array}$

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.

2.3.2. 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⁵.

FIGURE 7

Figure 7. Appearance of phase gradients and stochastic fluctuations in a chain of E–I modules. (A,B) Simulation of the stochastic rate model (FAT) for a chain with E → E connectivity, a space constant λ = 0.33 and network size N = 10⁴. (A) Instantaneous phase as obtained from the Hilbert transform of bandpass-filtered excitatory activity, where the (unfiltered) rates are shown in (B). (C) Distributions of local phase gradients ∇ϕ = ϕ_x+1 − ϕ_x for different noise strengths and space constants for E → E connectivity. Inset: distributions of gradients normalized to the respective mean value. (D,E) Phase and excitatory activity for the same parameters as in (A,B) but with E → E, I connectivity. (F) Same as in (C) but for E → E, I connectivity. (G,H) Power spectra of the spatial activity profiles (averaged over 5 s) for different module sizes and space constants λ, for E → E (G) and E → E, I (H) activity, respectively. The spectra obtained from spiking network simulations (N = 5, 000: L = 128; N = 10⁴: L = 32) are shown in light thick lines, spectra obtained from rate-model simulations are shown in dark thin lines (L = 256). (I) Comparison of the spectra obtained from simulations (FAT, solid lines) and theoretical predictions for E → 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 (138)), 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 124). 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 $\overline{\nabla \varphi }$ (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).

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 two-module 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⁴), 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.

2.3.3. Phase Waves and Propagation of Activity

We have shown that stochastic fluctuations in a chain of E–I modules create phase differences in the oscillations of E–I modules along the chain. In the presence of phase differences, activity in different modules peaks at different times and may produce “phase waves” (also described as traveling waves, e.g., in Wang, 2010, p1219), 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),

$\begin{array}{ll}{v}_{\varphi }=\frac{d}{\nabla \varphi },& (24)\end{array}$

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⁴ 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³ 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.

3. 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 rate-model 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 long-range 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 finite-size fluctuations and synchronization in order to describe the dynamics at the scale of a cortical column. Other studies have proposed more refined and complex ways to take into account finite-size fluctuations (see e.g., 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 non-exclusive 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.

4. Methods

4.1. Simulations and Numerical Analysis

4.1.1. Simulations of Spiking Networks

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

$\begin{array}{ll}{\tau }_m\frac{d}{dt}{V}_i=E_L-V_i+{\Delta }_T{e}^{\frac{V_i-V_T}{{\Delta }_T}}+I_{ext,i}+I_{syn,i},& (25)\end{array}$

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,

TABLE 2

Table 2. Summary of parameters used.

$\begin{array}{ll}{I}_{ext,i}=I_X^{ext}+{\sigma }_X\sqrt{{\tau }_m}{\xi }_i,& (26)\end{array}$

where ξ_i is a Gaussian white noise with zero mean and unit variance, $\langle {\xi }_i(t){\xi }_j(t^{\prime })\rangle ={\delta }_{ij}\delta (t-t^{\prime })$, and X ∈ {E, I} depending on whether neuron i is excitatory or inhibitory. Synaptic currents I_syn,i are given by

$\begin{array}{ll}{I}_{syn,i}={\tau }_m{\displaystyle \sum _j{J}_{ij}}{S}_j(t),& (27)\end{array}$

where the J_ij are the synaptic weights between post- and pre-synaptic neurons i and j, respectively, and $S_j(t)=\sum _k\delta (t-t_k^{(j)})$ is the spike train of neuron j with spike times $t_k^{(j)}$. 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_E^{ext},I_I^{ext}$ such that at steady state, excitatory and inhibitory neurons fire at an average rate of $r_E^{(s)}=5\text{Hz}$ and $r_I^{(s)}=10\text{Hz}$, respectively. We furthermore fixed the noise strengths σ_E, σ_I by imposing a total noise (external plus synaptic) of σ_tot = 10 mV, using ${\sigma }_{tot}^2={\sigma }_X^2+J_{XE}^2{N}_E{\tau }_m{r}_E^{(s)}+J_{XI}^2{N}_I{\tau }_m{r}_I^{(s)}$, X ∈ {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-to-all connectivity with lowered synaptic weights ${\tilde{J}}_{XE}=(1-f_{lr})J_{XE}$ for those synaptic connections within a given module that are supplemented by long-range excitation (X = E in the E → 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_lrN_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⁶) used in Figure 3, we also replaced the connections across modules with all-to-all connectivity, where we reduced the synaptic weights accordingly to ${\tilde{J}}_{XE}=f_{lr}{J}_{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 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 @ 3GHz, 32GB RAM) in a reasonable amount of time, the largest network simulations (large N, large L) required > 100GB 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.

4.1.2. 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

$\begin{array}{ll}{C}_{1j}(\tau )=\frac{1}{M}{\displaystyle \sum _{t_i}\delta }{r}_{1,E}(t_i)\delta r_{j,E}(t_i+\tau ),& (28)\end{array}$

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

4.1.3. 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 $\delta r=r-\bar{r}$, using a fifth-order Butterworth filter with lower and upper corner frequencies of 15 and 40 Hz, respectively. The Hilbert transform was obtained using the Python scipy package.

In a next step, we determined locally averaged phase gradients by averaging phase differences ∇ϕ_x(t_i) = ϕ_x+1(t_i) − ϕ_x(t_i) over five neighboring segments, i.e., $\overline{\nabla {\varphi }_x}(t)=\frac{1}{5}\sum _{k=0}^4\nabla {\varphi }_{x+k}(t_i)=({\varphi }_{x+5}-{\varphi }_x)/5$ (note that all phase differences are taken modulo a full oscillation period). For the detection of propagation events, we retained gradients where the absolute average phase gradient $|\overline{\nabla {\varphi }_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,

$\begin{array}{l}d(\Delta ,j;x,t_i)=\Vert {\varphi }_x(t_i)-{\varphi }_{x+\Delta }(t_{i+j}){\Vert }^2.\end{array}$

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:

$\begin{array}{l}{v}_{\varphi }(x,t_i)=\frac{1}{5}{\displaystyle \sum _{j=1}^5\frac{{\Delta }_{opt}(j;x,t_i)}{j\Delta t}}.\end{array}$

For reasons of computational efficiency, we restricted the search for an optimal shift Δ over shifts of ±15 modules.

4.2. 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,

$\begin{array}{l}{\tau }^{(OB)}(I)={\tau }_m{\Delta }_T\frac{{\Phi }_{\sigma }^{\prime }(I)}{{\Phi }_{\sigma }(I)}.\end{array}$

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 response ${\widehat{R}}_1(f)$ of EIF neurons (Richardson, 2007) for frequencies from 1Hz to 1 kHz in steps of 1Hz for all values of the input current on the grid of input currents I with step size dI = 0.1 mV, and fitted the modulus of the firing-rate response by $A/\sqrt{1+{(2\pi f{\tau }_{fit})}^2}$, which is the Fourier transformation of an exponential decay in the time domain (see 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.

4.3. 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_E^{ext},I_I^{ext}$ are chosen so as to impose the steady firing rates $r_E^{(s)},r_I^{(s)}$ for the excitatory and inhibitory populations. Namely,

$\begin{array}{l}{I}_E^{ext}=I_E^{(s)}-w_{EE}{r}_E^{(s)}+w_{EI}{r}_I^{(s)},\\ I_I^{ext}=I_I^{(s)}-w_{IE}{r}_E^{(s)},\end{array}$

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

We implemented the stochastic rate-model equations in Python using an Euler scheme with a time step dt = 0.01 ms. At each time t_i, we calculated the noisy firing rate r_X(t_i) by drawing a number of emitted spikes n_X(t_i) from a Poisson distribution with mean N_XΦ[I_X(t_i)]dt, and setting r_X(t_i) = n_X(t_i)/(N_Xdt), X ∈ E, I, in order to avoid negative firing rates. On an Intel Xeon E5 processor @ 3GHz, a simulation of about 10 s evolution of a network chain runs in about 1 min wall time, which proved sufficient for our purpose.

4.4. Oscillatory Instability for an E–I Module

The E–I module (3,4) either has its E and I populations stably firing at the chosen rates $r_E^{(s)},r_I^{(s)}$ (i.e., the chosen fixed point is stable) or departures from $r_E^{(s)},r_I^{(s)}$ are amplified and the rates cannot be maintained (i.e., the fixed point is unstable). These different cases are shown in parameter space in 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 ${\text{I}}^{(s)}=(I_E^{(s)},I_I^{(s)})$. Using a vector notation for the currents, I = (I_E, I_I), the perturbation δI = I − I^(s) evolves according to

$\begin{array}{ll}\frac{d}{dt}\delta \text{I}=L_{EI}·\delta \text{I},& (29)\end{array}$

where matrix multiplication is denoted by a dot and L_EI is the 2 × 2 linear stability matrix

$\begin{array}{rr}\hfill L_{EI}=\left(\begin{array}{cc}\hfill [-1+w_{EE}{\Phi }_E^{\prime }(I_E^{(s)})]/{\tau }_E(I_E^{(s)})& \hfill -w_{EI}{\Phi }_I^{\prime }(I_I^{(s)})/{\tau }_E(I_E^{(s)})\\ \hfill w_{IE}{\Phi }_E^{\prime }(I_E^{(s)})/{\tau }_I(I_I^{(s)})& \hfill -1/{\tau }_I(I_I^{(s)})\end{array}\right).& (30)\end{array}$

One can note that the derivatives ${\tau }_E^{\prime }(I_E^{(s)}),{\tau }_I^{\prime }(I_I^{(s)})$ do not appear at the linear level since these timescales multiply time derivatives which are themselves small quantities at the linear level. Thus, the adaptive character of τ_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

$\begin{array}{ll}{\tau }_E{\tau }_I{\kappa }^2+[{\tau }_E+{\tau }_I(1-\alpha )]\kappa +1-\alpha +\beta =0.& (31)\end{array}$

Here, we have defined the two positive constants

$\begin{array}{ll}\alpha ={\Phi }_E^{\prime }(I_E^{(s)})w_{EE}\text{ and }\beta ={\Phi }_I^{\prime }(I_I^{(s)}){\Phi }_E^{\prime }(I_E^{(s)})w_{EI}{w}_{IE}.& (32)\end{array}$

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

$\begin{array}{ll}\alpha <1+\frac{{\tau }_E}{{\tau }_I},& (33)\end{array}$

$\begin{array}{ll}\beta >\alpha -1.& (34)\end{array}$

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

$\begin{array}{ll}\alpha >1+\frac{{\tau }_E}{{\tau }_I},& (35)\end{array}$

$\begin{array}{ll}\beta >\frac{{\tau }_I}{4{\tau }_E}{\left(\alpha -1\right)}^2+\frac{1}{2}(\alpha -1)+\frac{{\tau }_E}{4{\tau }_I}=\frac{{\tau }_I}{4{\tau }_E}{\left(\alpha -1+\frac{{\tau }_E}{{\tau }_I}\right)}^2.& (36)\end{array}$

4.5. 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

$\begin{array}{ll}{I}_E^{ext}(t)=I_{E,0}^{ext}+I_E^{forc}(t),I_I^{ext}(t)=I_{I,0}^{ext}+I_I^{forc}(t).& (37)\end{array}$

Different choices of the forcing currents, $I_E^{forc}(t)$ and $I_I^{forc}(t)$, allow us to describe the influence of finite-size stochastic effects as well as the synchronization dynamics of coupled E–I modules. The computation of the drift and entrainment of oscillators is a classical topic (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 ${\text{I}}_0^{ext}$ alone, the E–I module follows the limit cycle ${\text{I}}^{(0)}(t+\varphi )=(I_E^{(0)}(t+\varphi ),I_I^{(0)}(t+\varphi ))$ of period T [i.e., I⁽⁰⁾(t + T) = I⁽⁰⁾(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 time-dependent 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⁽⁰⁾(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.

In order to determine how I^forc(t) entrains the phase ϕ, we compute perturbatively the correction δI(t + ϕ) to the zeroth-order evolution I⁽⁰⁾(t + ϕ). From Equations (3) and (4), we obtain for δI

$\begin{array}{ll}\frac{d\delta \text{I}(\zeta )}{d\zeta }={L}[{\text{I}}^{(0)}(\zeta )]·\delta \text{I}(\zeta )+\text{F}(\zeta )& (38)\end{array}$

with ζ = t + ϕ and

$\begin{array}{ll}\text{F}(t+\varphi )={\text{F}}^{forc}(t+\varphi )-\frac{d{\text{I}}^{(0)}(\zeta )}{d\zeta }\frac{d\varphi }{dt},& (39)\end{array}$

$\begin{array}{ll}{\text{F}}^{forc}(t)=\left(\frac{I_E^{forc}}{{\tau }_E[I_E^{(0)}(t)]},\frac{I_I^{forc}}{{\tau }_I[I_I^{(0)}(t)]}\right).& (40)\end{array}$

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

$\begin{array}{rr}\hfill {L}[{\text{I}}^{(0)}(t)]=\left(\begin{array}{cc}\hfill \frac{-1+w_{EE}{\Phi }_E^{\prime }[I_E^{(0)}(t)]-{\tau }_E^{\prime }[I_E^{(0)}(t)]{I^{\prime }}_E^{(0)}(t)}{{\tau }_E[I_E^{(0)}(t)]}& \hfill -\frac{w_{EI}{\Phi }_I^{\prime }[I_I^{(0)}(t)]}{{\tau }_E[I_E^{(0)}(t)]}\\ \hfill \frac{w_{IE}{\Phi }_E^{\prime }[I_E^{(0)}(t)]}{{\tau }_I[I_I^{(0)}(t)]}& \hfill -\frac{1+{\tau }_I^{\prime }[I_I^{(0)}(t)]{I^{\prime }}_I^{(0)}(t)}{{\tau }_I[I_I^{(0)}(t)]}\end{array}\right).& (41)\end{array}$

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

$\begin{array}{ll}\frac{d}{dt}M={L}[{\text{I}}^{(0)}(t)]·M& (42)\end{array}$

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 ${\text{e}}_1=({I_E^0}^{\prime }(0),{I_I^0}^{\prime }(0))$ is an eigenvector of M(T) with eigenvalue μ₁ = 1. The other eigenvector of M(T), e₂, has an eigenvalue μ₂ with |μ₂| < 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₁, e₂ along the limit cycle,

$\begin{array}{ll}\delta \text{I}(t)=x_1(t){\text{e}}_1(t)+x_2(t){\text{e}}_2(t),& (43)\end{array}$

$\begin{array}{ll}\text{F}(t)=f_1(t){\text{e}}_1(t)+f_2(t){\text{e}}_2(t)& (44)\end{array}$

with

$\begin{array}{ll}{\text{e}}_j(t)=M(t){\text{e}}_j,j=1,2,& (45)\end{array}$

for 0 ≤ t < T, with e_j(t) = e_j(t − T) for t ≥ T. Substitution in Equation (38) gives

$\begin{array}{ll}\frac{d{x}_j}{dt}=f_j(t),j=1,2.& (46)\end{array}$

That is after one turn,

$\begin{array}{ll}{x}_j((n+1)T)={\mu }_j\left[x_j(nT)+{\displaystyle {\int }_{nT}^{(n+1)T}d}t{f}_j(t)\right],& (47)\end{array}$

where we have used that M(T)e_j = μ_je_j. When |μ_j| < 1, the above recurrence relation implies that x_j(nT) 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(nT) tends toward the finite value $x_j^{*}$,

$\begin{array}{ll}{x}_j^{*}=\frac{{\mu }_j}{1-{\mu }_j}{\displaystyle {\int }_0^{T}d}t{f}_j(t).& (48)\end{array}$

However for μ₁ = 1, Equation (47) is an arithmetic series and x₁ 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

$\begin{array}{ll}{\displaystyle {\int }_t^{T+t}}du{f}_1(u)=0.& (49)\end{array}$

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

$\begin{array}{ll}\frac{d\varphi }{dt}=\frac{1}{T}{\displaystyle {\int }_t^{T+t}}du{\text{g}}_1(u)·{\text{F}}_1^{forc}(u),& (50)\end{array}$

where we have computed the component f₁(t) of F (Equation 39) with the help of the first vector g₁(t) = (g_{1, E}(t), g_{1, I}(t)) of the bi-orthogonal basis (g₁(t), g₂(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).

4.6. 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

$\begin{array}{ll}{\text{F}}^{forc}(t)=\left(\begin{array}{c}\frac{1}{{\tau }_E[I_E(t)]}\left[w_{EE}\sqrt{\frac{{\Phi }_E[I_E(t)]}{N_E}}{\xi }_E(t)-w_{EI}\sqrt{\frac{{\Phi }_I[I_I(t)]}{N_I}}{\xi }_I(t)\right]\\ \frac{1}{{\tau }_I[I_I(t)]}{w}_{IE}\sqrt{\frac{{\Phi }_E[I_E(t)]}{N_E}}{\xi }_E(t)\end{array}\right).& (51)\end{array}$

Substitution in Equation (50) and averaging over the noise gives a diffusive evolution for the phase ϕ (Kuramoto, 1984),

$\begin{array}{ll}\frac{d\varphi }{dt}=\eta (t),\langle \eta (t)\eta (t^{\prime })\rangle =D_N\delta (t-t^{\prime }).& (52)\end{array}$

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

$\begin{array}{ll}{D}_E=\frac{1}{T}{\displaystyle {\int }_0^T}dt{\left(\frac{g_{1,E}(t)w_{EE}}{{\tau }_E[I_E^{(0}(t)]}+\frac{g_{1,I}(t)w_{IE}}{{\tau }_I[I_I^{(0)}(t)]}\right)}^2{\Phi }_E[I_E^{(0)}(t)],& (53)\end{array}$

$\begin{array}{ll}{D}_I=\frac{1}{T}{\displaystyle {\int }_0^T}dt{\left(\frac{g_{1,E}(t)w_{EI}}{{\tau }_E[I_E^{(0)}(t)]}\right)}^2{\Phi }_I[I_I^{(0)}(t)].& (54)\end{array}$

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\text{ms}$, $D_I=2.0·10^3\text{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_E^{(0)}(t)$ is a periodic function of period T. It can be written in a Fourier series as

$\begin{array}{ll}{r}_E^{(0)}(t)={\displaystyle \sum _{n=-\infty }^{+\infty }{\tilde{r}}_n}\text{exp}\left(i\frac{2\pi nt}{T}\right).& (55)\end{array}$

For a finite but large number of neurons, the activity of the excitatory population can approximately be described by taking the dominant effect of phase diffusion into account, namely approximating r_E(t) by $r_E^{(0)}(t+\varphi (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,

$\begin{array}{l}{C}_{EE}(t)\simeq \langle r_E^{(0)}(s+\varphi (s))r_E^{(0)}(s+t+\varphi (s+t))\rangle -{\langle r_E^{(0)}(s)\rangle }^2\\ ={\displaystyle \sum _{n=-\infty }^{+\infty }}{\tilde{r}}_n{\tilde{r}}_{-n}\text{exp}\left(i\frac{2\pi nt}{T}\right)\text{exp}\left(-\frac{2{\pi }^2{n}^2}{T^2}\langle {[\varphi (s)-\varphi (s+t)]}^2\rangle \right)\\ -{\tilde{r}}_0^2=2{\displaystyle \sum _{n=1}^{+\infty }|{\tilde{r}}_n{|}^2}\text{cos}\left(\frac{2\pi nt}{T}\right)\text{exp}\left(-\frac{2{\pi }^2{n}^2}{T^2}{D}_{N}t\right),& (56)\end{array}$

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⁴ neurons (D_N = 2.5 ms). For a module with ten times more neurons (N = 10⁵), 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.

4.7. 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 ${\text{I}}^{(0)}(t)=(I_E^{(0)}(t),I_I^{(0)}(t))$ of period T with different phases ϕ₁ and ϕ₂,

$\begin{array}{ll}{\text{I}}_j(t)={\text{I}}^{(0)}(t+{\varphi }_j),j=1,2.& (57)\end{array}$

In order to determine how a small coupling f_lr makes the phases ϕ₁ and ϕ₂ 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₁ follows Equation (38) with the following forcing current (Equation 37) that arises from the coupling to the second E–I module,

$\begin{array}{ll}{\text{F}}^{forc}(t+{\varphi }_1)=f_{lr}\text{K}[{\text{I}}^{(0)}(t+{\varphi }_1),{\text{I}}^{(0)}(t+{\varphi }_2)].& (58)\end{array}$

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

$\begin{array}{l}\text{K}[{\text{I}}^{(0)}(t),{\text{I}}^{(0)}(t+{\varphi }_2-{\varphi }_1)]\\ =\left(\begin{array}{c}\frac{w_{EE}}{{\tau }_E[I_E^{(0)}(t)]}\left({\Phi }_E[I_E^{(0)}(t+{\varphi }_2-{\varphi }_1)]-{\Phi }_E[I_E^{(0)}(t)]\right)\\ 0\end{array}\right).& (59)\end{array}$

Equation (50) then determines the phase drift of the first E–I module to be

$\begin{array}{l}\frac{d{\varphi }_1}{dt}=\frac{f_{lr}}{T}{\displaystyle {\int }_0^T}dt{\text{g}}_1(t)·\text{K}[{\text{I}}^{(0)}(t),{\text{I}}^{(0)}(t+{\varphi }_2-{\varphi }_1)]\\ =\frac{f_{lr}{w}_{EE}}{T}{\displaystyle {\int }_0^T}dt\frac{g_{1,E}(t)}{{\tau }_E[I_E^{(0)}(t)]}({\Phi }_E[I_E^{(0)}(t+{\varphi }_2-{\varphi }_1)]\\ -{\Phi }_E[I_E^{(0)}(t)]).& (60)\end{array}$

The evolution of ϕ₂ is obtained by permutation of ϕ₁ and ϕ₂ in Equation (60). Subtraction of these two equations provides the evolution of the relative phase Δϕ = ϕ₁ − ϕ₂ of the two modules as given by Equation (13), with the explicit expression of the coupling function S_E being

$\begin{array}{l}{S}_E(\Delta \varphi )=\frac{w_{EE}}{T}{\displaystyle {\int }_0^T}dt\frac{g_{1,E}(t)}{{\tau }_E[I_E^{(0)}(t)]}({\Phi }_E[I_E^{(0)}(t-\Delta \varphi )]\\ -{\Phi }_E[I_E^{(0)}(t+\Delta \varphi )]).& (61)\end{array}$

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 Δϕ

$\begin{array}{ll}{S}_E(\Delta \varphi )\simeq -2\Delta \varphi D_{\varphi }^E,& (62)\end{array}$

with

$\begin{array}{ll}{D}_{\varphi }^E=\frac{w_{EE}}{T}{\displaystyle {\int }_0^T}dt\frac{g_{1,E}(t)}{{\tau }_E[I_E^{(0)}(t)]}{\Phi }_E^{\prime }[I_E^{(0)}(t)]\frac{d}{dt}{I}_E^{(0)}(t).& (63)\end{array}$

The constant $D_{\varphi }^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

$\begin{array}{l}\text{K}[{\text{I}}^{(0)}(t),{\text{I}}^{(0)}(t+{\varphi }_2-{\varphi }_1)]\\ =\left(\begin{array}{c}\frac{w_{EE}}{{\tau }_E[I_E^{(0)}(t)]}\left({\Phi }_E[I_E^{(0)}(t+{\varphi }_2-{\varphi }_1)]-{\Phi }_E[I_E^{(0)}(t)]\right)\\ \frac{w_{IE}}{{\tau }_I[I_I^{(0)}(t)]}\left({\Phi }_E[I_E^{(0)}(t+{\varphi }_2-{\varphi }_1)]-{\Phi }_E[I_E^{(0)}(t)]\right)\end{array}\right)& (64)\end{array}$

This gives in turn the coupling function for the phase difference Δϕ between the two modules,

$\begin{array}{l}{S}_{EI}(\Delta \varphi )\\ ={\displaystyle {\int }_0^T}\frac{dt}{T}\left(\frac{g_{1,E}(t)w_{EE}}{{\tau }_E[I_E^{(0)}(t)]}+\frac{g_{1,I}(t)w_{IE}}{{\tau }_I[I_I^{(0)}(t)]}\right)({\Phi }_E[I_E^{(0)}(t-\Delta \varphi )]\\ -{\Phi }_E[I_E^{(0)}(t+\Delta \varphi )]).& (65)\end{array}$

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

$\begin{array}{ll}{D}_{\varphi }^{EI}=\frac{1}{T}{\displaystyle {\int }_0^T}dt\left(\frac{g_{1,E}(t)w_{EE}}{{\tau }_E[I_E^{(0)}(t)]}+\frac{g_{1,I}(t)w_{IE}}{{\tau }_I[I_I^{(0)}(t)]}\right){\Phi }_E^{\prime }[I_E^{(0)}(t)]\frac{d}{dt}{I}_E^{(0)}(t).& (66)\end{array}$

In this case, the constant $D_{\varphi }^{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.

4.8. 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, ${\text{I}}_1(t)={\text{I}}_2(t)={\text{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 ${\text{I}}_j(t)={\text{I}}^{(0)}(t)+{\delta \text{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,

$\begin{array}{ll}\frac{d\delta {\text{I}}_1}{dt}={{L}}_{11}[{\text{I}}^{(0)}(t)]\delta {\text{I}}_1+{{L}}_{12}[{\text{I}}^{(0)}(t)]\delta {\text{I}}_2,& (67)\end{array}$

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₁(t), δI₂(t)) can be reduced to the study of the 2-dimensional evolution of symmetric (δI₁ = δI₂) and antisymmetric (δI₁ = −δI₂) 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,

$\begin{array}{ll}{{L}}_S(t)={{L}}_{11}+{{L}}_{12}={L}.& (68)\end{array}$

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)$,

$\begin{array}{ll}{{L}}_A[{\text{I}}^{(0)}(t)]={{L}}_{11}-{{L}}_{12}.& (69)\end{array}$

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)$,

$\begin{array}{ll}\frac{d}{dt}{M}_A={{L}}_A[{\text{I}}^{(0)}(t)]·M,M_A(0)=Id.& (70)\end{array}$

For E → E connectivity (Equations 11, 12), the matrices ${{L}}_{11}$ and ${{L}}_{12}$ read

$\begin{array}{l}{{L}}_{11}^E[{\text{I}}^{(0)}(t)]\\ =\left(\begin{array}{cc}\frac{-1+w_{EE}(1-f_{lr}){\Phi }_E^{\prime }[I_E^{(0)}(t)]-{\tau }_E^{\prime }[I_E^{(0)}(t)]{I^{\prime }}_E^{(0)}(t)}{{\tau }_E[I_E^{(0)}(t)]}& -\frac{w_{EI}{\Phi }_I^{\prime }[I_I^{(0)}(t)]}{{\tau }_E[I_E^{(0)}(t)]}\\ \frac{w_{IE}{\Phi }_E^{\prime }[I_E^{(0)}(t)]}{{\tau }_I[I_I^{(0)}(t)]}& -\frac{1+{\tau }_I^{\prime }[I_I^{(0)}(t)]{I^{\prime }}_I^{(0)}(t)}{{\tau }_I[I_I^{(0)}(t)]}\end{array}\right),& (71)\end{array}$

$\begin{array}{ll}{{L}}_{12}^E[{\text{I}}^{(0)}(t)]=\left(\begin{array}{cc}\frac{f_{lr}{w}_{EE}{\Phi }_E^{\prime }[I_E^{(0)}(t)]}{{\tau }_E[I_E^{(0)}(t)]}& 0\\ 0& 0\end{array}\right).& (72)\end{array}$

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

For E → E, I connectivity (Equations 11, 14), the calculation is identical but the resulting matrix ${{L}}_A^{EI}(t)$ reads

$\begin{array}{l}{{L}}_A^{EI}[{\text{I}}^{(0)}(t)]\\ =\left(\begin{array}{cc}\frac{-1+w_{EE}(1-2f_{lr}){\Phi }_E^{\prime }[I_E^{(0)}(t)]-{\tau }_E^{\prime }[I_E^{(0)}(t)]{I^{\prime }}_E^{(0)}(t)}{{\tau }_E[I_E^{(0)}(t)]}& -\frac{w_{EI}{\Phi }_I^{\prime }[I_I^{(0)}(t)]}{{\tau }_E[I_E^{(0)}(t)]}\\ \frac{(1-2f_{lr})w_{IE}{\Phi }_E^{\prime }[I_E^{(0)}(t)]}{{\tau }_I[I_I^{(0)}(t)]}& -\frac{1+{\tau }_I^{\prime }[I_I^{(0)}(t)]{I^{\prime }}_I^{(0)}(t)}{{\tau }_I[I_I^{(0)}(t)]}\end{array}\right).& (73)\end{array}$

The eigenvalues of the associated monodromy matrix $M_A^{EI}(T)$are found to be smaller than 1 irrespective of the value of f_lr (see 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_{\varphi }^E$ (Equation 63) and $D_{\varphi }^{EI}$ (Equation 66), respectively. This can also be seen directly by computing the largest eigenvalue ${\mu }_1^A$ 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 ${\mu }_1^A={\mu }_1=1$. For small f_lr → 0, ${\mu }_1^A$ 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

$\begin{array}{ll}{{L}}_A(t)={L}(t)+f_{lr}\delta {{L}}_A(t),M_A(T)=M(T)+f_{lr}\delta M_A,& (74)\end{array}$

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,

$\begin{array}{ll}\delta M_A=M(T){\displaystyle {\int }_0^T}du{M}^{-1}(u)\delta {{L}}_A(u)M(u).& (75)\end{array}$

First-order perturbation theory gives the development of the eigenvalue ${\mu }_1^A$ as a function of δM_A,

$\begin{array}{ll}{\mu }_1^A=1+f_{lr}{\text{g}}_1·\delta M_A{\text{e}}_1,& (76)\end{array}$

where g₁ the left eigenvector of M(T) associated to the eigenvalue 1, such that g₁ · e₁ = 1 [i.e., the first vector of the bi-orthogonal basis (g₁, g₂)]. Explicit formulas are obtained after the replacement of δM_A by its expression (75) together with the identity ${\text{g}}_1(u)={\text{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}}_A^{EI}$, and subsequently the expression for the largest eigenvalue of $M_A^{EI}(T)$,

$\begin{array}{ll}{\mu }_1^A-1=-2f_{lr}{D}_{\varphi }^{EI}T+{O}(f_{lr}^2).& (77)\end{array}$

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

4.9. 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.

4.9.1. 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

$\begin{array}{ll}\frac{d({\varphi }_1-{\varphi }_2)}{dt}=f_{lr}{S}_{EI}({\varphi }_1-{\varphi }_2)+{\eta }_{12}(t)\simeq -2f_{lr}{D}_{\varphi }^{EI}({\varphi }_1-{\varphi }_2)+{\eta }_{12}(t),& (78)\end{array}$

where in the second approximate equality we have replaced S_EI(Δϕ) by its linear approximation for small phase differences (Equation 15). The noise term η₁₂(t) is the difference of the independent finite-size noises for the two modules (52). Being a linear combination of white noises, η₁₂(t) is also a white noise, with an amplitude twice as large as for a single module,

$\begin{array}{ll}{\eta }_{12}(t)={\eta }_1(t)-{\eta }_2(t),\langle {\eta }_{12}(t){\eta }_{12}(t^{\prime })\rangle =2D_N\delta (t-t^{\prime }).& (79)\end{array}$

Equations (78) and (79) readily give that for a small noise amplitude ($D_N\ll f_{lr}{D}_{\varphi }^{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 $\delta r_E=(r_{1,E}-r_{2,E})/\sqrt{2}$ between the two modules. In the case of weak coupling, the firing rates can be approximated as $r_{j,E}(t)\simeq r_E^{(0)}(t)+{r^{\prime }}_E^{(0)}(t){\varphi }_j(t)$, j = 1, 2, as long as |ϕ_j| ≪ T. The firing rate fluctuations are then directly obtained from the phase differences as

$\begin{array}{l}{\langle \delta r_E{(t)}^2\rangle }_T=\frac{1}{2}\langle {\left({\varphi }_1-{\varphi }_2\right)}^2\rangle {\displaystyle {\int }_0^T\frac{dt}{T}}{\left[{\Phi }_E^{\prime }[I_E^{(0)}(t)]\frac{d{I}_E^{(0)}}{dt}(t)\right]}^2\\ =\frac{D_N}{4f_{lr}{D}_{\varphi }^{EI}}{\displaystyle {\int }_0^T\frac{dt}{T}}{\left[{\Phi }_E^{\prime }[I_E^{(0)}(t)]\frac{d{I}_E^{(0)}}{dt}(t)\right]}^2,& (80)\end{array}$

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

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

$\begin{array}{l}{\langle r_{1,E}(t)r_{2,E}(t)\rangle }_T\\ \simeq {\displaystyle {\int }_0^T\frac{dt}{T}}{\displaystyle \sum _{n,m}{\tilde{r}}_n}{\tilde{r}}_m\text{exp}\left(i\frac{2\pi (n+m)t}{T}\right)\langle \text{exp}\left(i\frac{2\pi [n{\varphi }_1(t)+m{\varphi }_2(t)]}{T}\right)\rangle \\ ={\displaystyle \sum _n|}{\tilde{r}}_n{|}^2\langle \text{exp}\left(i\frac{2\pi n}{T}[{\varphi }_1(t)-{\varphi }_2(t)]\right)\rangle \\ ={\displaystyle \sum _n|}{\tilde{r}}_n{|}^2\text{exp}\left(-\frac{2{\pi }^2{n}^2}{T}\langle {[{\varphi }_1(t)-{\varphi }_2(t)]}^2\rangle \right)& (81)\end{array}$

Here, we first used that $\frac{1}{T}{\int }_0^T\text{d}t{e}^{i2\pi (n+m)t/T}={\delta }_{n,-m}$ and that the relative phases ϕ₁(t), ϕ₂(t) vary slowly with time with respect to the oscillatory dynamics, i.e., can be considered constant over one oscillation period.

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

$\begin{array}{ll}{\langle \delta r_E{(t)}^2\rangle }_T=2{\displaystyle \sum _{n=1}^{+\infty }|}{\tilde{r}}_n{|}^2\left[1-\text{exp}\left(-\frac{2{\pi }^2{n}^2}{T}\langle {({\varphi }_1-{\varphi }_2)}^2\rangle \right)\right].& (82)\end{array}$

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

4.9.2. 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), ${\text{I}}_j(t)={\text{I}}^{(0)}(t)+{\delta \text{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

$\begin{array}{ll}\frac{d\delta {\text{I}}_1}{dt}={{L}}_{11}[{\text{I}}^{(0)}(t)]\delta {\text{I}}_1+{{L}}_{12}[{\text{I}}^{(0)}(t)]\delta {\text{I}}_2+{\text{B}}_1+f_{lr}{\text{B}}_{12},& (83)\end{array}$

with the stochastic terms

$\begin{array}{l}{\text{B}}_1=\left(\begin{array}{r}\hfill \frac{w_{EE}}{{\tau }_E[I_E^{(0)}(t)]}\\ \hfill \frac{w_{IE}}{{\tau }_I[I_I^{(0)}(t)]}\end{array}\right)(1-f_{lr})\sqrt{\frac{{\Phi }_E[I_E^{(0)}(t)]}{N_E}}{\xi }_{E,1}(t)\\ -\left(\begin{array}{c}\frac{w_{EI}}{{\tau }_E[I_E^{(0)}(t)]}\\ 0\end{array}\right)\sqrt{\frac{{\Phi }_I[I_I^{(0)}(t)]}{N_I}}{\xi }_{I,1}(t),\\ {\text{B}}_{12}=\left(\begin{array}{r}\hfill \frac{w_{EE}}{{\tau }_E[I_E^{(0)}(t)]}\\ \hfill \frac{w_{IE}}{{\tau }_I[I_I^{(0)}(t)]}\end{array}\right)\sqrt{\frac{{\Phi }_E[I_E^{(0)}(t)]}{N_E}}{\xi }_{E,2}(t)& (84)\end{array}$

Permuting indices 1 and 2 in Equations (83) and (84) gives the dynamics of the linear departure δI₂ 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,

$\begin{array}{ll}\delta {\text{I}}_A=\frac{1}{\sqrt{2}}(\delta {\text{I}}_1-\delta {\text{I}}_2).& (85)\end{array}$

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

$\begin{array}{ll}\frac{d\delta {\text{I}}_A}{dt}={{L}}_A^{EI}[{\text{I}}^{(0)}(t)]\delta {\text{I}}_A+{\text{B}}_A^{EI},& (86)\end{array}$

with the matrix ${{L}}_A^{EI}$ for E → E, I connectivity given by Equation (73). The stochastic term ${\text{B}}_A^{EI}$ is given by

$\begin{array}{l}{\text{B}}_A^{EI}=\left(\begin{array}{r}\hfill \frac{w_{EE}}{{\tau }_E[I_E^{(0)}(t)]}\\ \hfill \frac{w_{IE}}{{\tau }_I[I_I^{(0)}(t)]}\end{array}\right)(1-2f_{lr})\sqrt{\frac{{\Phi }_E[I_E^{(0)}(t)]}{N_E}}{\xi }_{E,A}(t)\\ -\left(\begin{array}{c}\frac{w_{EI}}{{\tau }_E[I_E^{(0)}(t)]}\\ 0\end{array}\right)\sqrt{\frac{{\Phi }_I[I_I^{(0)}(t)]}{N_I}}{\xi }_{I,A}(t)& (87)\end{array}$

with

$\begin{array}{ll}{\xi }_{E,A}(t)=\frac{1}{\sqrt{2}}[{\xi }_{E,1}(t)-{\xi }_{E,2}(t)],{\xi }_{I,A}(t)=\frac{1}{\sqrt{2}}[{\xi }_{I,1}(t)-{\xi }_{I,2}(t)].& (88)\end{array}$

In order to solve Equation (86), we introduce again the matrix $M_A^{EI}(t)$, as defined by Equation (70). The matrix $M_A^{EI}(T)$ has eigenvectors ${\text{e}}_1^A,{\text{e}}_2^A$ with eigenvalues ${\mu }_1^A$ and ${\mu }_2^A$. The two eigenvalues and eigenvectors are real when the coupling f_lr between the two modules is sufficiently small (f_lr ≲ 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 ${\text{e}}_1^A,{\text{e}}_2^A$ along the limit cycle,

$\begin{array}{ll}{\text{e}}_j^A(t)=M_A(t){\text{e}}_j^A,j=1,2,& (89)\end{array}$

for 0 ≤ t < T, with ${\text{e}}_j^A(t)={\text{e}}_j^A(t-T)$ for t ≥ T. Analogously to a single oscillatory E–I module subject to a forcing current (Equations 43, 44), we define

$\begin{array}{ll}\delta {\text{I}}_A(t)=x_1(t){\text{e}}_1^A(t)+x_2(t){\text{e}}_2^A(t),& (90)\end{array}$

$\begin{array}{ll}{\text{B}}_A(t)=b_1(t){\text{e}}_1^A(t)+b_2(t){\text{e}}_2^A(t).& (91)\end{array}$

Substitution of these expressions in the differential equation (86) simply gives

$\begin{array}{ll}\frac{d{x}_j}{dt}=b_j(t),j=1,2.& (92)\end{array}$

One readily obtains, for 0 ≤ t ≤ T,

$\begin{array}{l}{x}_j(t)={\displaystyle {\int }_{-\infty }^t}du{b}_j(u)={\displaystyle {\int }_0^t}du{\text{g}}_j^A(u)·{\text{B}}_A(u)\\ \begin{array}{l}+{\displaystyle \sum _{n=1}^{+\infty }}{\mu }_j^n{\displaystyle {\int }_0^t}du{\text{g}}_j^A(u)·{\text{B}}_A(u-nT).\end{array}& (93)\end{array}$

We note that

$\begin{array}{ll}\langle {\text{g}}_j^A(u)·{\text{B}}_A(u){\text{g}}_k^A(v)·{\text{B}}_A(v)\rangle =b_{jk}\delta (u-v),& (94)\end{array}$

with

$\begin{array}{l}{b}_{jk}={(1-2f_{lr})}^2\left(\frac{w_{EE}{g}_{j,E}^A(u)}{{\tau }_E[I_E^{(0)}(u)]}+\frac{w_{IE}{g}_{j,I}^A(u)}{{\tau }_I[I_I^{(0)}(u)]}\right)(\frac{w_{EE}{g}_{k,E}^A(u)}{{\tau }_E[I_E^{(0)}(u)]}\\ +\frac{w_{IE}{g}_{k,I}^A(u)}{{\tau }_I[I_I^{(0)}(u)]})\frac{{\Phi }_E[I_E^{(0)}(u)]}{N_E}\\ +\frac{w_{EI}^2{g}_{j,E}^A(u)g_{k,E}^A(u)}{{\tau }_E^2[I_E^{(0)}(u)]}\frac{{\Phi }_I[I_I^{(0)}(t)]}{N_I}& (95)\end{array}$

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

$\begin{array}{l}\langle \delta I_{A,E}{(t)}^2\rangle \\ ={\displaystyle \sum _{j,k=1,2}}{e}_{j,E}^A(t)e_{k,E}^A(t)\left[{\displaystyle {\int }_0^t}du{b}_{jk}(u)+\frac{{\mu }_j{\mu }_k}{1-{\mu }_j{\mu }_k}{\displaystyle {\int }_0^T}du{b}_{jk}(u)\right].& (96)\end{array}$

Since ${\text{e}}_j^A(T)={\mu }_j^A{\text{e}}_j^A,j=1,2$, the obtained expression is periodic in time, as it should.

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

$\begin{array}{l}{\langle \delta r_E{(t)}^2\rangle }_T=\frac{1}{T}{\displaystyle {\int }_0^T}dt{\displaystyle \sum _{j,k=1,2}}{\varphi }_E^{\prime }{\left[I_E^{(0)}(t)\right]}^2{e}_{j,E}^A(t)e_{k,E}^A(t)\\ \begin{array}{l}\times \left[{\displaystyle {\int }_0^t}du{b}_{jk}(u)+\frac{{\mu }_j{\mu }_k}{1-{\mu }_j{\mu }_k}{\displaystyle {\int }_0^T}du{b}_{jk}(u)\right].\end{array}& (97)\end{array}$

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, ${\mu }_1^A$ 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

$\begin{array}{l}{\langle \delta r_E{(t)}^2\rangle }_T\simeq \frac{1}{T}{\displaystyle {\int }_0^T}dt{\varphi }_E^{\prime }{\left[I_E^{(0)}(t)\right]}^2\frac{{(e_{1,E}^A(t))}^2}{4f_{lr}{D}_{\varphi }^{EI}T}{\displaystyle {\int }_0^T}du{b}_{11}(u)\\ \begin{array}{l}\simeq \frac{1}{T}{\displaystyle {\int }_0^T}dt{\varphi }_E^{\prime }{\left[I_E^{(0)}(t)\right]}^2\frac{{(e_{1,E}(t))}^2}{4f_{lr}{D}_{\varphi }^{EI}T}\left(\frac{D_I}{N_I}+{(1-2f_{lr})}^2\frac{D_E}{N_E}\right),\end{array}& (98)\end{array}$

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 ${\text{g}}_j^A$ and ${\text{e}}_1^A$ tend toward those of the single module limit cycle, g_j and e₁, 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,

$\begin{array}{ll}{\langle \delta r_E{(t)}^2\rangle }_T=\frac{1}{2}\langle {({\varphi }_1-{\varphi }_2)}^2\rangle {\displaystyle {\int }_0^T}\frac{dt}{T}{\left[{\varphi }_E^{\prime }\left[I_E^{(0)}(t)\right]\frac{d{I}_E^{(0)}}{dt}(t)\right]}^2.& (99)\end{array}$

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₁(t) is simply the velocity along the current limit cycle.

4.10. Stability Analysis of Full Synchronization for a Chain of Oscillatory E–I Modules

The fully synchronized oscillatory state ${\text{I}}_n(t)={\text{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 ${\text{I}}_n(t)={\text{I}}^{(0)}(t)+{\delta \text{I}}_n(t)$, n = 1, ⋯, L. The linear evolution of the perturbations δI_n is found to be

$\begin{array}{ll}\frac{d{\delta \text{I}}_n}{dt}={\displaystyle \sum _{m=1,\cdots ,L}}{{L}}_{nm}\left[{\text{I}}^{(0)}(t)\right]{\delta \text{I}}_m.& (100)\end{array}$

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

$\begin{array}{ll}\delta {\text{I}}_n(t)=\text{exp}(iqn)\delta {\tilde{\text{I}}}_q,n=1,\cdots ,L,& (101)\end{array}$

with L “wavevectors” q = 2πk/L, n = 0, ⋯, L − 1. This gives, for each q, the evolution

$\begin{array}{ll}\frac{d\delta {\tilde{\text{I}}}_q}{dt}={{L}}_q\left[{\text{I}}^{(0)}(t)\right]\delta {\tilde{\text{I}}}_q& (102)\end{array}$

with 2 × 2 matrices ${{L}}_q$ that we give below.

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

$\begin{array}{l}{{L}}_{nn}^E[{\text{I}}^{(0)}(t)]\\ =\left(\begin{array}{cc}\frac{-1+w_{EE}C(0){\Phi }_E^{\prime }[I_E^{(0)}(t)]-{\tau }_E^{\prime }[I_E^{(0)}(t)]{I^{\prime }}_E^{(0)}(t)}{{\tau }_E[I_E^{(0)}(t)]}& -\frac{w_{EI}{\Phi }_I^{\prime }[I_I^{(0)}(t)]}{{\tau }_E[I_E^{(0)}(t)]}\\ \frac{w_{IE}{\Phi }_E^{\prime }[I_E^{(0)}(t)]}{{\tau }_I[I_I^{(0)}(t)]}& -\frac{1+{{\tau }^{\prime }}_I[I_I^{(0)}(t)]{I^{\prime }}_I^{(0)}(t)}{{\tau }_I[I_I^{(0)}(t)]}\end{array}\right),& (103)\end{array}$

and for n ≠ m,

$\begin{array}{ll}{{L}}_{nm}^E[{\text{I}}^{(0)}(t)]=\left(\begin{array}{cc}\frac{w_{EE}C(n-m){\Phi }_E^{\prime }[I_E^{(0)}(t)]}{{\tau }_E[I_E^{(0)}(t)]}& 0\\ 0& 0\end{array}\right).& (104)\end{array}$

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

$\begin{array}{l}{{L}}_q^E[{\text{I}}^{(0)}(t)]\\ =\left(\begin{array}{cc}\frac{-1+w_{EE}{\tilde{C}}_{\lambda }(q){\Phi }_E^{\prime }[I_E^{(0)}(t)]-{\tau }_E^{\prime }[I_E^{(0)}(t)]{I^{\prime }}_E^{(0)}(t)}{{\tau }_E[I_E^{(0)}(t)]}& -\frac{w_{EI}{\Phi }_I^{\prime }[I_I^{(0)}(t)]}{{\tau }_E[I_E^{(0)}(t)]}\\ \frac{w_{IE}{\Phi }_E^{\prime }[I_E^{(0)}(t)]}{{\tau }_I[I_I^{(0)}(t)]}& -\frac{1+{{\tau }^{\prime }}_I[I_I^{(0)}(t)]{I^{\prime }}_I^{(0)}(t)}{{\tau }_I[I_I^{(0)}(t)]}\end{array}\right)& (105)\end{array}$

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

$\begin{array}{ll}{\tilde{C}}_{\lambda }(q)={\displaystyle \sum _{l=0,\cdots ,L-1}}\text{exp}(iql)C(l)=\frac{\text{cosh}(\lambda )-1}{\text{cosh}(\lambda )-cos(q)},& (106)\end{array}$

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

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

$\begin{array}{ll}\frac{d\delta {\tilde{\text{I}}}_q}{dt}=D_{\varphi }^{EI}\left[{\tilde{C}}_{\lambda }(q)-1\right]\delta {\tilde{\text{I}}}_q(t)\simeq -\frac{D_{\varphi }^{EI}{q}^2}{2\left[\text{cosh}(\lambda )-1\right]}\delta {\tilde{\text{I}}}_q(t),& (107)\end{array}$

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.

4.11. 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

$\begin{array}{l}{\displaystyle \sum _{m=1,\cdots ,L}}C(n,m)r_{E,m}={\Phi }_E[I_{E,n}(t)]+\\ {\displaystyle \sum _{m=1,\cdots ,L}}C(n,m)(\left\{{\Phi }_E[I_{E,m}(t)]-{\Phi }_E[I_{E,n}(t)]\right\}\\ \begin{array}{l}+\sqrt{\frac{{\Phi }_E[I_{E,m}(t)]}{N_E}}{\xi }_{E,m}(t)),\end{array}& (108)\end{array}$

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 ${\text{I}}^{(0)}(t+{\varphi }_n)=(I_E^{(0)}(t+{\varphi }_n),I_I^{(0)}(t+{\varphi }_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),

$\begin{array}{l}\frac{d\delta {\text{I}}_n(t+{\varphi }_n)}{dt}={L}[{\text{I}}^{(0)}(t+{\varphi }_n)]·\delta {\text{I}}_n(t+{\varphi }_n)+{\text{F}}_n^{forc}(t+{\varphi }_n)\\ -\frac{d{\text{I}}^{(0)}(t+{\varphi }_n)}{d{\varphi }_n}\frac{d{\varphi }_n}{dt}& (109)\end{array}$

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

$\begin{array}{ll}{\text{F}}_n^{forc}={\text{F}}_n^{f,mean}+{\text{F}}_n^{f,stoch}.& (110)\end{array}$

With the help of Equation (108), one obtains

$\begin{array}{l}{\text{F}}_n^{f,mean}=\left(\begin{array}{l}\frac{w_{EE}}{{\tau }_E[I_E^{(0)}(t+{\varphi }_n)]}\hfill \\ \frac{w_{IE}}{{\tau }_I[I_I^{(0)}(t+{\varphi }_n)]}\hfill \end{array}\right){\displaystyle \sum _{m=1}^{L}C}(n,m)\{{\Phi }_E[I_E^{(0}(t+{\varphi }_m)]\\ -{\Phi }_E[I_E^0(t+{\varphi }_n)]\}\\ \simeq \left(\begin{array}{l}\frac{w_{EE}}{{\tau }_E[I_E^{(0)}(t+{\varphi }_n)]}\hfill \\ \frac{w_{IE}}{{\tau }_I[I_I^{(0)}(t+{\varphi }_n)]}\hfill \end{array}\right){\displaystyle \sum _{m=1}^{L}C}(n,m)({\varphi }_m-{\varphi }_n){\Phi }_E^{\prime }[I_E^{(0)}(t\\ +{\varphi }_n)]\frac{d{I}_E^{(0)}(t+{\varphi }_n)}{dt}& (111)\end{array}$

and

$\begin{array}{l}{\text{F}}_n^{f,stoch}=\left(\begin{array}{c}\frac{w_{EE}}{{\tau }_E\left[I_E^{(0)}(t+{\varphi }_n)\right]}\\ \frac{w_{IE}}{{\tau }_I\left[I_I^{(0)}(t+{\varphi }_n)\right]}\end{array}\right){\displaystyle \sum _{m=1}^L}C(n,m)\sqrt{\frac{{\Phi }_E\left[I_E^{(0)}(t+{\varphi }_m)\right]}{N_E}}{\xi }_{E,m}(t)\\ \begin{array}{l}\simeq \left(\begin{array}{c}\frac{w_{EE}}{{\tau }_E\left[I_E^{(0)}(t+{\varphi }_n)\right]}\\ \frac{w_{IE}}{{\tau }_I\left[I_I^{(0)}(t+{\varphi }_n)\right]}\end{array}\right){\displaystyle \sum _{m=1}^L}C(n,m)\sqrt{\frac{{\Phi }_E\left[I_E^{(0)}(t+{\varphi }_n)\right]}{N_E}}{\xi }_{E,m}(t).\end{array}& (112)\end{array}$

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) provides coupled equations for the oscillation phases along the chain,

$\begin{array}{ll}\frac{d{\varphi }_n}{dt}=D_{\varphi }^{EI}{\displaystyle \sum _{m=1,\cdots ,L}}C(n,m)[{\varphi }_m(t)-{\varphi }_n(t)]+{\eta }_n(t),& (113)\end{array}$

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

$\begin{array}{ll}\langle {\eta }_n(t){\eta }_m(t^{\prime })\rangle =\delta (t-t^{\prime })\left[\frac{D_I}{N_I}{\delta }_{n,m}+\frac{D_E}{N_E}S(n-m)\right]& (114)\end{array}$

with

$\begin{array}{l}S(n-m)={\displaystyle \sum _{l=1,\cdots ,L}}C(n,l)C(m,l)\\ \begin{array}{l}=A_{\lambda }^2\text{exp}(-\lambda |n-m|)\left[\text{coth}(\lambda )+|n-m|\right],\end{array}& (115)\end{array}$

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.

4.12. 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,

$\begin{array}{ll}{\varphi }_n(t)=\frac{1}{\sqrt{L}}{\displaystyle \sum _{k=0,\cdots ,L-1}}{\tilde{\varphi }}_{q^{(k)}}(t)\text{exp}(i{q}^{(k)}n),& (116)\end{array}$

$\begin{array}{ll}{\eta }_n(t)=\frac{1}{\sqrt{L}}{\displaystyle \sum _{k=0,\cdots ,L-1}}{\tilde{\eta }}_{q^{(k)}}(t)\text{exp}(i{q}^{(k)}n),& (117)\end{array}$

where the L wavevectors are {q^(k) = 2πk/L, k = 0, ⋯, L − 1}. The stochastic terms ${\tilde{\eta }}_{q^{(k)}}$ are Gaussian as linear sums of the η_n(t) and their correlations are obtained from Equations (114, 115) as

$\begin{array}{ll}\langle {\tilde{\eta }}_{q^{(k)}}(t){\tilde{\eta }}_{q^{(l)}}(t^{\prime })\rangle =\delta (t-t^{\prime }){\delta }_{k+l\equiv 0(\text{mod}L)}\left[\frac{D_I}{N_I}+\frac{D_E}{N_E}{\tilde{C}}_{\lambda }^2(q^{(k)}),\right]& (118)\end{array}$

where we have made use of the identity

$\begin{array}{ll}{\displaystyle \sum _n}\text{exp}(iqn)S(n)={\tilde{C}}_{\lambda }^2(q).& (119)\end{array}$

The Fourier transform of the coupling function, ${\tilde{C}}_{\lambda }(q)$, is given in Equation (106). One can note that its normalization implies that ${\tilde{C}}_{\lambda }(q)\to 1$ in the long-wavelength limit q → 0.

Replacement of expressions (116) and (117) into Equation (113) gives

$\begin{array}{ll}\frac{d{\tilde{\varphi }}_{q^{(k)}}}{dt}=D_{\varphi }^{EI}\left[{\tilde{C}}_{\lambda }(q^{(k)})-1\right]{\tilde{\varphi }}_{q^{(k)}}(t)+{\tilde{\eta }}_{q^{(k)}}(t).& (120)\end{array}$

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 ${\tilde{\varphi }}_{q^{(k)}}(t)$,

$\begin{array}{ll}{\tilde{\varphi }}_{q^{(k)}}(t)={\displaystyle {\int }_{-\infty }^t}d{t}^{\prime }\text{exp}\left\{D_{\varphi }^{EI}\left[{\tilde{C}}_{\lambda }(q^{(k)})-1\right](t-t^{\prime })\right\}{\tilde{\eta }}_{q^{(k)}}(t^{\prime }).& (121)\end{array}$

Therefore, the mean square amplitude for the E–I module oscillation phases for a modulation of wavevector q^(k) ≠ 0 is given by

$\begin{array}{ll}\langle |{\tilde{\varphi }}_{q^{(k)}}(t){|}^2\rangle =\frac{D_I/N_I+D_E/N_E{\tilde{C}}_{\lambda }^2(q^{(k)})}{2D_{\varphi }^{EI}\left[1-{\tilde{C}}_{\lambda }(q^{(k)})\right]},& (122)\end{array}$

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 $\langle |{\tilde{\varphi }}_{q^{(k)}}(t){|}^2\rangle$ into an expression for the amplitudes of the firing rate fluctuations $\langle |{\tilde{r}}_{q^{(k)}}(t){|}^2\rangle$ for k ≠ 0, along the lines of the previous calculations. Namely, we find

$\begin{array}{l}{\langle |{\tilde{r}}_{q^{(k)}}(t){|}^2\rangle }_T=\frac{1}{L}{\displaystyle \sum _{n,m}}\text{exp}\left[i{q}^{(k)}(m-n)\right]{\langle r_{n,E}(t)r_{m,E}(t)\rangle }_T\\ \simeq \frac{1}{L}{\displaystyle \sum _{n,m=1}^L}\text{exp}\left[i{q}^{(k)}(m-n)\right]{\displaystyle \sum _{l=-\infty }^{+\infty }}|{\tilde{r}}_l{|}^2\\ \begin{array}{l}\text{exp}\left(-\frac{2{\pi }^2{l}^2}{T}\langle {\left[{\varphi }_n(t)-{\varphi }_m(t)\right]}^2\rangle \right)\end{array}& (123)\end{array}$

with

$\begin{array}{ll}\langle {[{\varphi }_n-{\varphi }_m]}^2\rangle =-2\langle {\varphi }_n{\varphi }_m\rangle =-2{\displaystyle \sum _{j=0}^{L-1}}\text{exp}\left[i{q}^{(j)}(m-n)\right]\langle |{\tilde{\varphi }}_{q^{(j)}}(t){|}^2\rangle .& (124)\end{array}$

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

$\begin{array}{ll}\langle |{\tilde{\varphi }}_q(t){|}^2\rangle ~\left[\frac{D_I}{N_I}+\frac{D_E}{N_E}\right]\frac{\text{cosh}(\lambda )-1}{D_{\varphi }^{EI}{q}^2},q\to 0.& (125)\end{array}$

In this limit, the phase equation (121) reduces to the classical Edwards-Wilkinson equation (Barabási and Stanley, 1995),

$\begin{array}{ll}{\partial }_t\varphi =D_{lw}{\partial }_{xx}\varphi +\eta (x,t),\langle \eta (x,t)\eta (x^{\prime },t^{\prime })\rangle =D_N\delta (x-x^{\prime })\delta (t-t^{\prime }),& (126)\end{array}$

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

$\begin{array}{ll}{D}_{lw}=\frac{D_{\varphi }^{EI}}{2\left[\text{cosh}(\lambda )-1\right]}.& (127)\end{array}$

4.12.2. 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}}_q^{EI}$ 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}}_q^{EI}$ 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).

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

$\begin{array}{ll}\frac{d\delta {\text{I}}_n}{dt}={\displaystyle \sum _{m=1,\cdots ,L}{{L}}_{nm}^{EI}}[{\text{I}}^{(0)}(t)]\delta {\text{I}}_m+{\text{B}}_n& (128)\end{array}$

with the stochastic terms B_n given by,

$\begin{array}{ll}{\text{B}}_n=\left(\begin{array}{c}\frac{1}{{\tau }_E[I_E^{(0)}(t)]}[w_{EE}\sqrt{\frac{{\Phi }_E[I_E^{(0)}(t)]}{N_E}}{\displaystyle {\sum }_{m}C}(n-m){\xi }_{E,m}(t)\\ -w_{EI}\sqrt{\frac{{\Phi }_I[I_I^{(0)}(t)]}{N_I}}{\xi }_{I,n}(t)]\\ \frac{1}{{\tau }_I[I_I^{(0)}(t)]}\left[w_{IE}\sqrt{\frac{{\Phi }_E[I_E^{(0)}(t)]}{N_E}}{\displaystyle {\sum }_{m}C}(n-m){\xi }_{E,m}(t)\right]\end{array}\right)& (129)\end{array}$

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

$\begin{array}{ll}\delta {\text{I}}_n(t)=\frac{1}{\sqrt{L}}{\displaystyle \sum _{k=0,\cdots ,L-1}\delta }{\tilde{\text{I}}}_{q^{(k)}}(t)\exp (i{q}^{(k)}n),& (130)\end{array}$

$\begin{array}{ll}{\xi }_{\nu ,n}(t)=\frac{1}{\sqrt{L}}{\displaystyle \sum _{k=0,\cdots ,L-1}}{\tilde{\xi }}_{\nu ,q^{(k)}}(t)\text{exp}(i{q}^{(k)}n),\nu =E,I.& (131)\end{array}$

The noise Fourier components obey,

$\begin{array}{ll}\langle {\tilde{\xi }}_{\nu ,q^{(k)}}(t){\tilde{\xi }}_{{\nu }^{\prime },q^{(l)}}(t^{\prime })\rangle =\delta (t-t^{\prime }){\delta }_{k+l\equiv 0(\text{mod}L)}{\delta }_{\nu ,{\nu }^{\prime }}& (132)\end{array}$

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

$\begin{array}{ll}\frac{\delta {\tilde{\text{I}}}_{q^{(k)}}}{dt}={{L}}_{q^{(k)}}^{EI}\left[{\text{I}}^{(0)}(t)\right]\delta {\tilde{\text{I}}}_{q^{(k)}}+\delta {\tilde{\text{B}}}_{q^{(k)}}(t).& (133)\end{array}$

The stochastic terms $\delta {\tilde{\text{B}}}_{q^{(k)}}$ are the Fourier components of the B_n (Equation 129),

$\begin{array}{ll}\delta {\tilde{\text{B}}}_{q^{(k)}}=\left(\begin{array}{c}\frac{1}{{\tau }_E[I_E^{(0)}(t)]}[w_{EE}\sqrt{\frac{{\Phi }_E[I_E^{(0)}(t)]}{N_E}}{\tilde{C}}_{\lambda }(q^{(k)}){\tilde{\xi }}_{E,q^{(k)}}(t)\\ -w_{EI}\sqrt{\frac{{\Phi }_I[I_I^{(0)}(t)]}{N_I}}{\tilde{\xi }}_{I,q^{(k)}}(t)]\\ \frac{1}{{\tau }_I[I_I^{(0)}(t)]}\left[w_{IE}\sqrt{\frac{{\Phi }_E[I_E^{(0)}(t)]}{N_E}}{\tilde{C}}_{\lambda }(q^{(k)}){\tilde{\xi }}_{E,q^{(k)}}(t)\right]\end{array}\right).& (134)\end{array}$

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 transform ${\tilde{C}}_{\lambda }(q)$ (Equation 106). The amplitudes of the Fourier mode modulations thus read, similarly to Equations (95, 96),

$\begin{array}{l}\langle |\delta {\tilde{I}}_{E,q^{(k)}}(t){|}^2\rangle ={\displaystyle \sum _{j,l=1,2}{e}_{j,E}^{q^{(k)}}}(t)e_{l,E}^{q^{(k)}}(t)[{\displaystyle {\int }_0^t}du{b}_{jl}(u)\\ +\frac{{\mu }_j{\mu }_l}{1-{\mu }_j{\mu }_l}{\displaystyle {\int }_0^T}du{b}_{jl}(u)],& (135)\end{array}$

where the vectors ${\text{e}}_j^q(t)$ are defined for the matrices ${{L}}_q^{EI}(t)$ as the vectors ${\text{e}}_j^A(t)$ for the matrix ${{L}}_A^{EI}(t)$.

The amplitudes of the rate fluctuations are again obtained by multiplying the expression above by ${\Phi }_E^{\prime }{\left[I_E^{(0)}(t)\right]}^2$, and averaged over the complete limit cycle we find

$\begin{array}{l}{\langle |{\tilde{r}}_{E,q^{(k)}}(t){|}^2\rangle }_T=\frac{1}{T}{\displaystyle {\int }_0^T}dt\{{\displaystyle \sum _{j,l=1,2}{\Phi }_E^{\prime }}{[I_E^{(0)}(t)]}^2{e}_{j,E}^{q^{(k)}}(t)e_{l,E}^{q^{(k)}}(t)\times \\ \left[{\displaystyle {\int }_0^t}du{b}_{jl}(u)+\frac{{\mu }_j{\mu }_l}{1-{\mu }_j{\mu }_l}{\displaystyle {\int }_0^T}du{b}_{jl}(u)\right]\}+\frac{{\langle r_E^{(0)}\rangle }_T}{N_E\Delta T}.& (136)\end{array}$

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.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Joran Deschamps for performing preliminary rate-model simulations and computations on synchronization of E–I modules during an M2 internship. We were also grateful to Frédéric Chavane and Thomas Brochier for instructive comments as well as to Josef Ladenbauer and Nicolas Brunel for drawing our attention, respectively to Augustin et al. (2017) and Rule et al. (2018). We also wish to thank Bard Ermentrout for drawing our attention to Borisyuk et al. (1995) which had escaped our notice.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2020.569644/full#supplementary-material

References

Adesnik, H., and Scanziani, M. (2010). Lateral competition for cortical space by layer-specific horizontal circuits. Nature 464:1155. doi: 10.1038/nature08935

PubMed Abstract | CrossRef Full Text | Google Scholar

Ali, R., Harris, J., and Ermentrout, B. (2016). Pattern formation in oscillatory media without lateral inhibition. Phys. Rev. E 94:012412. doi: 10.1103/PhysRevE.94.012412

PubMed Abstract | CrossRef Full Text | Google Scholar

Amari, S. (1977). Dynamics of pattern formation in lateral-inhibition type neural fields. Biol. Cybernet. 27, 77–87. doi: 10.1007/BF00337259

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashwin, P., Coombes, S., and Nicks, R. (2016). Mathematical frameworks for oscillatory network dynamics in neuroscience. J. Math. Neurosci. 6:2. doi: 10.1186/s13408-015-0033-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Augustin, M., Ladenbauer, J., Baumann, F., and Obermayer, K. (2017). Low-dimensional spike rate models derived from networks of adaptive integrate-and-fire neurons: comparison and implementation. PLoS Comput. Biol. 13:e1005545. doi: 10.1371/journal.pcbi.1005545

PubMed Abstract | CrossRef Full Text | Google Scholar

Badel, L., Lefort, S., Brette, R., Petersen, C. C., Gerstner, W., and Richardson, M. J. (2008). Dynamic iv curves are reliable predictors of naturalistic pyramidal-neuron voltage traces. J. Neurophysiol. 99, 656–666. doi: 10.1152/jn.01107.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Barabási, A.-L., and Stanley, H. E. (1995). Fractal Concepts in Surface Growth. Cambridge: Cambridge University Press.

Google Scholar

Battaglia, D., Brunel, N., and Hansel, D. (2007). Temporal decorrelation of collective oscillations in neural networks with local inhibition and long-range excitation. Phys. Rev. Lett. 99:238106. doi: 10.1103/PhysRevLett.99.238106

PubMed Abstract | CrossRef Full Text | Google Scholar

Battaglia, D., and Hansel, D. (2011). Synchronous chaos and broad band gamma rhythm in a minimal multi-layer model of primary visual cortex. PLoS Comput. Biol. 7:e1002176. doi: 10.1371/journal.pcbi.1002176

PubMed Abstract | CrossRef Full Text | Google Scholar

Bender, C. M., and Orszag, S. A. (2013). Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory. New York, NY: Springer Science & Business Media.

Google Scholar

Börgers, C., and Kopell, N. (2003). Synchronization in networks of excitatory and inhibitory neurons with sparse, random connectivity. Neural Comput. 15, 509–538. doi: 10.1162/089976603321192059

PubMed Abstract | CrossRef Full Text | Google Scholar

Borisyuk, G. N., Borisyuk, R. M., Khibnik, A. I., and Roose, D. (1995). Dynamics and bifurcations of two coupled neural oscillators with different connection types. Bull. Math. Biol. 57, 809–840. doi: 10.1016/S0092-8240(95)80002-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Bressloff, P. C. (2011). Spatiotemporal dynamics of continuum neural fields. J. Phys. A Math. Theor. 45:033001. doi: 10.1088/1751-8113/45/3/033001

CrossRef Full Text | Google Scholar

Brunel, N. (2000). Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J. Comput. Neurosci. 8, 183–208. doi: 10.1023/A:1008925309027

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunel, N., and Hakim, V. (1999). Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. Neural Comput. 11, 1621–1671. doi: 10.1162/089976699300016179

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunel, N., and Hakim, V. (2008). Sparsely synchronized neuronal oscillations. Chaos 18:015113. doi: 10.1063/1.2779858

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunel, N., and Wang, X.-J. (2003). What determines the frequency of fast network oscillations with irregular neural discharges? I. Synaptic dynamics and excitation-inhibition balance. J. Neurophysiol. 90, 415–430. doi: 10.1152/jn.01095.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Buice, M. A., and Chow, C. C. (2013). Beyond mean field theory: statistical field theory for neural networks. J. Stat. Mech. Theory Exp. 2013:P03003. doi: 10.1088/1742-5468/2013/03/P03003

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsaki, G. (2006). Rhythms of the Brain. Oxford: Oxford University Press.

Google Scholar

Capaday, C., Ethier, C., Brizzi, L., Sik, A., van Vreeswijk, C., and Gingras, D. (2009). On the nature of the intrinsic connectivity of the cat motor cortex: evidence for a recurrent neural network topology. J. Neurophysiol. 102, 2131–2141. doi: 10.1152/jn.91319.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Chavane, F., Sharon, D., Jancke, D., Marre, O., Frégnac, Y., and Grinvald, A. (2011). Lateral spread of orientation selectivity in v1 is controlled by intracortical cooperativity. Front. Syst. Neurosci. 5:4. doi: 10.3389/fnsys.2011.00004

PubMed Abstract | CrossRef Full Text | Google Scholar

Chichilnisky, E. (2001). A simple white noise analysis of neuronal light responses. Netw. Comput. Neural Syst. 12, 199–213. doi: 10.1080/713663221

PubMed Abstract | CrossRef Full Text | Google Scholar

Coombes, S. (2005). Waves, bumps, and patterns in neural field theories. Biol. Cybernet. 93, 91–108. doi: 10.1007/s00422-005-0574-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Cross, M. C., and Hohenberg, P. C. (1993). Pattern formation outside of equilibrium. Rev. Mod. Phys. 65:851. doi: 10.1103/RevModPhys.65.851

CrossRef Full Text | Google Scholar

Defelipe, J., Markram, H., and Rockland, K. S. (2012). The neocortical column. Front. Neuroanat. 6:22. doi: 10.3389/978-2-88919-042-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Denker, M., Zehl, L., Kilavik, B. E., Diesmann, M., Brochier, T., Riehle, A., et al. (2018). LFP beta amplitude is linked to mesoscopic spatio-temporal phase patterns. Sci. Rep. 8:5200. doi: 10.1038/s41598-018-22990-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Dumont, G., Ermentrout, G. B., and Gutkin, B. (2017a). Macroscopic phase-resetting curves for spiking neural networks. Phys. Rev. E 96:042311. doi: 10.1103/PhysRevE.96.042311

PubMed Abstract | CrossRef Full Text | Google Scholar

Dumont, G., and Gutkin, B. (2019). Macroscopic phase resetting-curves determine oscillatory coherence and signal transfer in inter-coupled neural circuits. PLoS Comput. Biol. 15:e1007019. doi: 10.1371/journal.pcbi.1007019

PubMed Abstract | CrossRef Full Text | Google Scholar

Dumont, G., Payeur, A., and Longtin, A. (2017b). A stochastic-field description of finite-size spiking neural networks. PLoS Comput. Biol. 13:e1005691. doi: 10.1371/journal.pcbi.1005691

PubMed Abstract | CrossRef Full Text | Google Scholar

Engel, A. K., and Fries, P. (2010). Beta-band oscillations—signalling the status quo? Curr. Opin. Neurobiol. 20, 156–165. doi: 10.1016/j.conb.2010.02.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, B. (1994). Reduction of conductance-based models with slow synapses to neural nets. Neural Comput. 6, 679–695. doi: 10.1162/neco.1994.6.4.679

CrossRef Full Text | Google Scholar

Ermentrout, B. (1998). Neural networks as spatio-temporal pattern-forming systems. Rep. Prog. Phys. 61:353. doi: 10.1088/0034-4885/61/4/002

CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Kleinfeld, D. (2001). Traveling electrical waves in cortex: insights from phase dynamics and speculation on a computational role. Neuron 29, 33–44. doi: 10.1016/S0896-6273(01)00178-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Kopell, N. (1991). Multiple pulse interactions and averaging in systems of coupled neural oscillators. J. Math. Biol. 29, 195–217. doi: 10.1007/BF00160535

CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Terman, D. H. (2010). Mathematical Foundations of Neuroscience, Vol. 35. New York, NY: Springer Science & Business Media.

Google Scholar

Fourcaud-Trocmé, N., Hansel, D., Van Vreeswijk, C., and Brunel, N. (2003). How spike generation mechanisms determine the neuronal response to fluctuating inputs. J. Neurosci. 23, 11628–11640. doi: 10.1523/JNEUROSCI.23-37-11628.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Fries, P. (2009). Neuronal gamma-band synchronization as a fundamental process in cortical computation. Annu. Rev. Neurosci. 32, 209–224. doi: 10.1146/annurev.neuro.051508.135603

PubMed Abstract | CrossRef Full Text | Google Scholar

Gardiner, C. (1985). Stochastic Methods. Springer Series in Synergetics. Berlin: Springer-Verlag.

Google Scholar

Goldobin, D. S., Teramae, J., Nakao, H., and Ermentrout, G. B. (2010). Dynamics of limit-cycle oscillators subject to general noise. Phys. Rev. Lett. 105:154101. doi: 10.1103/PhysRevLett.105.154101

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodman, D. F., and Brette, R. (2009). The brian simulator. Front. Neurosci. 3:192. doi: 10.3389/neuro.01.026.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Hao, Y., Riehle, A., and Brochier, T. G. (2016). Mapping horizontal spread of activity in monkey motor cortex using single pulse microstimulation. Front. Neural Circuits 10:104. doi: 10.3389/fncir.2016.00104

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoppensteadt, F. C., and Izhikevich, E. M. (2012). Weakly Connected Neural Networks, Vol. 126. New York, NY: Springer Science & Business Media.

Google Scholar

Huang, X., Elyada, Y. M., Bosking, W. H., Walker, T., and Fitzpatrick, D. (2014). Optogenetic assessment of horizontal interactions in primary visual cortex. J. Neurosci. 34, 4976–4990. doi: 10.1523/JNEUROSCI.4116-13.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, X., Xu, W., Liang, J., Takagaki, K., Gao, X., and Wu, J. (2010). Spiral wave dynamics in neocortex. Neuron 68, 978–990. doi: 10.1016/j.neuron.2010.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Huntley, G. W., and Jones, E. G. (1991). Relationship of intrinsic connections to forelimb movement representations in monkey motor cortex: a correlative anatomic and physiological study. J. Neurophysiol. 66, 390–413. doi: 10.1152/jn.1991.66.2.390

PubMed Abstract | CrossRef Full Text | Google Scholar

Ince, E. L. (1956). Ordinary Differential Equations. New York, NY: Dover Publications, Inc.

Google Scholar

Kriener, B., Helias, M., Rotter, S., Diesmann, M., and Einevoll, G. T. (2014). How pattern formation in ring networks of excitatory and inhibitory spiking neurons depends on the input current regime. Front. Comput. Neurosci. 7:187. doi: 10.3389/fncom.2013.00187

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuramoto, Y. (1984). Chemical Oscillations, Waves, and Turbulence, Vol. 19. Berlin: Springer-Verlag.

Google Scholar

Levy, R. B., and Reyes, A. D. (2012). Spatial profile of excitatory and inhibitory synaptic connectivity in mouse primary auditory cortex. J. Neurosci. 32, 5609–5619. doi: 10.1523/JNEUROSCI.5158-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Montbrió, E., Pazó, D., and Roxin, A. (2015). Macroscopic description for networks of spiking neurons. Phys. Rev. X 5:021028. doi: 10.1103/PhysRevX.5.021028

CrossRef Full Text | Google Scholar

Mountcastle, V. B. (1997). The columnar organization of the neocortex. Brain 120, 701–722. doi: 10.1093/brain/120.4.701

PubMed Abstract | CrossRef Full Text | Google Scholar

Muller, L., Chavane, F., Reynolds, J., and Sejnowski, T. J. (2018). Cortical travelling waves: mechanisms and computational principles. Nat. Rev. Neurosci. 19, 255–268. doi: 10.1038/nrn.2018.20

PubMed Abstract | CrossRef Full Text | Google Scholar

Nauhaus, I., Busse, L., Carandini, M., and Ringach, D. L. (2009). Stimulus contrast modulates functional connectivity in visual cortex. Nat. Neurosci. 12:70. doi: 10.1038/nn.2232

PubMed Abstract | CrossRef Full Text | Google Scholar

Ostojic, S., and Brunel, N. (2011). From spiking neuron models to linear-nonlinear models. PLoS Comput. Biol. 7:e1001056. doi: 10.1371/journal.pcbi.1001056

PubMed Abstract | CrossRef Full Text | Google Scholar

Palmigiano, A., Geisel, T., Wolf, F., and Battaglia, D. (2017). Flexible information routing by transient synchrony. Nat. Neurosci. 20, 1014–1022. doi: 10.1038/nn.4569

PubMed Abstract | CrossRef Full Text | Google Scholar

Pikovsky, A., Kurths, J., Rosenblum, M., and Kurths, J. (2003). Synchronization: A Universal Concept in Nonlinear Sciences, Vol. 12. Cambridge: Cambridge University Press.

Google Scholar

Prechtl, J., Cohen, L., Pesaran, B., Mitra, P., and Kleinfeld, D. (1997). Visual stimuli induce waves of electrical activity in turtle cortex. Proc. Natl. Acad. Sci. U.S.A. 94, 7621–7626. doi: 10.1073/pnas.94.14.7621

PubMed Abstract | CrossRef Full Text | Google Scholar

Richardson, M. (2007). Firing-rate response of linear and nonlinear integrate-and-fire neurons to modulated current-based and conductance-based synaptic drive. Phys. Rev. E 76:021919. doi: 10.1103/PhysRevE.76.021919

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, J. A., Gollo, L. L., Abeysuriya, R. G., Roberts, G., Mitchell, P. B., Woolrich, M. W., et al. (2019). Metastable brain waves. Nat. Commun. 10:1056. doi: 10.1038/s41467-019-08999-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenbaum, R., and Doiron, B. (2014). Balanced networks of spiking neurons with spatially dependent recurrent connections. Phys. Rev. X 4:021039. doi: 10.1103/PhysRevX.4.021039

CrossRef Full Text | Google Scholar

Rosenbaum, R., Smith, M. A., Kohn, A., Rubin, J. E., and Doiron, B. (2017). The spatial structure of correlated neuronal variability. Nat. Neurosci. 20:107. doi: 10.1038/nn.4433

PubMed Abstract | CrossRef Full Text | Google Scholar

Roxin, A., and Compte, A. (2016). Oscillations in the bistable regime of neuronal networks. Phys. Rev. E 94:012410. doi: 10.1103/PhysRevE.94.012410

PubMed Abstract | CrossRef Full Text | Google Scholar

Rubino, D., Robbins, K. A., and Hatsopoulos, N. G. (2006). Propagating waves mediate information transfer in the motor cortex. Nat. Neurosci. 9:1549. doi: 10.1038/nn1802

PubMed Abstract | CrossRef Full Text | Google Scholar

Rule, M. E., Vargas-Irwin, C., Donoghue, J. P., and Truccolo, W. (2018). Phase reorganization leads to transient β-lfp spatial wave patterns in motor cortex during steady-state movement preparation. J. Neurophysiol. 119, 2212–2228. doi: 10.1152/jn.00525.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanes, J. N., and Donoghue, J. P. (1993). Oscillations in local field potentials of the primate motor cortex during voluntary movement. Proc. Natl. Acad. Sci. U.S.A. 90, 4470–4474. doi: 10.1073/pnas.90.10.4470

PubMed Abstract | CrossRef Full Text | Google Scholar

Schaffer, E. S., Ostojic, S., and Abbott, L. F. (2013). A complex-valued firing-rate model that approximates the dynamics of spiking networks. PLoS Comput. Biol. 9:e1003301. doi: 10.1371/journal.pcbi.1003301

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwalger, T., Deger, M., and Gerstner, W. (2017). Towards a theory of cortical columns: from spiking neurons to interacting neural populations of finite size. PLoS Comput. Biol. 13:e1005507. doi: 10.1371/journal.pcbi.1005507

PubMed Abstract | CrossRef Full Text | Google Scholar

Senk, J., Korvasová, K., Schuecker, J., Hagen, E., Tetzlaff, T., Diesmann, M., et al. (2018). Conditions for traveling waves in spiking neural networks. arXiv 1801.06046.

Google Scholar

Shriki, O., Hansel, D., and Sompolinsky, H. (2003). Rate models for conductance-based cortical neuronal networks. Neural Comput. 15, 1809–1841. doi: 10.1162/08997660360675053

PubMed Abstract | CrossRef Full Text | Google Scholar

Takahashi, K., Kim, S., Coleman, T. P., Brown, K. A., Suminski, A. J., Best, M. D., et al. (2015). Large-scale spatiotemporal spike patterning consistent with wave propagation in motor cortex. Nat. Commun. 6:7169. doi: 10.1038/ncomms8169

PubMed Abstract | CrossRef Full Text | Google Scholar

Takahashi, K., Saleh, M., Penn, R. D., and Hatsopoulos, N. G. (2011). Propagating waves in human motor cortex. Front. Hum. Neurosci. 5:40. doi: 10.3389/fnhum.2011.00040

PubMed Abstract | CrossRef Full Text | Google Scholar

Veit, J., Hakim, R., Jadi, M. P., Sejnowski, T. J., and Adesnik, H. (2017). Cortical gamma band synchronization through somatostatin interneurons. Nat. Neurosci. 20, 951–959. doi: 10.1038/nn.4562

PubMed Abstract | CrossRef Full Text | Google Scholar

Vreeswijk, C. v., and Sompolinsky, H. (1998). Chaotic balanced state in a model of cortical circuits. Neural Comput. 10, 1321–1371. doi: 10.1162/089976698300017214

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X.-J. (2010). Neurophysiological and computational principles of cortical rhythms in cognition. Physiol. Rev. 90, 1195–1268. doi: 10.1152/physrev.00035.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, H. R., and Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1–24. doi: 10.1016/S0006-3495(72)86068-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshimura, K., and Arai, K. (2008). Phase reduction of stochastic limit cycle oscillators. Phys. Rev. Lett. 101:154101. doi: 10.1103/PhysRevLett.101.154101

PubMed Abstract | CrossRef Full Text | Google Scholar