Beta oscillations and waves in motor cortex can be accounted for by the interplay of spatially structured connectivity and fluctuating inputs
Abstract
The beta rhythm (13–30 Hz) is a prominent brain rhythm. Recordings in primates during instructeddelay reaching tasks have shown that different types of traveling waves of oscillatory activity are associated with episodes of beta oscillations in motor cortex during movement preparation. We propose here a simple model of motor cortex based on local excitatoryinhibitory neuronal populations coupled by longrange excitation, where additionally inputs to the motor cortex from other neural structures are represented by stochastic inputs on the different model populations. We show that the model accurately reproduces the statistics of recording data when these external inputs are correlated on a short time scale (25 ms) and have two different components, one that targets the motor cortex locally and another one that targets it in a global and synchronized way. The model reproduces the distribution of beta burst durations, the proportion of the different observed wave types, and wave speeds, which we show not to be linked to axonal propagation speed. When the longrange connectivity or the local input targets are anisotropic, traveling waves are found to preferentially propagate along the axis where connectivity decays the fastest. Different from previously proposed mechanistic explanations, the model suggests that traveling waves in motor cortex are the reflection of the dephasing by external inputs, putatively of thalamic origin, of an oscillatory activity that would otherwise be spatially synchronized by recurrent connectivity.
Editor's evaluation
This article makes a valuable contribution to the field. Here the authors have developed a convincing model to characterize the generation of motor cortex β oscillations. Using the model, the authors are able to recapitulate several properties observed experimentally. Given the longstanding interest in motor cortical β oscillations and current interest in traveling waves, this article will be of significant interest to the neuroscience community.
https://doi.org/10.7554/eLife.81446.sa0Introduction
Neural rhythms are one of the most obvious feature of neural dynamics (Buzsáki, 2006). They have been recorded for more than 90 years in human (Berger, 1929) and they are a daily tool for the diagnosis of neurological dysfunction. Classic studies have shown that neural rhythms depend on neural structures and the behavioral state of the animal (Adrian, 1935). Examples include the theta rhythm, the gamma rhythm in the visual cortex, or the fast 160–200 Hz cerebellar rhythm. However, in spite of their ubiquity and common diagnostic use, neural rhythms remain a somewhat mysterious feature of the brain dynamics. It remains to better understand how they are created and what they are a reflection of (Wang, 2010).
The beta rhythm consists of oscillations in the 13–30 Hz range. The motor cortex was found to be one of its main locations in early recordings in human subjects (Jasper and Penfield, 1949; Penfield, 1954). It was recorded in cats during motionless focused attention, when fixating a mouse behind a transparent screen (Bouyer et al., 1987). Subsequent studies in monkeys trained to perform an instructeddelaytask (Murthy and Fetz, 1992; Sanes and Donoghue, 1993) showed that beta oscillations develop in the motor and premotor cortex during the movement preparatory period and wane during the movement itself, in agreement with earlier observations in human (Jasper and Penfield, 1949; Penfield, 1954). It was also noted that beta oscillations did not have a constant amplitude in time but rather appeared as synchronized bursts of a few cycles of oscillations, often synchronized on electrodes 1–2 mm apart. More recently, it was observed using multielectrode arrays that beta oscillations can come as planar (Rubino et al., 2006) or more complex (Rule et al., 2018; Denker et al., 2018) waves propagating horizontally on the motor cortex. Our aim in this paper is to develop a mechanistic framework for these observations and to compare it to available recordings in monkeys performing an instructed delayed reachtograsp task (Denker et al., 2018; Brochier et al., 2018). The model is based on recurrent interaction between local populations of excitatory (E) and inhibitory (I) populations of neurons coupled by longer range excitation. The local EI module is wellknown to exhibit oscillatory dynamics when the recurrent interaction within the E population is sufficiently strong. Longrange excitation between these local modules has been shown in previous works to synchronize the oscillatory dynamics of local EI modules (Ermentrout and Kopell, 1991; Kulkarni et al., 2020). Comparison with recordings leads us to propose that the whole network is close to the oscillation threshold and that fluctuating inputs from other cortical areas or other structures, such as the thalamus, power the burst of beta oscillations. We show that such a spatially extended network submitted to both local and global external inputs exhibits waves of different types that closely resemble those recorded in monkey motor cortex. Our analysis makes testable, specific predictions on the structure of external inputs. More generally, it highlights the dynamical interplay between fluctuating inputs and intrinsic dynamics shaped by spatially structured connectivity, which is worthy of further experimental investigation.
Data and modelling assumptions
We ground our modeling by comparing it to recordings obtained during a delayed reachtograsp 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 nonoscillating 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 EI 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 multielectrode 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).
Results
EI modules connected by longrange excitatory connections exhibit oscillatory activity at beta frequency
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 longrange 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álezBurgos 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, $D\mathbf{x}I\mathbf{y}$, between the activity $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 timevarying inputs and fluctuations due to the finite number of neurons in each module are then addressed in the following sections.
The network of EI 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 fI 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 EI 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 selfinhibition $\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 steadystate 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 inhibitionstabilized 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 $\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 nonoscillatory regime when external inputs vary. These networks that operate close to the bifurcation line thus appear as promising candidates to describe the waxingandwaning 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.
Fluctuation of external inputs and finitesize stochasticity produces model LFP signals statistically similar to the recorded ones
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 postexecution 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, trialaveraged 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.
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 EI module, the population of neurons in a 200 μm neighborhood of each electrode, is finite, of about 2 × 10^{4} 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 ‘proxyLFP’ 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 nonsinusoidal 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 timevarying inputs coming from other cortical areas or extracortical 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 $1c$ (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 singleelectrode LFP trace and with the crosscorrelation between the LFP traces of different electrodes.
The singleelectrode power spectra and the singleelectrode 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 ${\nu}_{ext}>0.2$ Hz. For the SN model, it should also be small enough (${\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 nonoscillatory 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 singleelectrode characteristics of the recording, we turn to the equaltime 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 singleelectrode 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 singleelectrode 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 spatiotemporal characteristics of the LFP signals, namely, traveling waves of activity of different types and their characteristics.
Traveling waves of different types
Several groups have reported traveling waves on the motor cortex during the preparatory phase of the movement in instructeddelay 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 spatiotemporal 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.
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 oscillationpromoting effect and the wave speed decreases with increasing input amplitude. This leads in the SN model for ${\nu}_{ext}$ in the range $13$ 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 selfgenerated and the oscillations between different modules are wellsynchronized 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 spatiotemporal patterns arise from the opposing effects of synchronization by longrange 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 crosscorrelations 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 spatiotemporal 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 crosscorrelation 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.10.2$ mm^{2} 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 longrange 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 spatiotemporal 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 proxyLFP 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 longrange connectivity, is crucial for the appearance of planar waves. This conclusion is further supported by varying the connectivity range as described below.
Properties of planar waves and characteristics of their inputs
The presence of planar waves is one remarkable feature of the multielectrode 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 nonmyelinated horizontal fibers (GonzálezBurgos 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 longrange 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 longrange 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 crosscorrelation plots in Figure 3—figure supplement 12a. This confirms that longrange 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.$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 halfduration coincide, and, where spacedependent 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 proxyLFP (Figure 4b) clearly show the planar wave propagation. Similar mean spatiotemporal phase maps are obtained by distinguishing three groups of planar wave velocities before averaging (Figure 4—figure supplement 1).
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 selfsustained 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).
Synaptic connection anisotropy and planar wave propagation direction
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 longrange 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.
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 ydirection than in the xdirection, 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 yaxis are clearly different. There is a large predominance of waves observed to propagate along the xaxis, namely, along the axis where the connectivity is shorter range. Longerrange connectivity promotes longerrange synchronization along the yaxis and synchronized states. The weaker connectivity along the xaxis allows for the formation of the larger phase gradients along the xaxis 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 crosscorrelation 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 ydirection (${l}_{y}^{(I)}=1$ or 400 μm in real units) than in the xdirection (${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 xaxis 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 nonhuman primates in the rostrocaudal direction. Rubino et al., 2006 found that planar waves predominantly propagated along the same rostrocaudal 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 mediolateral dimension, as the model suggests.
Discussion
We have proposed and studied here a simplified model of the motor cortex based on local recurrent connections coupling excitatory and inhibitory neuronal populations together with longrange excitatory connections targeting both populations. This builds on previous works describing beta oscillations in motor cortex as arising from recurrent interactions between excitatory and inhibitory populations (Pavlides et al., 2015; Sherman et al., 2016; Chen et al., 2020). Of note, the model of a single EI module of Sherman et al., 2016 includes two layers of hundred multicompartment spiking neurons to approximately take into account the layered structure of the cortex and allow a detailed comparison with MEG signals. In a complementary fashion, our model uses a simpler description of an EI module but goes beyond these previous works by taking into account longrange interactions and analyzing the resulting spatiotemporal patterns.
The wellknown oscillatory instability of coupled EI networks is not qualitatively modified by the presence of longrange connections and the model displays an instability leading to sustained oscillatory behavior. These oscillations occur at beta frequency for adequate synaptic parameters. Comparison with available data has led us to conclude that the motor cortex operates close to this beta oscillatory instability line under strong fluctuating inputs from other areas. This constantly displaces the operating point of the dynamics and leads to waxing and waning of beta amplitude oscillations. A similar proximity of an oscillatory regime has been proposed for the bursts of beta oscillations in the basal ganglia using different models (Rubchinsky et al., 2012; Bahuguna et al., 2020). Analysis of this oscillatory behavior in time and across electrodes has led us to infer the characteristics of the appropriate external inputs. The observed power enhancement on the low frequencyside of the beta peak requires that their time correlation is not too short on a 10 ms scale while the duration of the observed oscillation bursts requires that it is not too long on the same time scale. In an analogous way, the external input amplitude should be, on the one hand, sufficiently large to significantly modify the LFP characteristics above the background of intrinsic local fluctuations arising from the finite number of neurons sampled by each electrode. On the other hand, it should be small enough to preserve the harmonicity of beta oscillations, as indicated by the absence of a significant secondary 50–60 Hz peak in the LFP spectra. The spatial correlation of beta oscillations at the millimeter distance scale has furthermore led us to suggest that the motor cortex receives inputs that target local areas but also synchronous inputs that target it more globally. Under these conditions, traveling waves of different kinds are observed that resemble those recorded in vivo both in their repartition between the different wave types and their speeds of propagation. Characterization of the inputs onto motor cortex during movement preparation is needed to test these theoretical predictions.
Propagating waves of neural activity have been observed in different contexts (Ermentrout and Kleinfeld, 2001; Muller et al., 2018). Some are unrelated to neural oscillations such as subthreshold waves (Bringuier et al., 1999) or propagation of spiking activity in visual cortex (Davis et al., 2020). Others are based on oscillatory activity with proposed mechanisms relying on a gradient of frequencies (Ermentrout et al., 1998) or structural sequences of activation (Muller et al., 2016). The traveling waves and the mechanism here proposed to underlie them are quite different. They are obviously based on the existence of oscillations since the traveling waves are a reflection of different phases of oscillation at different positions on motor cortex. The model we have developed here includes structural connectivity that tends to synchronize oscillations at different positions, but the stochastic external inputs play an essential role in creating dephasings between different locations. The wave propagation speeds are found to be rather low and distributed in the few tens of cm/s. The mean speed of about 30 cm/s is comparable to the propagation speed along nonmyelinated horizontal fibers (GonzálezBurgos et al., 2000; Girard et al., 2001). However, we have found that it is actually independent of it, as a comparable distribution of traveling speeds is found in a model with no propagation delays. This stands, for instance, in sharp contrast with a recent model of the propagation of spiking activity in visual cortex (Davis et al., 2021), which does not involve oscillatory activity but also takes into account strong local fluctuations of neural activity. In the present context, the observed speeds of the oscillatory waves of activity are set by the oscillatory frequency and by the dephasings produced by the external inputs in the recorded area.
We have proposed that external inputs have both local and global components. The existence of a global component appears consistent with the presence of global inhibition in motor cortex during movement preparation (Greenhouse et al., 2015). The model ‘external inputs’ could represent direct inputs from different neural structures, including other cortical areas such as frontal and parietal cortex (Davare et al., 2011), and extracortical ones like the thalamus (Hooks et al., 2013), or a mixture of the two. The known thalamocortical connectivity (Jones, 2001) makes the thalamus a privileged candidate for the origin of, at least, part of these external inputs. Indeed the described diffused connectivity from calbindinpositive matrix neurons could be the source of our global inputs while core parvalbuminpositive neurons could be the source of our local inputs. This requires further experiments assessing the origin of synaptic inputs and their influence on waves and wave types. We have moreover suggested that the observed anisotropy of wave propagation is linked to the anisotropy of local input targets, with the specific prediction that traveling waves of oscillatory activity are predominantly observed orthogonal to the axis where they are more spread out. This also needs to be tested in further experiments.
The model suggests that inputs to the motor cortex play an important role in the generation of beta oscillations and of spatiotemporal patterns during movement preparation. The experimental measurement of these inputs and their correlations at different locations thus appears an important research direction. Do synchronized inputs target distant locations of motor cortex? Does a a shift of these global inputs accompany the birth of planar waves as we have found in the model? Perturbation experiments, for instance, by optogenetic means, would also bring important information. Can sustained beta oscillations be induced by mild inhibition of pyramidal cells and interneurons? The model furthermore suggests that such perturbations applied in a transient way would enhance the appearance of planar waves.
Finally, the role of traveling waves during movement preparation remains an important question that needs to be clarified. We have shown that a stochastic description of the inputs to the motor cortex is sufficient to account for the recordings. However, the produced traveling waves are correlated with the particular inputs that the motor cortex receives. Indeed, traveling waves have been shown to carry information about the subsequent movement (Rubino et al., 2006). Movement preparation is presently viewed, in a dynamical systems perspective, as taking place in a subspace orthogonal to the dynamics of the movement itself and producing the proper initial condition for it (Kaufman et al., 2014; Wang et al., 2018; Zimnik and Churchland, 2021; Inagaki et al., 2022; BachschmidRomano et al., 2022). Beta oscillations and their associated timescale are however absent from this description. Further work is needed to refine it and include beta oscillations and waves and obtain a more comprehensive description of motor cortex dynamics.
Methods
Recording data
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.
Data analysis
An overview of our data analysis protocol is provided in Figure 3—figure supplement 1. The different steps are detailed below.
Data filtering
The index of each electrode in the data (Figure 2—figure supplement 2) was matched to the $(x,y)$position of the electrode in the data that we use in the following. The few missing electrode signals were replaced by the average of the neighboring electrode signals. The followings steps were applied to the so obtained $N=100$ electrode signals. The LFP signal of each electrode was bandpass filtered in the 13–30 Hz range using a thirdorder Butterworth filter (signal.filtfilt function of the Python package scipy). The filtered signal of each electrode was then zscore normalized and Hilberttransformed, using the scipy Hilbert transform function signal.hilbert to obtain its instantaneous amplitude ${A}_{xy}(t)$ and phase ${\mathrm{\Phi}}_{xy}(t)$ (Figure 2—figure supplement 2). The instantaneous maps ${A}_{xy}(t)$ and ${\mathrm{\Phi}}_{xy}(t)$ were used for beta burst and wave pattern detection, respectively.
Beta burst analysis
For each electrode amplitude time trace ${A}_{xy}(t)$, the percentiles of the amplitude distribution were determined (Tinkhauser et al., 2017). The beta burst threshold was set at the $75th$ percentile. A burst onset was defined as a time point at which the analytical amplitude exceeded the threshold and its termination as the earliest following time point at which the amplitude fell below the threshold. The time difference between these two time points was defined as the burst duration. The amplitude of a burst was defined as the average amplitude during the burst duration. The procedure is illustrated in Figure 2—figure supplement 2.
Spatiotemporal pattern analysis and wave classification
The spatiotemporal patterns of oscillatory activity were classified with the help of the phase map ${\mathrm{\Phi}}_{xy}(t)$ following the steps of Denker et al., 2018 with minor modifications. First, the phase gradient map ${\mathrm{\Gamma}}_{xy}(t)$ and its normalized version, the phase directionality map ${\mathrm{\Delta}}_{xy}(t)$, were obtained. Specifically the phase gradient map was obtained by averaging the phase differences between a point and its existing next and nextnearest neighbors along the x and yaxis,
where the vector has been written as a complex number and ${N}_{x}$, ${N}_{y}$ count the actual number of neighbors along x and y, respectively (${N}_{x}+{N}_{y}$ may differ from 8 for electrodes near the array boundaries). Normalization of the gradients to unit length serves to produce the ‘phase directionality map’ (Denker et al., 2018),
The alignment of the unit vectors ${\mathrm{\Delta}}_{xy}(t)$ at time $t$ was quantified using the norm ${\sigma}_{g}$ of their mean, the ‘circular variance of phase directionality’ (Denker et al., 2018), with
with $N$ the total number of electrodes.
The data was classified into successive 1 ms spatial ‘patterns.’ Patterns with wellaligned ${\mathrm{\Delta}}_{xy}(t)$ were classified as planar wave patterns using the criterion ${\sigma}_{g}(t)>0.5$ and provided the first category of the classification.
The remaining patterns were further classified using additional measures. Patterns were next categorized as ‘radial waves’ (or not) by considering critical points (i.e., maxima, minima or saddle points) of the phase map ${\mathrm{\Phi}}_{xy}(t)$ (Rule et al., 2018). A smoothed phase gradient map, the ‘gradient coherence map,’ was defined following Denker et al., 2018,
(Note that the sum is taken over existing neighboring electrodes, the number of which is given by ${N}_{xy}\le 9$.) Minima, maxima, and saddle points were obtained by locating the sign changes in the local gradient coherence map (Perry and Chong, 1987). Maps with exactly one critical point were classified as radial waves, thus providing the second category of our classification (these patterns could be subdivided in further subcategories [Rule et al., 2018] but this was not pursued).
The synchronous patterns were finally identified using the ‘circular variance of phases’ (Denker et al., 2018), ${\sigma}_{p}$, analogous to the wellknown Kuramoto order parameter,
with $N$ the total number of electrodes in the array. Patterns with ${\sigma}_{p}(t)>0.85$, indicating tight synchronization of the oscillations of all electrodes, were classified as ‘synchronized,’ the third classification category. All remaining patterns were assigned to the ‘random’ fourth category.
After all successive patterns were classified, a planar, a radial, or a synchronized wave episode was registered when at least six consecutive frame beared such a label. In other words, planar, radial, or synchronized wave episodes were requested to last for at least 6 ms to be declared as such. All remaining patterns were registered as random waves.
Alignment and averaging of planar waves
To compute average quantities relating to planar waves like average phase maps or input strength, we aligned events in time in two different ways, either (i) with respect to the phase of the underlying oscillation or (ii) with respect to the onset of a planar wave. For the phase alignment, we determined the average phase across the array at half the duration of any given event, which we denote ${\overline{\varphi}}_{\mathrm{center}}^{(k)}$ for event $k$ in the following. We then converted this phase to a time with respect to a reference frequency f_{ref} according to ${t}_{\mathrm{center}}^{(k)}={\overline{\varphi}}_{\mathrm{center}}^{(k)}/(2\pi {f}_{\mathrm{ref}})$, where ${f}_{\mathrm{ref}}={\sum}_{k}{f}_{k}/K$ is the mean frequency over all $K$ events of a given type, with f_{k} being given by $\u27e8d\varphi /dt\u27e9/(2/\pi )$ averaged over space and time during the event. We then aligned all events in time by centering each event $k$ on time ${t}_{\mathrm{center}}^{(k)}$ rounded to the closest time bin. With this choice, time $t=0$ corresponds to an average phase $\overline{\varphi}=0$. For the alignment relative to wave onset, we simply aligned all planar wave events with respect to the first frame classified as planar wave.
To furthermore compute averages over planar waves without loosing spatial information related to wave direction, but not restrict ourselves to the analysis of waves that propagate along a single direction, we considered all wave events that propagate along any of the four principal directions (0, $\pi /2$, $\pi $, $3\pi /2$) of the grid, which are statistically equivalent. We specifically selected waves with a mean gradient vector, averaged over units and wave duration, that points within a $\pm \pi /8$ interval centered on the respective direction. (Note that waves with a mean gradient outside the selected sector are not taken into account.) In a second step, we then rotated units by (0, $3\pi /2$, $\pi $, $\pi /2$) for events of the four principal directions, respectively, such that eventually all gradients pointed to the $\pm \pi /8$ interval centered on 0, corresponding to a leftward traveling wave (opposite to the mean gradient). As all four principal directions are equivalent due to the symmetry of the system, this allowed us to improve the statistics of our measurements without introducing any bias because of the discretized grid.
Model
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 EI modules are connected by longrange excitatory connections with a connection probability decaying with the distance between modules. The neural activity of each EI module is represented at the level of its excitatory and inhibitory neuronal populations in the ratemodel 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 integrateandfire (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 EI 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 fI 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 twodimensional network of EI modules coupled by longrange 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 $\u27e8{\xi}_{A}(\mathbf{x},t){\xi}_{B}({\mathbf{x}}^{\prime},{t}^{\prime})\u27e9={\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 EI 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 timedependent, targeting the excitatory (resp. inhibitory) population of the EI module at position $\mathbf{x}$,
With our normalization choice $\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 OU 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, $\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 postsynaptic 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 nonmyelinated fibers in the motor cortex (GonzálezBurgos 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}$.
Theoretical analysis
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 EI 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.
Steady state
We consider first the steady state of the deterministic network in which excitatory populations and inhibitory populations are firing at constant rates ${r}_{E}^{s}$ and ${r}_{I}^{s}$, independently of time and module position, with
In this state, the synaptic currents in the different populations are obviously also independent of time and position. Given the normalization conditions (Equations 15; 18), they simply read,
From the model definition (Equations 8; 9), these firing rates are produced by constant external inputs of magnitudes,
Bifurcation lines
The stability of the steady firing state can be assessed by imposing the constant external currents (Equation 25) and computing the dynamics of small perturbations around the steady state. To this end, we linearize Equations 8–17 and look for solutions that are oscillatory in space and exponential in time,
with similar expansions for the other variables. Substitution in the explicit formulas (Equation 14), or in the corresponding differential equations (Equations 16; 17), provides the expression of the synaptic current in term of the module activities,
where ${\stackrel{~}{S}}_{E}(\sigma )$ and ${\stackrel{~}{S}}_{I}(\sigma )$ are the Laplace transforms of ${S}_{E}(t)$ and ${S}_{I}(t)$,
Substitution of Equation 27 in the linearized Equations 8; 9 gives
The matrix ${\stackrel{~}{L}}_{EI}(\mathbf{q},\sigma )$ reads,
where the function $C(\mathbf{q},\sigma )$ is the Fourier transform of the coupling function with propagation delays taken into account,
In Equation 30 and in the following, the shorthand notation ${\tau}_{E}$ (resp. ${\tau}_{I}$) is used for ${\tau}_{E}({I}_{E}^{s})$ (resp. ${\tau}_{I}({I}_{I}^{s})$). The existence of a nontrivial solution of Equation 29 requires that the determinant, $W(\mathbf{q},\sigma )$, of the matrix ${\stackrel{~}{L}}_{EI}(\mathbf{q},\sigma )$ vanishes, that is,
with the functions ${\stackrel{~}{T}}_{E}(\sigma ),{\stackrel{~}{T}}_{I}(\sigma )$, defined by,
The constants $\alpha $ and $\gamma $ respectively measure the gain of monosynaptic recurrent excitation and inhibition while $\beta $ measures the gain of disynaptic recurrent inhibition,
where ${\mathrm{\Phi}}_{A}^{\prime}(I)$ denotes the derivative of ${\mathrm{\Phi}}_{A}$ with respect to $I$.
The oscillatory instability, or ‘Hopf bifurcation,’ line corresponds to the parameters for which the growth rate is purely imaginary, $\sigma =i\omega $. It is obtained in parametric form, with $\alpha $ and $\beta $ as functions of the frequency $\omega $ (and of the recurrent inhibition $\gamma $) by separating the real and imaginary parts of Equation 32 and solving the resulting linear equations for $\alpha $ and $\beta $. This gives
The instability first appears at long wavelengths, namely, at $\mathbf{q}=0$ on a large enough lattice. The expressions (35) with $\mathbf{q}=0$ have been used to draw the oscillatory instability lines and the frequency dependence on parameters shown in Figure 1b and Figure 1—figure supplement 4.
Besides this oscillatory instability, there is a possible loss of stability towards a highfiring rate state. It is obtained when the growth rate $\sigma $ changes with parameter variation from being real negative to real positive, that is, for parameters such that $W(\mathbf{q},0)=0$. Since at the ‘real’ instability threshold, the instability growth rate vanishes, the instability line is independent of the synaptic current kinetics. One indeed checks from Equations 28; 33 that ${\stackrel{~}{T}}_{E}(0)=1$. Thus, Equation 32 gives for the real instability line for a given wavevector $\mathbf{q}$,
Since $C(\mathbf{q},0)<1$ when $\mathbf{q}$ does not vanish, the instability appears first at $\mathbf{q}=0$ with $C(0,0)=1$ when all the modules are in the exact same state.
The expressions (Equation 35, Equation 36) with $\mathbf{q}=0$ have been used to draw the diagram of Figure 1b. The diagram of Figure 1f is similarly obtained by computing the external inputs corresponding to different steady discharges, for fixed synaptic parameters, and assessing the stability of these states from their position with respect to the stability lines (Equation 35 and Equation 36). We have taken the synaptic strength of external inputs such that the network operating point crosses the oscillatory bifurcation line when the external input amplitude varies:
Auto and crosscorrelations of module activities and power spectrum
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 ${\stackrel{~}{L}}_{EI}(\mathbf{q},i\omega )$ is given in Equation 30. The two components of the stochastic forcing term $\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 $W(\mathbf{q},i\omega )$ defined in Equation 32 and ${V}_{E}(\mathbf{q},i\omega )$ given by,
The Fourier components of the input fluctuations read
with the white noise averages,
The shorthand notation ${\delta}^{2}(\mathbf{q})$ has been used for the twodimensional $\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 crosscorrelation is a real function, as it should, since ${S}_{EE}^{ext}$ and ${S}_{EE}^{N}$ are related to their complex conjugates by ${S}_{EE}^{ext\phantom{\rule{thinmathspace}{0ex}}\ast}(\mathbf{q},\omega )={S}_{EE}^{ext}(\mathbf{q},\omega )$ and ${S}_{EE}^{N\phantom{\rule{thinmathspace}{0ex}}\ast}(\mathbf{q},\omega )={S}_{EE}^{N}(\mathbf{q},\omega )$, a symmetry inherited from the function $C(\mathbf{q},\omega )$ (Equation 31).
Taking coincident points (i.e., $\mathbf{x}={\mathbf{x}}^{\mathrm{\prime}}$) in the current–current correlation (Equation 50) gives the ${I}_{E}$ current autocorrelation for a local module in the network. Remembering that the autocorrelation is the Fourier transform of the spectrum, this provides the spectrum $S(\omega )$ of a local module excitatory current time series:
We take the local module excitatory current as a proxy for the LFP. Equations 50; 51 have been used to draw the theoretical lines for current–current correlations and power spectra in Figure 2, Figure 2—figure supplement 1 and Figure 2—figure supplement 3.
Simulations
The mathematical model defined in section ‘Model’ was numerically simulated. We used a 28 × 28 grid of EI modules. Two external layers of with EI populations at their fixed points were added as boundary conditions. Measurements were only performed in the center square 10 × 10 grid to minimize boundary effects. A sketch of the grid is shown in Figure 1—figure supplement 2. Model distributions and averages in all figures were obtained by performing 20 independent network simulations of 10 s simulated time each.
The reference parameters for all simulations are provided in Table 1.
Simulations were performed with a custom C program using a firstorder Euler–Maruyama integration method, with a time step $dt=0.01$ ms. Python programs were used for data analysis and to draw the figures. These computer codes are available at https://github.com/LKANG777/BetaOscillation (copy archived at Kang, 2023). All simulations were performed on ECNU computer clusters.
Data availability
The source codes for this manuscript are available on GitHub at https://github.com/LKANG777/BetaOscillation (copy archived at Kang, 2023).

GNodeMassively parallel recordings in macaque motor cortex during an instructed delayed reachtograsp task.https://doi.org/10.12751/gnode.f83565
References

Mathematical frameworks for oscillatory network dynamics in neuroscienceJournal of Mathematical Neuroscience 6:2.https://doi.org/10.1186/s1340801500336

Uncoupling the roles of firing rates and spike bursts in shaping the stngpe beta band oscillationsPLOS Computational Biology 16:e1007748.https://doi.org/10.1371/journal.pcbi.1007748

Über das elektrenkephalogramm des menschenArchiv Für Psychiatrie Und Nervenkrankheiten 87:527–570.https://doi.org/10.1007/BF01797193

BookRhythms of the brainOxford University Press.https://doi.org/10.1093/acprof:oso/9780195301069.001.0001

On the nature of the intrinsic connectivity of the cat motor cortex: Evidence for a recurrent neural network topologyJournal of Neurophysiology 102:2131–2141.https://doi.org/10.1152/jn.91319.2008

Emergence of beta oscillations of a resonance model for Parkinson’s diseaseNeural Plasticity 2020:8824760.https://doi.org/10.1155/2020/8824760

On the complexity of resting state spiking activity in monkey motor cortexCerebral Cortex Communications 2:tgab033.https://doi.org/10.1093/texcom/tgab033

Interactions between areas of the cortical grasping networkCurrent Opinion in Neurobiology 21:565–570.https://doi.org/10.1016/j.conb.2011.05.021

BookTheoretical Neuroscience: Computational and Mathematical Modeling of Neural SystemsMIT press.

Multiple pulse interactions and averaging in systems of coupled neural oscillatorsJournal of Mathematical Biology 29:195–217.https://doi.org/10.1007/BF00160535

Minimal model of oscillations and waves in the Limax olfactory lobe with tests of the model’s predictive powerJournal of Neurophysiology 79:2677–2689.https://doi.org/10.1152/jn.1998.79.5.2677

BookMathematical foundations of neuroscienceNew York, NY: Springer Science & Business Media.https://doi.org/10.1007/9780387877082

Contributions of intrinsic membrane dynamics to fast network oscillations with irregular neuronal dischargesJournal of Neurophysiology 94:4344–4361.https://doi.org/10.1152/jn.00510.2004

Feedforward and feedback connections between areas V1 and V2 of the monkey have similar rapid conduction velocitiesJournal of Neurophysiology 85:1328–1331.https://doi.org/10.1152/jn.2001.85.3.1328

Nonspecific inhibition of the motor system during response preparationThe Journal of Neuroscience 35:10675–10684.https://doi.org/10.1523/JNEUROSCI.143615.2015

Mapping horizontal spread of activity in monkey motor cortex using single pulse microstimulationFrontiers in Neural Circuits 10:104.https://doi.org/10.3389/fncir.2016.00104

Organization of cortical and thalamic input to pyramidal neurons in mouse motor cortexThe Journal of Neuroscience 33:748–760.https://doi.org/10.1523/JNEUROSCI.433812.2013

Electrocorticograms in man: Effect of voluntary movement upon the electrical activity of the precentral gyrusArchiv Für Psychiatrie Und Nervenkrankheiten 183:163–174.https://doi.org/10.1007/BF01062488

The thalamic matrix and thalamocortical synchronyTrends in Neurosciences 24:595–601.https://doi.org/10.1016/s01662236(00)019226

SoftwareBetaOscillation, version swh:1:rev:0cc496ef13b5db7a51b3a1d2c7b06403916bacb7Software Heritage.

Cortical activity in the null space: Permitting preparation without movementNature Neuroscience 17:440–448.https://doi.org/10.1038/nn.3643

Synchronization, stochasticity, and phase waves in neuronal networks with spatiallystructured connectivityFrontiers in Computational Neuroscience 14:569644.https://doi.org/10.3389/fncom.2020.569644

LongRange horizontal connections between supragranular pyramidal cells in the extrastriate visual cortex of the ratThe Journal of Comparative Neurology 344:543–558.https://doi.org/10.1002/cne.903440405

Cortical travelling waves: Mechanisms and computational principlesNature Reviews. Neuroscience 19:255–268.https://doi.org/10.1038/nrn.2018.20

From spiking neuron models to linearnonlinear modelsPLOS Computational Biology 7:e1001056.https://doi.org/10.1371/journal.pcbi.1001056

A description of eddying motions and flow patterns using criticalpoint conceptsAnnual Review of Fluid Mechanics 19:125–155.https://doi.org/10.1146/annurev.fl.19.010187.001013

Natural frequencies of human corticothalamic circuitsThe Journal of Neuroscience 29:7679–7685.https://doi.org/10.1523/JNEUROSCI.044509.2009

Intermittent neural synchronization in Parkinson’s diseaseNonlinear Dynamics 68:329–346.https://doi.org/10.1007/s110710110223z

Propagating waves mediate information transfer in the motor cortexNature Neuroscience 9:1549–1557.https://doi.org/10.1038/nn1802

Paradoxical effects of external modulation of inhibitory interneuronsThe Journal of Neuroscience 17:4382–4388.https://doi.org/10.1523/JNEUROSCI.171104382.1997

Neurophysiological and computational principles of cortical rhythms in cognitionPhysiological Reviews 90:1195–1268.https://doi.org/10.1152/physrev.00035.2008

Flexible timing by temporal scaling of cortical responsesNature Neuroscience 21:102–110.https://doi.org/10.1038/s4159301700286

Independent generation of sequence elements by motor cortexNature Neuroscience 24:412–424.https://doi.org/10.1038/s41593021007985
Decision letter

Nicole C SwannReviewing Editor; University of Oregon, United States

Laura L ColginSenior Editor; University of Texas at Austin, United States

Thomas KnöpfelReviewer; Imperial College London, United Kingdom

Nicholas G HatsopoulosReviewer; University of Chicago, United States
Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Decision letter after peer review:
Thank you for submitting your article "Βeta oscillations and waves in motor cortex can be accounted for by the interplay of spatiallystructured connectivity and fluctuating inputs" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Laura Colgin as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Thomas Knopfel (Reviewer #1); Nicholas G Hatsopoulos (Reviewer #3).
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential revisions:
(1) When introducing the multielectrode arrays as well as when describing the model layout, spatial dimensions should be stated upfront.
(2) The abbreviations for the wave types (Pla., Syn., Rad., Ran. ) need to be expanded and introduced appropriately.
(3) Regarding the "external inputs": the design of the model would allow for interpreting the external inputs as inputs from other cortical areas or a mixed corticalextracortical input. Similarly, the global and nonglobal inputs may have different mixtures of cortical and extracortical components. This should be stated explicitly.
(4) The authors repetitively state a relatively slow propagation speed of 0.3 m/s along unmyelinated horizontal axons (with no supporting reference). This value is correct for horizontal cerebellar parallel fibers (which are 0.1 μm in diameter) but to my knowledge not for the thicker horizontal longrange axons in the cerebral cortex. This should be discussed/modified.
(5) Could you show the spatiotemporal structure of an average of traveling waves (average across same type and perhaps a limited range of feature parameters)?
(6) Anatomical data on connectivity (type, spatial decay and anisotropies) should be derived from literature and these empirical data should be used to consider the plausibility of the corresponding model assumptions.
(7) Can you derive tractable tests to evaluate the proposed model against alternative published models and can you propose new tractable experimental research questions based on the insights from your study?
(8) It would be helpful to have an intuitive explanation on what are the parameters affecting the β frequency of oscillations in the main text.
(9) In the Methods, the authors state that they consider the local module of excitatory current as a proxy for the LFP; this should be stated in the main text too, for clarity.
(10) The authors mention that reference (Davis et al., 2021) "also relies on external inputs". Isn't that model without an external drive, only relying on selfsustained dynamics?
(11) On lines 325327, the authors state that the SN model leads to a proportion of planar wave states of a few percent. Figures 3l shows that the proportion of planar waves is actually quite large, >60%.
Reviewer #2 (Recommendations for the authors):
Related work
l 178179: "These networks that operate close to the bifurcation line thus appear as promising candidates to describe waxingandwaning of β oscillations seen in recording data"
The authors may be interested to know that similar analysis/observations have been made in the EI subcircuit of the basal ganglia (STNGPe) loop by two studies maybe indicating a common mechanism of β bursts across cortical and subcortical circuits.
 Rubchinsky LL, Park C, Worth RM. Intermittent neural synchronization in Parkinson's disease. Nonlinear Dynamics. 2012; p. 329346. pmid:22582010 –
 Bahuguna J, Sahasranamam A, Kumar A (2020) Uncoupling the roles of firing rates and spike bursts in shaping the STNGPe β band oscillations. PLOS Computational Biology 16(3): e1007748. https://doi.org/10.1371/journal.pcbi.1007748
https://doi.org/10.7554/eLife.81446.sa1Author response
Essential revisions:
(1) When introducing the multielectrode arrays as well as when describing the model layout, spatial dimensions should be stated upfront.
The dimensions of the multielectrode array and the simulation grid are now described at the end of the Introduction (lines 103104) and in the first paragraph of the Results section (lines 119120). We also refer there to Figure 1Figure supp.2 depicting the simulation grid.
(2) The abbreviations for the wave types (Pla., Syn., Rad., Ran. ) need to be expanded and introduced appropriately.
We have defined the abbreviations in the text (lines 325326) and in the caption of Figure 3 where they appear for the first time.
(3) Regarding the "external inputs": the design of the model would allow for interpreting the external inputs as inputs from other cortical areas or a mixed corticalextracortical input. Similarly, the global and nonglobal inputs may have different mixtures of cortical and extracortical components. This should be stated explicitly.
We have now stated this explicitly in the Results first section (lines 245246) and in the discussion (lines 610).
(4) The authors repetitively state a relatively slow propagation speed of 0.3 m/s along unmyelinated horizontal axons (with no supporting reference). This value is correct for horizontal cerebellar parallel fibers (which are 0.1 μm in diameter) but to my knowledge not for the thicker horizontal longrange axons in the cerebral cortex. This should be discussed/modified.
We apologize for having omitted to give references supporting our choice of conduction velocity. The conduction velocity along horizontal unmyelinated axons have been measured in different works with results in the range 0.10.6m/s. Most of the works we know in the cortex of monkeys are performed in visual cortex. We actually based our choice on Girard, Hup´e & Bullier, J Neurophysiol (2001) who find a median conduction velocity 0.3 m/s along horizontal unmyelinated axon at physiological temperature. In Guillermo GonzalezBurgos et al., Cerebral Cortex (2000) the mean horizontal axonal conduction velocity was measured in monkey prefrontal cortex to be 0.14 m/s at room temperature (20^{◦} − 24^{◦}C). It seems reasonable to attribute the difference, as the authors do, to the difference in experimental temperature.
Using voltage sensitive dyes, Grinvald et al., J Neurosci (1994) found the speed of propagation along lateral connections (“apparent conduction velocity” which may underestimate the true conduction velocity) to be in the range of 0.09−0.25 m/s. Finally, in the rat visual cortex, Lohmann & Rorig, J Comp Neurol (1994) report a mean conduction velocity of 0.28 m/s and in rat motor cortex, Aronadiou & Keller, J Physiol (1993) report a conduction velocity of 0.1 m/s. However, it is unclear to us whether data obtained in rats could be directly applied to monkey motor cortex. Given these data, it seems that a choice of conduction velocity of 0.3 m/s is reasonable and not an underestimate. We have now cited (lines 126127) the references pertaining to monkey cortex which seem the most relevant for our study. Nonetheless, we are aware of the possibility that the conduction velocity could be larger in monkey motor cortex. That is the reason why we investigated the other extreme case in which propagation delays are actually negligible. This leads us to conclude that the precise value of the velocity does not impact our conclusions.
(5) Could you show the spatiotemporal structure of an average of traveling waves (average across same type and perhaps a limited range of feature parameters)?
In order to address this suggestion and related questions by the referees, we have aligned different planar wave events using their computed spatiotemporal phase map, both for simulations and monkey data. With the help of this alignment, we have computed the average over events of the LFP (or its proxy in simulations) and of the synchronization parameter σ_{p}. Additionally for the simulations, we have computed the averages of the local inputs as well as of the global inputs, which are not available in the recording data. The structure of the planar traveling waves clearly appear in the average phase and (proxy) LFP maps. We see that the synchronization parameter rises concomitantly with the appearance of the planar waves both in simulations and in experiments. Remarkably, in simulations, this corresponds to an excursion of the global input away from its mean value, which persists during the wave episode. We did not observe any significant average features in the local inputs when averaged over events (not shown), presumably because there is a large variability of local inputs that can contribute to the appearance of a planar wave. The results of this analysis are described in the new section “Properties of planar waves an characterization of their inputs” and shown in the new Figure 4 and Figure 4Figure supp. 1, and in Figure 4Figure supp. 2 for the experimental recordings.
In order to try and go beyond this purely correlative observations, we have also performed simulations in which transient steps of global inputs have been injected, with a magnitude similar to the observed excursion. This indeed results in a significant increase of planar wave episodes, as shown in Figure 4eg.
(6) Anatomical data on connectivity (type, spatial decay and anisotropies) should be derived from literature and these empirical data should be used to consider the plausibility of the corresponding model assumptions.
Precise anatomical data on connectivity in monkey motor cortex is not abundant. We base our choices (lines 122125) on Gatter and Powell, Brain (1978), Huntley and Jones, J Neurophys (1991), Hao et al., Front Neural Circ (2016) and on the data of Capaday et al., J Neurophys (2009) on cat motor cortex. We thank referee 3 for pointing out the data on anisotropic connectivity in Gatter and Powell, Brain (1978) that we now take into account and discuss in the revised section “Synaptic connection anisotropy and planar wave propagation direction”.
(7) Can you derive tractable tests to evaluate the proposed model against alternative published models and can you propose new tractable experimental research questions based on the insights from your study?
This is of course an important question. We are not aware of other models specifically aimed at describing waves during episodes of β oscillations in the motor cortex. More generally, several models have been developed to describe the propagation of neural activity but among those, few are targeted to traveling waves of oscillatory activity, namely oscillations at different spatial location with a dephasing that depends on space. The ones we know and cite are based on a gradient of excitability/natural oscillation frequency, like for instance the study of Ermentrout et al., J Neurophys (1998) aimed at describing traveling oscillatory activity in the Limax olfactory lobe, or on inputs targeting different locations with different structural delays, as the model developed by Muller et al., eLife (2016) to describe rotating waves during sleep spindles. That the oscillatory waves of β oscillations in motor cortex are linked to such structural features appears unlikely to us since propagation appears to take place in different directions and particularly along the two opposite directions of the same axis. This led us to propose the model in the manuscript. Since the model relies on inputs from other cortical or subcortical areas, measuring these would offer important tests of the model. We assume that there are global synchronized inputs to the motor cortex. Showing that this is indeed the case would be a first important test of the model. We now show in the revised manuscript (new Figure 4) that shortly before the appearance of a planar wave, there is a shift in the strength of the global inputs. Experimental observation of this effect would bring further support to our model.
Optogenetic perturbations could provide further tests. Our model predicts that changing the strength of the global inputs would allow one to put the motor cortex in a regime of sustained oscillations at β frequency. This could in principle be tested by retroviral injections and optogenetic techniques. As described in point (5) above, we also show in the revised manuscript that creating a diminution of the global external inputs similar to the spontaneous decrease observed in our simulations at the onset of planar waves, increases the probability of appearance of a planar wave (new Figure 4). This could also be experimentally attempted. Finally, on the anatomical side the model suggests that local inputs preferentially target the cortex in the mediolateral direction.
We have added a paragraph to the Discussion section (lines 619627) to discuss these experimental research questions suggested by our study.
(8) It would be helpful to have an intuitive explanation on what are the parameters affecting the β frequency of oscillations in the main text.
We have added a paragraph (lines 152162) to explain that increasing recurrent excitation slows down the oscillation, that increasing recurrent inhibition among inhibitory interneurones increases the oscillation frequency and to discuss the influence of the kinetics of the synaptic current. We have added Figure 1Figure supp. 4 to show how these parameters affect the oscillation frequency and change.
(9) In the Methods, the authors state that they consider the local module of excitatory current as a proxy for the LFP; this should be stated in the main text too, for clarity.
This has now been stated in the main text (lines 230232). We now refer to the model signal as the proxyLFP to avoid any misunderstanding.
(10) The authors mention that reference (Davis et al., 2021) "also relies on external inputs". Isn't that model without an external drive, only relying on selfsustained dynamics?
This has been corrected. We are sorry for the misformulation. We meant to write “also take into account strong local fluctuations of neural activity”.
(11) On lines 325327, the authors state that the SN model leads to a proportion of planar wave states of a few percent. Figures 3l shows that the proportion of planar waves is actually quite large, >60%.
In Figure 3l, the proportion of planar corresponds to the green sector, it is equal to 2.7% which we indeed describe as a few percent in the text. The 60 % blue sector corresponds to the proportion of synchronized states. The color code is indicated in the insert. The abbreviations (Pla, Syn,..) have now been defined in the caption as well as in the main text (lines 325326).
https://doi.org/10.7554/eLife.81446.sa2Article and author information
Author details
Funding
China Scholarship Council (Graduate Student Fellowship)
 Ling Kang
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We are very grateful to M Denker for extensive discussions about the analysis in Denker et al., 2018 and for allowing us to precisely compare it with our own. It is a pleasure to thank Nicolas Brunel and David Hansel for thoughtprovoking comments, Nicolas Brunel and Ludovica BachschmidRomano for sharing their unpublished work on the motor cortex and Betsy Herbert for her work and discussions during an internship on a related project. We also wish to acknowledge exchanges with Thomas Brochier and Alexa Riehle, and thank Nicholas Hatsopoulos for pointing out to us the data on anisotropic connectivity in Gatter et al., 1978.
Senior Editor
 Laura L Colgin, University of Texas at Austin, United States
Reviewing Editor
 Nicole C Swann, University of Oregon, United States
Reviewers
 Thomas Knöpfel, Imperial College London, United Kingdom
 Nicholas G Hatsopoulos, University of Chicago, United States
Version history
 Preprint posted: June 17, 2022 (view preprint)
 Received: June 27, 2022
 Accepted: March 2, 2023
 Accepted Manuscript published: March 14, 2023 (version 1)
 Version of Record published: April 18, 2023 (version 2)
Copyright
© 2023, Kang et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 795
 Page views

 128
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Neuroscience
The functional complementarity of the vestibuloocular reflex (VOR) and optokinetic reflex (OKR) allows for optimal combined gaze stabilization responses (CGR) in light. While sensory substitution has been reported following complete vestibular loss, the capacity of the central vestibular system to compensate for partial peripheral vestibular loss remains to be determined. Here, we first demonstrate the efficacy of a 6week subchronic ototoxic protocol in inducing transient and partial vestibular loss which equally affects the canal and otolithdependent VORs. Immunostaining of hair cells in the vestibular sensory epithelia revealed that organspecific alteration of type I, but not type II, hair cells correlates with functional impairments. The decrease in VOR performance is paralleled with an increase in the gain of the OKR occurring in a specific range of frequencies where VOR normally dominates gaze stabilization, compatible with a sensory substitution process. Comparison of unimodal OKR or VOR versus bimodal CGR revealed that visuovestibular interactions remain reduced despite a significant recovery in the VOR. Modeling and sweepbased analysis revealed that the differential capacity to optimally combine OKR and VOR correlates with the reproducibility of the VOR responses. Overall, these results shed light on the multisensory reweighting occurring in pathologies with fluctuating peripheral vestibular malfunction.

 Neuroscience
Genuinely new discovery transcends existing knowledge. Despite this, many analyses in systems neuroscience neglect to test new speculative hypotheses against benchmark empirical facts. Some of these analyses inadvertently use circular reasoning to present existing knowledge as new discovery. Here, I discuss that this problem can confound key results and estimate that it has affected more than three thousand studies in network neuroscience over the last decade. I suggest that future studies can reduce this problem by limiting the use of speculative evidence, integrating existing knowledge into benchmark models, and rigorously testing proposed discoveries against these models. I conclude with a summary of practical challenges and recommendations.