We ground our modeling by comparing it to recordings obtained during a delayed reach-to-grasp task in two macaque monkeys (Denker et al., 2018) that have been made publicly available (Brochier et al., 2018). As in other recordings in similar tasks (Murthy and Fetz, 1992; Sanes and Donoghue, 1993; Rubino et al., 2006), beta oscillations are prominent during the 1 s waiting time, the movement preparatory period, and wane during the movement itself, as shown in Figure 1—figure supplement 1a and c. We thus focus on the recordings during the preparatory period in the following.

In order to develop a model of beta oscillations, assumptions have to be made on the source and mechanism of their generation. We explicitly list below our main ones and their rationale, since they differ from those made in some previous works.

Debate exists regarding the origin of beta oscillations in cortex. On the one hand, the basal ganglia display prominent beta oscillations during different phases of movement (Leventhal et al., 2012). Therefore, one view is that beta oscillations in the motor cortex are simply conveyed to the motor cortex from other structures such as the basal ganglia. On the other hand, sources of beta oscillations have been identified in the cortex (Jensen et al., 2005; Sherman et al., 2016). That different cortical regions have intrinsic oscillation frequencies matching those of their prominent rhythms is also supported by transcranial magnetic stimulation (TMS) perturbation studies. Specifically, the premotor area/supplementary motor area 6 has been found to resonate at ∼30 Hz after stimulation by a short TMS pulse (Rosanova et al., 2009). This seems difficult to explain if the rhythm frequency originates from a distant structure. We thus adopt the intermediate view, previously advocated by Sherman et al., 2016, that beta oscillations are generated by recurrent interactions in the motor cortex but are strongly modulated by inputs from other cortical areas or other structures, notably thalamic ones.

Assuming that beta oscillations originate within M1, the question arises of the main neuronal populations and the recurrent interactions underlying rhythm generation. Two natural candidates are recurrent interactions between interneurons, or recurrent interactions between excitatory and inhibitory populations. Here, one important observation is that during the preparatory period of movement when beta oscillations are prominent, neurons do not fire periodically. As shown in Figure 1—figure supplement 1b, isolated units from the recordings of Brochier et al., 2018 display large CVs. Beta oscillations thus appear to be a collective phenomenon arising from sparse synchronization (Brunel and Hakim, 2008) of different non-oscillating units. In such a regime, oscillations can emerge from recurrent interaction between interneurons, when inhibition is sufficiently strong. However, their frequencies are mainly controlled by the kinetics of synaptic transmission and tend to be in the upper gamma range or higher (Brunel and Wang, 2003; Geisler et al., 2005; de Solages et al., 2008). Lower frequencies arise in networks of E-I units, in which each of the two populations inhibits itself via a slower disynaptic loop. This leads us to assume that beta oscillations are sparsely synchronized oscillations arising from recurrent interaction between excitatory and inhibitory populations. Our aim is to account for recordings obtained from 4 mm × 4 mm multi-electrode arrays with 10 × 10 electrodes. Therefore, the spatial structure of the connectivity needs to be taken into account. We assume that each electrode records the activity of a local neuronal population. For computational tractability, we describe this population activity in a classic way by its firing rate (Wilson and Cowan, 1972; Dayan and Abbott, 2005), but we choose a particular firing rate formulation (see ‘Methods’) that was shown to agree well with direct simulations of spiking network models in previous works (Kulkarni et al., 2020; Ostojic and Brunel, 2011; Augustin et al., 2017). We assume that neurons under different electrodes are mainly connected by excitatory connections and that the connection probability decreases with distance (Huntley and Jones, 1991; Capaday et al., 2009; Hao et al., 2016).

The previous considerations lead us to study the model schematically depicted in Figure 1a, which generalizes one previously analyzed (Kulkarni et al., 2020). It consists of multiple recurrently coupled modules of excitatory (E) and inhibitory (I) neuron populations. These local modules are further coupled by long-range excitatory connections. The different modules are placed on a square 24 × 24 grid with the central 10 × 10 ones corresponding to the recording electrodes (see Figure 1—figure supplement 2). The decrease in synaptic connection probability with distance is described by the function $C(\mathbf{x})$, chosen in the numerical simulations to be Gaussian with range $l\sim 1$ mm, based on anatomical data and microstimulation studies (Huntley and Jones, 1991; Capaday et al., 2009; Hao et al., 2016). Notably, Capaday et al., 2009 report a uniform distribution of boutons along horizontal axons, not a patchy one. A relatively slow propagation speed (∼30 cm/s) along unmyelinated horizontal axons in the monkey cortex has been reported in previous works (González-Burgos et al., 2000; Girard et al., 2001). A comparable speed has been found for horizontal connections in the rat visual cortex (Lohmann and Rörig, 1994). We have taken this into account by introducing a delay proportional to distance, ${\displaystyle D|\mathbf{x}I\mathbf{y}|}$, between the activity ${\displaystyle r_E}$ in excitatory population at location $\mathbf{y}$ and the corresponding postsynaptic depolarization in the E and I population at location $\mathbf{x}$. However, since data are predominantly available for the visual cortex and could be different for the motor cortex, we also compare the results in the following with the case of fast propagation and negligible propagation delays. The model is further detailed in ‘Methods.’ We begin by examining the dynamics of this network when external inputs are constant and the number of neurons in each module is very large. The effects of time-varying inputs and fluctuations due to the finite number of neurons in each module are then addressed in the following sections.

Figure 1 with 4 supplements see all

Download asset Open asset

The network of E-I modules, described by Equations 8–17, has different dynamical regimes as a function of the synaptic weights $w_{EE},w_{IE},w_{EI},w_{II}$ between the excitatory and inhibitory neural populations. This has been extensively studied in the rate model framework (Ermentrout and Terman, 2010; Ashwin et al., 2016) since the pioneering work of Wilson and Cowan, 1972. In the particular case of the present model (‘Methods,’ see also Figure 1—figure supplement 3 for the f-I curve and adaptive time scale used), the different regimes are depicted in Figure 1b, generalizing previous results (Kulkarni et al., 2020) to take into account the finite kinetics of synaptic current and propagation delays. Essential parameters controlling stability are the strength of recurrent excitation, the strength of feedback inhibition through the disynaptic E-I loop, and the strength of autoinhibition of interneurons on themselves, as respectively measured by parameters $\alpha ,\beta$, and $\gamma$ (Equation 34). With other parameters fixed, a steady firing rate state in which excitatory and inhibitory populations fire at $r_E^s$ and $r_I^s$ is unstable when the strength of recurrent disynaptic self-inhibition $\beta$ becomes larger than a threshold value. Mathematical analysis (see ‘Methods’) provides the bifurcation line that separates the steady state regimes from the oscillatory ones. Its exact position depends on propagation delays as depicted in Figure 1b for zero delay and for the propagation delay, here chosen, of 1.3 ms between neighboring modules, see also ‘Methods’ Equation 35. The frequency of oscillations that arise when crossing the bifurcation line at a particular point is also shown in Figure 1b. When recurrent excitation increases, the oscillation frequency decreases. As seen in Figure 1b, propagation delays displace the bifurcation lines but do not change much how the oscillation frequency depends on recurrent excitation. The positions of the bifurcation lines in the $(\alpha ,\beta )$ parameter space also depend on the strength $\gamma$ of recurrent inhibition between interneurons as shown in Figure 1—figure supplement 4a. The oscillation frequency on the bifurcation line also increases as $\gamma$ becomes larger (Figure 1—figure supplement 4). The latency, rise time, and especially the decay time of the synaptic currents are other parameters that influence the oscillatory bifurcation as shown in Figure 1—figure supplement 4d–k. Quite generally, the oscillation frequency decreases when the synaptic time constants increase.

For the parameters of Figure 1 (see Table 1), the dynamics of the network is illustrated for two values of recurrent inhibition, one above (parameters ON, ‘oscillating network’) and one below (parameters SN, ‘steady network’) the bifurcation line point corresponding to beta oscillation frequency (∼20 Hz). As illustrated in Figure 1c and d, for recurrent inhibition stronger than the critical value, the ON model oscillates regularly. For a lower value of recurrent inhibition (SN model), the firing rates of the excitatory and inhibitory populations steadily fire at constant rates. When perturbed away from these rates, the firing rates relax to their stationary values in an oscillatory fashion, transiently exhibiting beta oscillations.

For fixed synaptic coupling, the steady-state firing rates of the excitatory and inhibitory neuronal populations depend on the strengths of the external inputs on those two populations. This is quantitatively illustrated in Figure 1e for synaptic parameters of model SN. Interestingly, when the external input on the inhibitory neurons is increased, their steady firing rate and the steady firing rate of the excitatory population both decrease for a large range of the external input current. (The firing rate of the inhibitory population starts to increase again for very large inputs, when the firing rate of the excitatory population tends to vanish.) This ‘paradoxical suppression of inhibition’ (Tsodyks et al., 1997) is typical of inhibition-stabilized networks and is observed in many cortical areas (Ozeki et al., 2009; Sanzeni et al., 2020) including motor cortex. It arises in the SN model because it is only when recurrent excitation is large, as measured by the parameter α, that its frequency of oscillation is in the beta range (Figure 1b). Namely all networks with ${\displaystyle \alpha >1}$ cannot stably fire at moderate rates without being stabilized by inhibition.

When the external input strengths are changed on the excitatory and inhibitory populations, the steady discharge rates of the two populations change accordingly. They can enter into the parameter regime where a steady discharge is actually unstable and is replaced by an oscillatory one, as shown in Figure 1f. We consider in the present model excitatory inputs to the motor cortex coming from a distal area. The respective strengths of the external inputs, $I_E^{ext}$ and $I_I^{ext}$, on the excitatory and inhibitory populations are proportional to the respective strengths of their synapses, $w_E^{ext}$ and $w_I^{ext}$. Thus, $I_E^{ext}$ and $I_I^{ext}$ vary with the activity of the distal area on a line in the ($I_E^{ext},I_I^{ext}$) plane as depicted in Figure 1f. The respective strengths of $w_E^{ext}$ and $w_I^{ext}$ are taken in the model (Equation 37) such that the network enters the oscillatory regime when the external inputs are decreased as shown in Figure 1f. This conforms to the early observation (Jasper and Penfield, 1949; Penfield, 1954) that external inputs suppress beta oscillations in motor cortex.

When the external input strength on the SN model varies in time, its operating point moves relative to the bifurcation line and correspondingly, beta oscillations are dampened at varying rates. This leads to waxing and waning of their amplitudes as shown in Figure 1g and h. This is also true for the ON model that stands close to the bifurcation line and can move into the non-oscillatory regime when external inputs vary. These networks that operate close to the bifurcation line thus appear as promising candidates to describe the waxing-and-waning of beta oscillations seen in recording data. This leads us to try and compare beta oscillations produced in the model by inputs varying on an appropriate time scale, together with intrinsic fluctuations arising from the finite number of neurons in each module, to those seen in electrophysiological recordings.

Early reports (Murthy and Fetz, 1992) already stressed that beta oscillations were not of constant amplitude but modulated on different time scales. As confirmed in several later studies (see, e.g., Feingold et al., 2015), the average power of beta oscillations changes on a second time scale with different phases of movement, for example, preparation, execution, and post-execution periods. Within individual trials, the amplitude of beta oscillations fluctuates on much shorter time scales with brief bursts of elevated amplitudes of ∼100 ms duration. Examples of LFP spectrograms from two trials of Brochier et al., 2018 are shown in Figure 2a and b and Figure 2—figure supplement 1a and b, respectively, for monkeys L and N. Whereas single trials exhibit short bursts of beta oscillations, trial-averaged spectrograms (Figure 2c and Figure 2—figure supplement 1c) and LFP power spectra (Figure 2d and Figure 2—figure supplement 1d) only show the average increase of oscillation in the beta band and its modulation on a 1 s behavioral time scale (Figure 2c and Figure 2—figure supplement 1c). A more quantitative characterization of the brief high oscillation amplitude events, the beta bursts, here taken to be the 75th percentile of the beta oscillation amplitude distribution (Figure 2e and Figure 2—figure supplement 1e, see also Figure 2—figure supplement 2), is provided by their duration distributions shown in Figure 2f and Figure 2—figure supplement 1f.

Figure 2 with 3 supplements see all

Download asset Open asset

Can the model introduced in the previous section account for these data? At least two sources can be considered for these transient bursts of oscillations fluctuating from trial to trial. The first one is that the number of neurons in each E-I module, the population of neurons in a 200 μm neighborhood of each electrode, is finite, of about 2 × 10⁴ cells given the cell density in cortex. This by itself produces stochastic fluctuations of each module’s activity. These intrinsic fluctuations can be described in a simple way by adding to each neuronal population’s mean firing rate a stochastic component that is inversely proportional to the square root of its population size (see ‘Methods’; Kulkarni et al., 2020; Brunel and Hakim, 1999). The effects of these intrinsic fluctuations for models ON and SN (Figure 1b) can be obtained both mathematically and through computer simulations. The results are displayed in Figure 2—figure supplement 3 for varying numbers $N$ of neurons per module around our estimated value $N=2\times {10}^4$. In this figure and following ones, we plot the excitatory current in the excitatory population as a simple proxy for the LFP in the model. We refer to it as the model ‘proxy-LFP’ in the following. For the parameters SN (Figure 2—figure supplement 3a and b), the model LFP power spectrum amplitude decreases with the number of neurons, but in the whole range $N=2\times {10}^3-2\times {10}^5$, the power spectrum shape is invariant and accurately coincides with the one obtained mathematically by a linear analysis (Equation 51). Compared to the experimental LFP power spectra (Figure 2d and Figure 2—figure supplement 1d), the peak around beta frequencies is less pronounced; furthermore, the power spectrum misses the trough and enhancement for frequencies lower than beta frequency observed in the recordings. For the (oscillatory) ON model (Figure 2—figure supplement 3j and k), the power spectra are more strongly peaked than the experimental ones, with in addition a notable peak at the harmonic frequency of ∼50 Hz coming from the non-sinusoidal shape of the oscillations (Figure 1c and d). This indicates that in either the SN or the ON model, intrinsic fluctuations are not sufficient to account for the experimental spectra of Figure 2d and Figure 2—figure supplement 1d.

Besides intrinsic dynamical fluctuations, it is likely that fluctuations in the motor cortex arise from time-varying inputs coming from other cortical areas or extra-cortical structures, or a mixture of the two. In absence of specific data, we model these as fluctuating inputs with a finite correlation time, more precisely as Ornstein–Uhlenbeck processes with a relaxation time ${\tau }_{ext}$ (‘Methods’) and refer to them simply as ‘external inputs.’ The network dynamics also depend on the way these external inputs target the motor cortex. Are they essentially independently targeting local areas or do they target the motor cortex in a more globally synchronized manner? Again in the absence of specific measurements, we suppose that the external inputs comprise a mixture of global inputs that provide a fraction $c$ of the power, and independent local inputs that provide the remaining fraction $1-c$ (see ‘Methods).

With this simple description, external inputs are characterized by three parameters: their amplitude ${\nu }_{ext}$, their correlation time ${\tau }_{ext}$, and the fraction $c$. We set out to determine these parameters by comparison with the single-electrode LFP trace and with the cross-correlation between the LFP traces of different electrodes.

The single-electrode power spectra and the single-electrode oscillation depend little on the fraction $c$, we thus consider them first (Figure 2—figure supplement 3).

The inclusion of external fluctuations increases the lower frequency part of the power spectra. It also comparatively decreases the power at frequencies higher than beta frequencies as soon as the fluctuation correlation time is in the few tens of millisecond range as shown for models SN (Figure 2—figure supplement 3c–f) and ON (Figure 2—figure supplement 3l–o). The amplitude of external fluctuations should be large enough for their effect to be of larger magnitude than the one produced by intrinsic fluctuations namely ${\displaystyle {\nu }_{ext}>0.2}$ Hz. For the SN model, it should also be small enough (${\displaystyle {\nu }_{ext}<3}$ Hz), not to produce a secondary power enhancement at double the beta frequency as shown in Figure 2—figure supplement 3c. In this range of amplitudes, the power spectrum shape is independent of the external fluctuation amplitude ${\nu }_{ext}$ after normalization and corresponds to the mathematical expression (Ozeki et al., 2009; Figure 2—figure supplement 3d). For the ON model, a sharp secondary peak at twice the beta frequency is present at low external fluctuation amplitude (Figure 2—figure supplement 3l). When the external fluctuations are stronger, this peak is blurred and becomes a wide power enhancement around twice the beta frequency for ${\nu }_{ext}\sim 3$ Hz (Figure 2—figure supplement 3f and o), comparable to the one seen for the SN model at the same external fluctuation amplitude (Figure 2—figure supplement 3c). The external amplitude fluctuations also produce the power spectrum trough and enhancement observed at lower frequencies in experimental recordings (Figure 2d and Figure 2—figure supplement 1d) when their correlation time is not too short (Figure 2—figure supplement 3e–f and n–o).

Examination of the burst durations and amplitudes provides further information and constraints on the parameters. For both SN and ON models, the burst mean duration grows with the input correlation time ${\tau }_{ext}$ (Figure 2—figure supplement 3g and p) while it is weakly dependent on the input fluctuation amplitude ${\nu }_{ext}$ (Figure 2—figure supplement 3i and r). For the SN model, the input fluctuation amplitude ${\nu }_{ext}$ does not strongly affect the burst amplitude and duration distributions (Figure 2—figure supplement 3h–i). Their shapes are moreover close to the experimentally observed ones (Figure 2e and f). On the contrary, for the ON model, the shape of the burst amplitude distribution strongly depends on the amplitude ${\nu }_{ext}$ (Figure 2—figure supplement 3q). It is only for large ${\nu }_{ext}$ that the shape is close to the one obtained for the SN model and resembles the experimental one. The coincidence of the SN and ON distributions, as for the power spectra, is indeed expected when the fluctuations are large enough as compared to the difference between their reference parameters. Since the ON network appears potentially relevant only for large fluctuations, when the distinction before the two network set points is not meaningful, we focus in the following on a network at the set point SN, that is fluctuating around a non-oscillatory set point.

Figure 2i–n and Figure 2—figure supplement 1i–n show the results of model simulations for the SN model with a relaxation time ${\tau }_{ext}=25$ ms and ${\nu }_{ext}=3$ Hz. Both LFP power spectra and beta burst duration as well as amplitude distributions closely resemble the experimental ones.

Having reproduced the single-electrode characteristics of the recording, we turn to the equal-time correlation between signals on different electrodes. As shown in Figure 2g and Figure 2—figure supplement 1g, the correlation between two neighboring electrodes is lower than the autocorrelation of a single-electrode LFP, namely, its variance, but is comparable to the correlation between more distant electrodes. Namely, part of the LFP signal is strongly synchronized between distant electrodes. This is not the case in the model if the external inputs on different modules are uncorrelated, that is, for $c=0$. At the other extreme, when $c=1$, the signals are much more correlated than in the data. As shown in Figure 2o, it is only when the external inputs are almost equally split between local and global contributions that the model correlations are comparable to the experimental ones. For monkey L (N), the value $c=0.4$ ($c=0.3$) is found to provide the best match. (Note that for the comparison with monkey N, we also use slightly changed synaptic couplings SN’ that lead to a slightly lower peak beta oscillation frequency in accordance with experimental data, see Table 1 for retained parameter values.) The autocorrelation of the single-electrode LFP (Figure 2h and Figure 2—figure supplement 1h for monkeys L and N, respectively) is also well accounted for by the model (Figure 2p and Figure 2—figure supplement 1p for parameters SN and SN’, respectively).

Figure 2 and Figure 2—figure supplement 1 provide a comparison between the experimental data for the two monkeys and the model at points SN and SN’, respectively, with the determined characteristics of the external inputs. The model clearly accounts well for the data characteristics, as quantified by the power spectra, burst duration, and amplitude distributions and correlation between different electrodes.

This leads us to investigate how the model compares to the data for other spatio-temporal characteristics of the LFP signals, namely, traveling waves of activity of different types and their characteristics.

Several groups have reported traveling waves on the motor cortex during the preparatory phase of the movement in instructed-delay reaching tasks. In the early study of Rubino et al., 2006, planar waves of beta oscillations in multielectrode LFP recordings were described. Later works (Rule et al., 2018; Denker et al., 2018) have classified the spatio-temporal beta oscillations recorded by multielectrode arrays into different states. Here, we use the data of Brochier et al., 2018 and closely follow the classification scheme proposed by these authors (Denker et al., 2018), distinguishing periods with planar (Pla.) or radial (Rad.) waves, as well as globally synchronized (Syn.) episodes and more random (Ran.) appearing states. The classification criteria for these four states are based on the instantaneous phase and phase gradient spatial distributions as precisely described in ‘Methods.’ In short, negligible phase delays between electrodes characterize the synchronized state. In the other states, the oscillatory activities of some electrodes have significant phase differences and give rise to spatial phase gradients. The phase gradients are tightly aligned for planar waves, pointing inward or outward for radial waves, or are more disordered in the remaining ‘random’ category. We apply the exact same analysis to the experimental recordings and our model simulation data (Figure 3—figure supplement 1 and ‘Methods’).

Figure 3a provides one example of a planar wave in monkey L recordings (Brochier et al., 2018) with the narrow distribution of phase gradient directions shown in Figure 3b. An example of a radial wave is shown in Figure 3c with its outward pointing phase gradients shown in Figure 3d. Analogous examples for monkey N are shown in Figure 3—figure supplement 2.

Figure 3 with 12 supplements see all

Download asset Open asset

For a traveling wave of oscillatory activity at frequency $f$, the local phase velocity is directly related to the inverse of the local phase gradient $\nabla \varphi$. Following previous works (Rubino et al., 2006; Rule et al., 2018; Denker et al., 2018), we define the average wave speed in the presence of phase gradients varying from electrode to electrode as

where the denominator is the gradient magnitude averaged over all $N$ electrode positions. We refer to $v$ simply as the wave speed. The distribution of observed planar wave speeds is shown in Figure 3e. The median planar wave speed is about 30 cm/s but it should be noted that the distribution includes also events with much smaller velocities. The proportion of the different wave types is depicted in Figure 3f.

Can these data be accounted for by our model network? As discussed in the previous section, the single electrode LFP data and their correlations determine the correlation time (${\tau }_{ext}$) of the fluctuating inputs and their repartition ($c$) between global and local inputs, but do not tightly constrain their amplitude. We thus performed a series of model simulations with different input amplitudes ${\nu }_{ext}$ for the SN model as summarized in Figure 3—figure supplement 3a. The four wave states can be observed for most ${\nu }_{ext}$. Examples of planar and radial waves are provided in Figure 3g–j. However, the proportion of the different wave types and the planar wave speeds greatly vary with ${\nu }_{ext}$ (Figure 3—figure supplement 3). In the SN model, the input fluctuations promote oscillations. The stronger the amplitude ${\nu }_{ext}$, the more developed the oscillations in neural activity and the better they synchronize. Thus, increasing the noise amplitude first increases the proportion of synchronized activity and the propagation speed of planar waves. At higher noise amplitudes, the desynchronizing effect of noise counteracts its oscillation-promoting effect and the wave speed decreases with increasing input amplitude. This leads in the SN model for ${\nu }_{ext}$ in the range $1-3$ Hz to a proportion of a planar wave states of a few percent with a planar wave speed of about 30 cm/s, as observed in the experimental recordings (Figure 3e–f). As a consequence, the proportion of planar and radial waves as well as synchronized activity increases with ${\nu }_{ext}$. The effect of increasing the amplitude of the input fluctuations is quite different for the ON model (Figure 3—figure supplement 3b). In this case, the oscillatory activity is self-generated and the oscillations between different modules are well-synchronized without fluctuating inputs. The local inputs tend to desynchronize the different modules and create oscillatory phase differences between them. The proportion of fully synchronized states thus decreases with increasing ${\nu }_{ext}$. The speed of planar waves that is inversely proportional to the magnitude of phase gradients (Equation 1) also decreases. For an input fluctuation amplitude ${\nu }_{ext}=3$ Hz, the SN model agrees well with the distribution of observed wave speeds and the repartition between the different wave types for monkey L as shown in Figure 3, although radial waves seem somewhat underrepresented in the model network. This is also the case for monkey N as shown in Figure 3—figure supplement 2. The detailed statistics of the duration, wave speeds, oscillation amplitude, and frequency for the four wave states are detailed in Figure 3—figure supplement 4 and Figure 3—figure supplement 5 for monkey L recordings and the SN model and in Figure 3—figure supplement 6 and Figure 3—figure supplement 7 for monkey N recordings and the SN’ model. The SN model simulated in Figure 3 and Figure 3—figure supplement 5 is taken as our reference model in the following and referred to as such. All its parameters are provided in Table 1.

The model and the simulation results suggest that the observed spatio-temporal patterns arise from the opposing effects of synchronization by long-range excitatory connectivity and global inputs and desynchronization by local inputs and the intrinsic stochasticity arising from the finite number of neurons in each module.

Before proceeding, it appears worth discussing and scrutinizing different aspects of our simulations. At a technical level, the simulation results come from the central 10 × 10 modules of a larger 24 × 24 grid with fixed boundary conditions (Figure 1—figure supplement 2 and ‘Methods’). In order to test the effect of the boundary conditions on the results, we performed a set of simulations and measurements with the same setting but with periodic boundary conditions. As shown in Figure 3—figure supplement 8, this did not significantly change the auto- and cross-correlations of module activities, nor the distributions of different wave types or of planar wave speeds. Therefore, one can conclude that the results are not dependent on the implemented boundary condition.

At a more conceptual level, one can wonder whether all the model components are necessary. A first natural question is whether intrinsic stochasticity would be enough to create spatio-temporal patterns similar to experiments without the necessity for desynchronization coming from local inputs. In order to address this possibility, we performed a set of simulations without local inputs. We increased the intrinsic noise as compared to our reference model ($N=20000$) by lowering the numbers of neurons per module. As shown in Figure 3—figure supplement 9, for a number of neurons of $N=2000$ (or higher) per module, the dynamics is too synchronized to account for the experimental data, the cross-correlation of the activity of different modules is very high even for distant modules and the speed of propagation of planar waves is also quite high. If the amplitude of the fluctuating global inputs to the network is decreased to ease the desynchronization of different modules, the network stays too close to the operating point SN and planar waves are no longer seen (not shown). Results comparable to our reference model can only be obtained when the number of neurons per module is lowered to $N=200$ (Figure 3—figure supplement 9a–d). This of course would be quite a low number for the $0.1-0.2$ mm² area of motor cortex that each module is meant to represent, even if restricted to the cortical layers (e.g., II/III) where the dynamics could mainly take place. We found in a previous work (Kulkarni et al., 2020) that our reduced description of local noise (Equation 10 and ‘Methods’) agreed well with simulations of spiking networks with dense local connectivity and sparse long-range connectivity. Whether our description of intrinsic stochasticity underestimates the fluctuations in a more sparsely connected network would have to be checked by large direct numerical simulations of spiking networks. Even if this was the case, the very low number of neurons in a module that we found to be necessary to reproduce the experimental recordings would still support the role of local inputs in creating dephasing between different local regions of motor cortex and planar waves.

A second legitimate question concerns the role of spatial correlations, produced by the connectivity spatial structure, in the observed spatio-temporal dynamics. Could the results be simply the product of random fluctuations in the activities of different modules? In order to address this question, we simulated the model of Figure 3 and shuffled the proxy-LFP signals of the different modules after the simulation. We then reclassified the different types of waves in the shuffled signals. As a comparison, we followed the same procedure for the recordings of monkey L, namely, we classified the occurence of different types of waves after shuffling the LFPs recorded on different electrodes. The results are shown in Figure 3—figure supplement 10 and very similar for the simulations and the experimental recordings. Shuffling has little effect on the proportion of synchronized and random waves. This is intuitively understandable since synchronized signals remain synchronized after shuffling, and shuffling cannot synchronize unsynchronized signals. More interestingly, shuffling entirely suppresses the appearance of planar waves. Therefore, the spatial correlation of different electrode signals, a reflection of the underlying long-range connectivity, is crucial for the appearance of planar waves. This conclusion is further supported by varying the connectivity range as described below.

The presence of planar waves is one remarkable feature of the multi-electrode recordings. It is thus worth analyzing some of their properties in the model.

The mean speeds of planar waves in the recordings and in the model are of a few tens of cm/s (Figure 3e and k and Figure 3—figure supplement 2e and k), which is comparable to the propagation speed along non-myelinated horizontal fibers (González-Burgos et al., 2000; Girard et al., 2001). One might thus think that both are tightly linked. This is actually not the case. An equally good agreement with the recording data can be found in the model in the absence of propagation delay, that is, for $D=0$. The corresponding wave pattern statistics are shown in Figure 3—figure supplement 11, and agree as well as the model with delay with the data for monkey L (Figures 2 and 3).

One can wonder how the range $l$ of long-range excitatory connections influences planar wave appearance and properties. This is more easily addressed in the absence of propagation delay ($D=0$) since, in this case, with other parameters fixed, the operating point of the network (Figure 1b) does not change when $l$ is varied, which is not the case when propagation delays are significant. The results of simulations for different $l$ in the model of Figure 3—figure supplement 11 are shown in Figure 3—figure supplement 12.

First, when long-range excitatory connections are absent ($l=0$), both synchronized states and planar waves are absent and replaced by random states (Figure 3—figure supplement 12a). This is linked to the significant decrease in the correlation of the activities of different modules, shown by the cross-correlation plots in Figure 3—figure supplement 12a. This confirms that long-range excitatory connections are required for planar wave existence.

Second, when $l$ is halved (Figure 3—figure supplement 12b) as compared to the reference model of Figure 3 (i.e., $l=1$ in module distance units or 400μm in real units), both synchronized states and planar waves exist but in reduced proportions as compared to the reference case. Correlatively, the proportion of random states increases. The speed of planar waves is also strongly reduced.

Finally, when $l$ is increased by 50% (Figure 3—figure supplement 12c) as compared to the reference model, synchronization between distant modules is increased, the synchronized regime prevails and planar wave occurrence is reduced.

Synchronized and random states are much more prevalent than planar waves and the number of these two types of patterns is comparable ($N=3906$ synchronized events vs.${\displaystyle N=4293}$ random events in 72 simulations of 10 s duration each). It is interesting that planar waves predominantly arise from synchronized states and moreover that synchronized states tend to persist in the background of planar waves. Namely, of 421 planar waves in SN model simulations that started after a synchronized state, 364 (resp. 57) were followed by a synchronized state (resp. random state), clearly a much greater proportion than if the ending state was drawn at random ($p=3\times {10}^{-57}$). In contrast, of the 161 planar waves that started from a random state, 92 (69) ended in a random (synchronized) state, a proportion compatible with random draw ($p=0.22$).

It is difficult to decipher the structure of the inputs associated to planar waves from the existing experimental recordings. The model allows one to more easily address this question for the simulations. To this end, different planar waves episodes in the reference model of Figure 3 were averaged after aligning their oscillatory phase maps, as described in ‘Methods.’ In brief, the different events were shifted in time, so that to make their phases at at half-duration coincide, and, where space-dependent maps of planar waves were considered, rotated in space to make the wave direction coincide. This was similarly performed for the experimental recordings. The structure of the mean planar wave obtained by averaging planar wave events propagating along the principal four axis of the grid is shown in Figure 4. The average phase (Figure 4a) and proxy-LFP (Figure 4b) clearly show the planar wave propagation. Similar mean spatio-temporal phase maps are obtained by distinguishing three groups of planar wave velocities before averaging (Figure 4—figure supplement 1).

Figure 4 with 2 supplements see all

Download asset Open asset

We furthermore aligned planar wave events in time based on wave onset. This temporal alignment procedure can be used to compute the average parameter ${\sigma }_p$, quantifying synchronization between different module activities (‘circular variance of phase,’ see Equation 6 in ‘Methods’), associated to the average planar wave around wave onset. As shown in Figure 4c, ${\sigma }_p$ rises just before the planar wave to reach a high synchronization value. Similar phase and LFP profiles are obtained with experimental data (Figure 4—figure supplement 2a and b) together with the rise of the synchronization parameter ${\sigma }_p$ (Figure 4—figure supplement 2c).

The average inputs on the modules, which are unavailable in the recordings, can be computed in the model. The global input to the whole network is found to be shifted to negative values before the appearance of a planar wave (Figure 4d), effectively reducing the drive to both the excitatory and the inhibitory populations of all modules and pushing the network towards its self-sustained oscillating regime (see Figure 1f). When the local inputs are averaged in the same manner, no significant feature emerges. This is presumably due to the large variability of local input patterns that can give rise to a planar wave.

In order to go beyond correlation and test whether a shift of global inputs can help create planar waves, we performed simulations of the reference model in which, in addition to the fluctuating local and global inputs (Figure 4e), steps of current were injected in all modules of the network, as depicted in Figure 4f and g. The current steps were chosen of a duration of 50 ms and of an amplitude comparable to the global input shift associated to planar waves. The injection of these current steps resulted in a 50% increase of the number of planar wave episodes (Figure 4f). This was not simply a result of the perturbation of the network since injection of steps with the opposite sign of the injected current on the contrary strongly reduced the appearance of planar waves (Figure 4g).

Our model network provides a satisfactory description of LFP correlations and wave statistics in the data of Brochier et al., 2018 when averaged over directions. The long-range synaptic connection probability that we implemented decreases with distance and is independent of orientation. With this choice, we checked that the direction of observed planar waves does not significantly depend on orientation as shown in Figure 5a–c (namely, the anisotropy created by our choice of a rectangular grid is weak). However, it was reported by Rubino et al., 2006 that waves tend to propagate along preferred axes on the cortical surface with respect to sulcal landmarks and that the orientation of the preferred axes of propagation also vary between different regions, for example, between primary motor cortex and dorsal premotor cortex. For a given preferred axis, propagation in both directions was observed, for example, both rightward and leftward propagating waves were recorded. In order to investigate the possible origin of this observation, we assessed the effect of introducing some anisotropy in our model. We considered two simple possibilities, namely, that intrinsic connectivity was anisotropic or that local inputs target the motor cortex in an anisotropic manner (Figure 5). We consider them in turn.

Figure 5 with 1 supplement see all

Download asset Open asset

First, in order to analyze the effect of anisotropic connectivity in our model, we simulated a modified version of the model in which the connectivity was of twice longer range in the y-direction than in the x-direction, as illustrated in Figure 5d. This did not change much the distribution of the observed wave types (Figure 5e). In contrast, Figure 5f shows the measured distribution of the directions along which planar waves propagate when connections are anisotropic. The fractions of waves propagating along the x- and y-axis are clearly different. There is a large predominance of waves observed to propagate along the x-axis, namely, along the axis where the connectivity is shorter range. Longer-range connectivity promotes longer-range synchronization along the y-axis and synchronized states. The weaker connectivity along the x-axis allows for the formation of the larger phase gradients along the x-axis and planar waves.

Second, we wished to analyze the effect of a possible anisotropy of the local input spread. In our reference model, local inputs target only one module. In order to be able to introduce anisotropy, we first generalized the local input description. We considered the case where local inputs predominantly target a given module at position $\mathbf{x}$ but also the neighboring ones, at positions $\mathbf{y}$, with smaller weights as described by a kernel function $G(\mathbf{x},\mathbf{y})$ (Equation 22 and ‘Methods’). As for connectivity, $G$ was taken to be Gaussian. Simulation results for an isotropic $G$ of range $l^{(I)}$ are shown in Figure 5g–i and Figure 5—figure supplement 1. When $l^{(I)}$ is large enough so that local inputs significantly target neighboring modules, it has a strong synchronizing effect as seen both in the large proportion of synchronized states and the high cross-correlation of the activity of different modules (Figure 5—figure supplement 1a and b for $l^{(I)}=2$ in grid units or 800 μm in real units and Figure 5—figure supplement 1c and d for $l^{(I)}=1$ or 400 μm in real units). For local inputs of narrower footprints (Figure 5g–i and Figure 5—figure supplement 1e and l for $l^{(I)}=0.5$ or 200 μm in real units), the results are analogous to the ones obtained with our reference model, as intuitively expected.

With this generalized description of local inputs, we can now consider an anisotropic kernel $G(\mathbf{x},\mathbf{y})$ (Equation 22) so that local input targeting is more spread out in one direction than in the orthogonal one. This is represented in Figure 5j for the case that we simulated where local inputs are more spread out in the y-direction ($l_y^{(I)}=1$ or 400 μm in real units) than in the x-direction ($l_x^{(I)}=0.5$ or 200 μm in real units). This anisotropic targeting of local inputs produces a clear anisotropy in the distribution of planar wave propagation directions, as shown in Figure 5l. Planar waves preferentially propagate along the x-axis where the inputs are confined to a narrower region and gradients in the phase of the oscillatory activity are more easily produced.

Thus, both anisotropy in connectivity and in the targets of local inputs can produce an anisotropy in wave speed direction. However, Gatter et al., 1978 reported a more spatially extended connectivity in layers 2/3 and 5 of motor cortex in non-human primates in the rostro-caudal direction. Rubino et al., 2006 found that planar waves predominantly propagated along the same rostro-caudal axis. Given our results, this argues against connectivity anisotropy being the dominant factor determining planar wave speed direction. It remains to be investigated whether local input targets are more spread out along the medio-lateral dimension, as the model suggests.

We analyzed the data that was made publicly available and described in detail in Brochier et al., 2018. We content ourselves here to describe their main features for the convenience of the reader. The data were obtained from two macaque monkeys (L and N) trained to perform one of four movements at a GO signal. In brief, in each trial, the animal was presented for 300 ms with a visual cue that provided partial information on the movement to be performed. It had to wait for 1 s, before the GO signal that also brought the missing information about the rewarded movement. The different phases of a trial are depicted in Figure 1—figure supplement 1. The neural activity was recorded during the task with a 10 × 10 square multielectrode Utah array, with neighboring electrodes separated by 400 μm, implanted in the primary motor cortex (M1) or premotor cortex contralateral to the active arm. We make use of two published recording sessions (one for monkey L of 11:49 min/135 correct trials, one for monkey N of 16:43 min/141 correct trials) as fully detailed in Brochier et al., 2018.

An overview of our data analysis protocol is provided in Figure 3—figure supplement 1. The different steps are detailed below.

The motor cortex is described as a collection of recurrently connected excitatory (E) and inhibitory (I) neuronal populations of linear size $\sim 400\mu m$, comparable to the one of a cortical column. These E-I modules are connected by long-range excitatory connections with a connection probability decaying with the distance between modules. The neural activity of each E-I module is represented at the level of its excitatory and inhibitory neuronal populations in the rate-model framework (Wilson and Cowan, 1972; Dayan and Abbott, 2005). This eases numerical simulations that would otherwise be extremely demanding for a comparable spatially structured network of a few hundred modules with a few tens of thousand spiking neurons each. Additionally, the rate model description also eases analytical computations. We choose a rate model description with an adaptive time scale (Ostojic and Brunel, 2011; Augustin et al., 2017). It was shown to quantitatively describe networks of stochastically spiking exponential integrate-and-fire (EIF) neurons, either uncoupled and receiving identical noisy inputs (Ostojic and Brunel, 2011), or coupled by recurrent excitation (Augustin et al., 2017), as well as sparsely synchronized oscillations of recurrently coupled spiking E-I modules of EIF neurons (Kulkarni et al., 2020). The adaptive time scale rate model describes the activity of a population of EIF neurons as

with ${\mathrm{\Phi }}_{\sigma }[I]$ the f-I curve of an EIF neuron with white noise current input of mean $I$ and noise strength $\sigma$ (see Kulkarni et al., 2020 for a detailed description). The time scale $\tau (I)$ is referred to as ’adaptive’ because it varies with the current $I$. The function $\tau (I)$ is chosen here, as in Kulkarni et al., 2020 and precisely described there, to best fit the response to oscillatory inputs in a wide frequency range of 1 Hz to 1 kHz, of an EIF neuron subjected also to a white noise current of mean $I$ and strength $\sigma$. The tabulated functions ${\mathrm{\Phi }}_{\sigma }[I]$ and $\tau (I)$ are plotted in Figure 1—figure supplement 3. They are given in Kulkarni et al., 2020 and are also provided here as Source data for the convenience of the reader.

We generalize the rate model description Equation 7 to a two-dimensional network of E-I modules coupled by long-range excitation,

where the index denotes the module position on a rectangular grid. The firing rates are related to the currents by

The ${\xi }_A(\mathbf{x},t)$ are independent unit amplitude white noises for each neuronal population in each module $⟨{\xi }_A(\mathbf{x},t){\xi }_B({\mathbf{x}}^{\prime },t^{\prime })⟩={\delta }_{A,B}{\delta }_{\mathbf{x},{\mathbf{x}}^{\prime }}\delta (t-t^{\prime })$ and Ito’s prescription is used to define Equation 10 (Kulkarni et al., 2020). These stochastic terms account, in the firing rate framework, for the fluctuations due to the finite numbers $N_E$ of excitatory neurons and $N_I$ of inhibitory in each E-I module (Brunel and Hakim, 1999). For the numerical computations, both ${\mathrm{\Phi }}_E$ and ${\mathrm{\Phi }}_I$ are taken equal to the function ${\mathrm{\Phi }}_{\sigma }[I]$ provided in the Source data and plotted in Figure 1—figure supplement 3a.

The currents $I_E^{ext}(\mathbf{x},t)$ (resp. $I_I^{ext}(\mathbf{x},t)$) represent inputs external to the motor cortex, possibly position and time-dependent, targeting the excitatory (resp. inhibitory) population of the E-I module at position $\mathbf{x}$,

With our normalization choice ${\displaystyle {\nu }_{ext}}$, has the dimension of a frequency and can be interpreted as the discharge rate amplitude of the external inputs. We model the fluctuations of external inputs as stochastic O-U processes with global and independent components,

The currents $I_{EE}^{syn}(\mathbf{x},t)$, $I_{EI}^{syn}(\mathbf{x},t)$ (resp. $I_{IE}^{syn}(\mathbf{x},t)$, $I_{II}^{syn}(\mathbf{x},t)$) represent the recurrent excitatory and inhibitory inputs on the excitatory (resp. inhibitory) population of the module at position $\mathbf{x}$. The recurrent inputs depend on the firing rates of the neuronal populations of the different motor cortex modules and on the kinetics of the different synapses. Namely,

The kinetic kernels $S_E(t),S_I(t)$ include synaptic current rise times, ${\tau }_r^E,{\tau }_r^I$, decay times, ${\tau }_d^E,{\tau }_d^I$, and latencies, ${\tau }_l^E,{\tau }_l^I$,

where $\theta (t)$ denotes the Heaviside function, ${\displaystyle \theta (t)=1,t>0}$ and 0 otherwise. The kernels $S_A,A\in \{E,I\}$ are normalized such that:

We have supposed, for simplicity, that the kinetics of the synaptic currents depend on their excitatory or inhibitory character but not on the nature of their post-synaptic targets (e.g., $I_{EI}^{syn}$ and $I_{II}^{syn}$ have the same kinetics). Instead of using the kinetic kernels $S_A(t),A\in \{E,I\}$, one can equivalently compute the synaptic currents by introducing supplementary variables $J_{AE},J_{AI}$,

In Equation 13 or Equation 16, the probability that an excitatory neuron at position $\mathbf{y}$ targets a neuron at position $\mathbf{x}$ is represented by the function $C(\mathbf{x}-\mathbf{y})$ normalized such that

The delay $D|\mathbf{x}-\mathbf{y}|$ accounts for the signal propagation time along horizontal non-myelinated fibers in the motor cortex (González-Burgos et al., 2000; Girard et al., 2001). The mathematical expressions below are written for a general function $C(\mathbf{x})$. For all the numerical computations, except those shown in Figure 5d–f, the following isotropic Gaussian function is taken,

In Figure 5d–f, we consider instead the anisotropic function,

with $\mathbf{x}=(x,y)$.

In Figure 5, an anisotropic distribution of the local input targets is considered. This is modeled by replacing Equation 12 by

The kernel $G$ is taken to be Gaussian,

where the normalization $Z$ is chosen such as to have the same local variance of the noise as in Equation 12. In Figure 5g–i and Figure 5—figure supplement 1, an isotropic kernel $G$ is taken with $l_x^I=l_y^I=l^I$.

Our analysis starts by considering the deterministic version of the model described by Equations 8–17, that is, the limit $N_E\to \mathrm{\infty }$, $N_I\to \mathrm{\infty }$, in which the amplitudes of the stochastic terms vanish in Equation 10. All module populations are supposed to fire steadily in time at constant rates $r_E^s$ and $r_I^s$. This requires external inputs that are constant in time, that is, with ${\sigma }_E^{ext}={\sigma }_I^{ext}=0$ (Equation 11), and of specific magnitudes that we first determine.

We then assess the stability of this steady state. This provides the bifurcation diagrams of Figure 1. When the steady state is stable, we compute the effect of fluctuations arising both from the finite number of neurons in each E-I module and from the time variation of the external inputs. This provides the analytic curves for the power spectra and correlations of currents, shown in Figure 2, Figure 2—figure supplement 1, and Figure 2—figure supplement 3.

We consider the network described by Equations 8 and 9 with the kinetics of the synaptic currents given by Equation 14. We include the stochastic effects arising from the finite number of neurons in each module by using the stochastic description (Equation 10) of the instantaneous firing rates of the excitatory and inhibitory module neuronal populations. We also include external input fluctuations as described by Equations 11; 12.

We treat these two kinds of stochastic effects as perturbations of the steady dynamics and fully characterize the stochastic dynamics of the network at the linear level.

Linearizing the currents around their stationary values gives

Translation invariance leads us to search for $\delta I_E$ and $\delta I_I$ in Fourier space:

with an analogous expansion for $\delta I_I$. The linearized Equations 8; 9 then read, with a vectorial notation

where the 2 × 2 matrix ${\tilde{L}}_{EI}(\mathbf{q},i\omega )$ is given in Equation 30. The two components of the stochastic forcing term ${\displaystyle \mathbf{F}(\mathbf{q},i\omega )}$ read

where ${\sigma }_A^{ext}$ is given by Equations 11; 37. Solution of the linear system (40) provides the expression of the Fourier components of the currents. For the excitatory current, one obtains

with ${\displaystyle W(\mathbf{q},i\omega )}$ defined in Equation 32 and ${\displaystyle V_E(\mathbf{q},i\omega )}$ given by,

The Fourier components of the input fluctuations read

with the white noise averages,

The short-hand notation ${\displaystyle {\delta }^2(\mathbf{q})}$ has been used for the two-dimensional $\delta$-function $\delta (q_x)\delta (q_y)$. Similarly, one has for the finite size noise averages:

The current–current correlation function is obtained by averaging the product of the currents (Equation 42) over the noise with the help of Equations 45; 46. One obtains,

with

and

This provides the expression of the current–current correlation in real space:

One can check that the cross-correlation is a real function, as it should, since $S_{EE}^{ext}$ and $S_{EE}^N$ are related to their complex conjugates by ${\displaystyle S_{EE}^{ext\ast }(\mathbf{q},\omega )=S_{EE}^{ext}(-\mathbf{q},-\omega )}$ and ${\displaystyle S_{EE}^{N\ast }(\mathbf{q},\omega )=S_{EE}^N(-\mathbf{q},-\omega )}$, a symmetry inherited from the function ${\displaystyle C(\mathbf{q},\omega )}$ (Equation 31).