A mean field model for movement induced changes in the beta rhythm

In electrophysiological recordings of the brain, the transition from high amplitude to low amplitude signals are most likely caused by a change in the synchrony of underlying neuronal population firing patterns. Classic examples of such modulations are the strong stimulus-related oscillatory phenomena known as the movement related beta decrease (MRBD) and post-movement beta rebound (PMBR). A sharp decrease in neural oscillatory power is observed during movement (MRBD) followed by an increase above baseline on movement cessation (PMBR). MRBD and PMBR represent important neuroscientific phenomena which have been shown to have clinical relevance. Here, we present a parsimonious model for the dynamics of synchrony within a synaptically coupled spiking network that is able to replicate a human MEG power spectrogram showing the evolution from MRBD to PMBR. Importantly, the high-dimensional spiking model has an exact mean field description in terms of four ordinary differential equations that allows considerable insight to be obtained into the cause of the experimentally observed time-lag from movement termination to the onset of PMBR (∼ 0.5 s), as well as the subsequent long duration of PMBR (∼ 1 − 10 s). Our model represents the first to predict these commonly observed and robust phenomena and represents a key step in their understanding, in health and disease.


Introduction
The modelling of brain rhythms is now a well established and vibrant part of computational neuroscience. Ever since the first recordings of the human electroencephalogram (EEG) in 1924 by Berger (1929) electrophysiological brain recordings have been shown to be dominated by oscillations (rhythmic activity in cell assemblies) across a wide range of temporal scales and scientists have sought to develop large scale models to describe the five main frequency bands of delta (1 − 4 Hz), theta (4 − 8 Hz), alpha (8 − 13 Hz), beta (13 − 30 Hz) and gamma . Moreover, it has long been known, since the early works of Andrews (1936, 1938), that different brain rhythms can be localised to specific areas of the brain, and that these rhythms can be functionally distinct. For example, they showed that the beta rhythm present in the vicinity of the central sulcus was not affected by the presentation of a weak visual stimulus which suppressed the alpha rhythm, recorded from the occipital lobe. Given the challenge of modelling such complex behaviour it is perhaps no surprise that the long and industrious history of brain modelling has delivered more than one tool for this job. For issues that relate to spike times and their synchrony we can appeal to conductance based modelling, large scale network simulations and theories for understanding coupled oscillators, as recently surveyed in Ashwin et al. (2016). For questions that relate to understanding the coarse grained activity of either synaptic currents, mean membrane potentials or population firing rates, it is more natural to appeal to neural mass models, as reviewed in Coombes (2010). Indeed the latter have proven especially fruitful in providing large scale descriptions of how neural activity evolves over both space as well as time (Coombes et al. 2014;Pinotsis et al. 2014). However, these two approaches are dangerously close to creating a dichotomy so that there is no ideal computational modelling framework for understanding the role of spike-timing in generating localised brain rhythms.
A case in point that challenges the modelling tools currently available to us is the work of Jasper and Penfield (1949) who showed that beta rhythms generated from the motor cortex are suppressed during voluntary movement. This phenomenon is known as movement related beta decrease (MRBD). It wasn't until some years later that the post-movement beta rebound (PMBR) (a temporary rise in amplitude of beta oscillations following movement cessation) was discovered (Riehle and Vaadia 2004;Pfurtscheller et al. 1996;Jurkiewicz et al. 2006). MRBD usually lasts for approximately 0.5 seconds and PMBR can last for up to several seconds. MRBD and PMBR are extremely robust, with clear amplitude changes in individual subjects and trials (Pfurtscheller and Lopes da Silva 1999). Interestingly, similar effects have been seen in studies where the subject is asked to think about moving, without carrying out the movement (Schnitzler et al. 1997;Pfurtscheller et al. 2005). These beta band modulations are believed to be caused by changes of synchrony within a relatively localised region of motor cortex (Stancák and Pfurtscheller 1995). Hence, MRBD is regarded as a special case of eventrelated desynchronisation (ERD) and PMBR a special case of event-related synchronisation (ERS).
Multiple papers have employed a large number of carefully controlled paradigms, in humans and animals to further investigate beta rebound phenomena and their modulations by tasks (see Cheyne 2013;Kilavik et al. 2013 for reviews). However, despite the robust nature of the beta task induced decrease and post stimulus rebound, the effect itself is relatively poorly understood and, at the time of writing, there has been, to our knowledge, no computational model capable of describing the beta rebound. In general, high amplitude beta oscillations are thought to reflect inhibition (Cassim et al. 2001;Gaetz et al. 2011), a hypothesis supported by quantifiable relationships between beta amplitude and local concentrations of the inhibitory neurotransmitter gamma aminobutyric acid (GABA) (Gaetz et al. 2011;Hall et al. 2011;Jensen et al. 2005;Muthukumaraswamy et al. 2013). This means that the observed MRBD might reflect an increase in processing during movement planning and execution and the PMBR might reflect active inhibition of neuronal networks post movement (Alegre et al. 2008;Solis-Escalante et al. 2012). An alternative, but not mutually exclusive hypothesis which has been proposed by Donner and Siegel (2011) (also outlined in Liddle et al. 2016) is that the beta signal, in part, represents long range integration across multiple brain regions (see also Liddle et al. 2013). Indeed this is a hypothesis supported by some evidence suggesting that large scale distributed network connectivity is mediated by beta oscillations Hall et al. 2014;Hipp et al. 2012).
To describe beta rebound we are faced with modelling a mesoscopic brain scale and in particular the changes of synchrony within a population of say 10 6−7 excitatory pyramidal cells and their associated inhibitory interneurons. A neural mass model would be ideal for this scale, if the question of interest related to rate rather than spike, which suggests instead a simulation of a spiking neural network model. Unfortunately the latter can be notoriously hard to gain insight from for very large numbers of neurons. Ideally we would have access to a statistical neurodynamics providing a bridge between the two levels of description. This is an open mathematical problem. However, recent progress in obtaining a mean field reduction for a very specific choice of large scale spiking model has been made, and is ideally suited as a basis for breaking the dichotomy noted above. The single neuron model of choice being either a θ -neuron (So et al. 2014;Luke et al. 2013) or a (formally equivalent) quadratic integrate-and-fire (QIF) neuron (Pazó and Montbrió 1009), and the coupling being global and mediated by pulses (namely instantaneous synapses). Given the dense connections of connections in cortex on small scales (Klinshov et al. 2014) the global coupling assumption is not so restrictive for our purposes, though the assumption of fast synapses should be relaxed to incorporate more realistic post synaptic responses. This is precisely the issue we address here to develop a model capable of explaining MRBD and PMBR.
In what follows, Section 2 gives a recapitulation of cortical rebound, illustrated with newly acquired magnetoencephalography (MEG) data, along with some recently published results. A candidate large scale computational model is described in Section 3, utilising realistic synaptic conductance changes. The model is cast in both voltage and phase variable so that it can be understood both as a QIF network, and also as a phase-oscillator network so that a connection to Kuramoto type networks can be made. Importantly we develop an exact low dimensional mean field description capable of capturing the behaviour of a globally coupled network in the limit of a large number of neurons. The macroscopic variables of interest are now firing rate and mean membrane potential (for the voltage description) or the Kuramoto measure of synchrony (for the phase description). Importantly we show in Section 4 that the response of the mean field model to stimulation leads to spectrograms with all of the key features observed in MRBD and PMBR. Finally in Section 5 we emphasise the benefits of this new type of neural mass model, capable of tracking not only changes in firing rate but also coherence within a population, in describing cortical rebound, as well as discuss natural extensions to our initial single population approach.

MRBD and PMBR: a recapitulation
Beta band modulation is robust across subjects; occurring during internally and externally cued movements, as well as during somatosensory stimulation. To demonstrate this, for somatosensory stimulation, we carried out a series of median nerve stimulation experiments on two healthy participants. The participant's median nerve was stimulated at the wrist using a constant current stimulator and the neurological response was measured using MEG (for more details on the experimental design see Appendix A). The experiment was repeated on separate days to determine how reproducible the results were. Figure 1 shows the relative time-frequency spectrograms for two such experiments for each participant, where baseline activity has been subtracted. The top line represents the data for participant 1 and the bottom line represents participant 2. For each case there is a 10-20% decrease in power for ∼ 0.5 s, demonstrating MRBD. At ∼ 0.5 s there is a 60-100% increase in power, exemplifying PMBR. Although the comparison between participants shows dissimilarities in the shape and length of PMBR, the strength and timing of both the MRBD and PMBR are comparable. Importantly, the similarity between each participant's time-frequency spectrogram on separate days is compelling.
Recent work, reviewed in Brookes et al. (2011Brookes et al. ( , 2012 and Robson et al. (2015), has begun to show the potential importance of beta band modulation. For example, Fig. 2 (recreated with permission from Robson et al. 2015) shows relative time-frequency spectra depicting the changes in neural oscillations in sensorimotor cortex in response to a cued finger movement task. The left hand panel (a) shows the case for healthy individuals. The time-frequency spectrum is extracted from a location of interest in left primary motor cortex. Notice that in the beta band, the MRBD and PMBR are observed clearly. The right hand panel (b) shows the case for patients with schizophrenia. Note the  Robson et al. 2015.) significant reduction in PMBR. Furthermore, this same study showed that the magnitude of the beta rebound correlated significantly with the severity of ongoing symptoms of schizophrenia, thus highlighting direct clinical relevance to the measurement. This is just one example of how beta band oscillations have been identified as a potential biomarker of disease; other examples include Parkinson's Disease (Timmermann and Florin 2012). In addition, the robustness of MRBD and PMBR has meant that they have also been used in neuroscience applications ranging from brain computer interfaces (Pfurtscheller and Solis-Escalante 2009) to markers of neural plasticity (Gaetz et al. 2010;Mary et al. 2015). It is also noteworthy that the beta band power loss and rebound, whilst commonly thought of as being observable in the sensorimotor cortex, is not a sole property of the sensorimotor system. For example, Fig. 3 shows instances of observation of very similar phenomena in other cortical areas. Figure 3a shows the time-frequency oscillatory dynamics of a network of brain areas encapsulating bilateral insula cortex, throughout a cognitive task (Liddle et al. 2016;Brookes et al. 2012). The task itself involves presentation of a series of visual stimuli; some stimuli are relevant to the task, others irrelevant. Subjects were asked to respond if the relevant stimuli match some predetermined condition. Note here that only the non-target stimuli are shown (meaning that the subjects did not actually make a response). In the relevant condition, clear beta modulation is observed with a decrease in amplitude followed by a rebound above baseline. Furthermore, this effect was also shown to be abnormal in schizophrenia, again demonstrating its clinical relevance. Figure 3b shows the case for simple sensory stimulation of the visual cortex . Here, subjects were asked to passively view a drifting visual grating; the figure shows the envelope of beta band oscillations throughout the task. Note again the clear structure with a loss in beta amplitude during stimulation and an increase on stimulus cessation. These represent two simple examples which show that the beta band effect is not simply a property of the sensorimotor system, but rather is a ubiquitous effect that is observed robustly across many cortical regions. The above indicates that stimulus related beta power loss and post stimulus rebound are generally observable effects, seen in many cortical areas, during both sensory and cognitive tasks. Further, the reduction of rebound in disease has been robustly demonstrated. Thus, the generation of new mathematical models from which we can accurately predict task induced beta band dynamics are of interest.

A mean field model for spiking networks
There are a now a plethora of single neuron models for describing the spiking dynamics of cortical cells, many of which are extensions of the basic Hodgkin-Huxley model to incorporate nonlinear ionic currents that allow low frequency firing in response to constant current injection. Importantly mathematical neuroscience has identified a number of mechanisms that can generate 'f-I' curves with this property, with perhaps the most well known being the saddle-node on an invariant circle (SNIC) bifurcation (Ermentrout and Kopell 1986). This has led to the formulation of the elegant θ -neuron model (or Ermentrout-Kopell canonical model) which can mimic the firing and response properties of a cortical cell with a purely one dimensional dynamical system evolving on a circle. In a certain limit this is formally equivalent to the quadratic integrate-and-fire (QIF) model, also designed explicitly for understanding the generation of low firing rates in cortex (Latham et al. 2000). Given the simplicity of these models they are a natural candidate for cortical network studies, not only because they are computationally cheap, but because there is more chance to develop a statistical neurodynamics for such models than their biophysically complicated conductance based counterparts. Indeed mean field models for globally pulse-coupled networks have recently been developed by Luke et al. (2013) for θ -neuron models, and by Montbrió et al. (2015) for QIF models. To make these models more relevant to the interpretation of brain imaging signals, and in particular EEG and MEG, it is vital to augment the networks with more realistic forms of synaptic interaction and to move away from the overly restrictive assumption of fast pulsatile synaptic currents.
We consider a network of N QIF neurons each with a voltage v i , for i = 1, . . . , N, evolving according according to the following set of coupled ordinary differential equations (ODEs): Here C is a capacitance, η i is a constant current drive, κ is a proportionality constant, which from now on (without loss of generality) will be set to unity, and I i is the synaptic current, where g(t) represents a common time-dependent synaptic drive, which we shall take to arise through global coupling. This acts to push the voltage toward the synaptic reversal potential v syn . If the synaptic current is positive (negative) we say that the synapse is excitatory (inhibitory). The QIF network has discontinuous trajectories since whenever v i reaches a threshold value v th it is reset to the value v reset . This firing condition is also used to define an implicit set of firing times according to v i (T m i ) = v th , where T m i denotes the mth firing time of the ith neuron. These in turn can be used to generate a set of conductance changes for the ith neuron that we write in the form m∈Z s(t − T m i ), where s(t) is a fixed temporal filter. For a globally coupled network, with strength of synapse k/N, the total synaptic conductance change at each neuron is then For fast pulsatile interactions we may set s(t) = δ(t), where δ is a Dirac-delta function. For a more realistic form describing a normalised post synaptic potential (PSP) with an exponential decay we may set s(t) = αe −αt (t), whilst for a more general PSP with both a rise and fall time we would set s( Here (t) is a Heaviside step function included to enforce causality, and the parameters α, α 1,2 are decay rates. Exploiting the fact that in all these cases s(t) is the Green's function of The first example shows pulsatile coupling, where no synaptic filtering has been applied to the incoming spike. The second type of filter shown is an exponentially decaying function, which accounts for the slow decay of the instantaneous pulse. It does not however account for the time it take a synapse to process the incoming action potential, increasing instantaneously to the maximum value as soon as the spike arrives. The last example takes into account this synaptic processing delay, increasing smoothly to its peak value and then decaying exponentially back to zero. Note that (t) represents the Heaviside function a linear differential operator Q then we may also write g as the solution to the ODE system For the corresponding operator Q to the choice of s see Table 1. For the rest of this paper we shall work with the choice s(t) = α 2 te −αt , describing the so-called α-function. This can be obtained from the difference of exponentials form described above in the limit α 1,2 → α, so that the corresponding differential operator Q is Thus Eqs.
(1)-(5) define our spiking model. This is further illustrated in Fig. 4 for an all-to-all coupled neural network, with an inset showing the single neuron dynamics. Each neuron generates a train of spikes described by a sequence of Dirac-delta functions that are then filtered with the kernel s(t) to generate a synaptic current according to Eq.
(2). From this one may in principle use Maxwell's equations to determine the magnetic field that would underly an MEG-like signal. However, for simplicity we shall take the average network current to be a proxy for this physiological signal. This is given explicitly by which describes the average membrane potential. For a single uncoupled neuron (k = 0) it is a simple matter to show that, when v th → ∞ and v reset → −∞, the frequency of oscillation is given by 2 √ η i /C. For further simplicity we shall work with these choices for the values of v th and v reset . From now on we shall choose η i to be random variables drawn from a Lorentzian distribution L(η) with where η 0 is the centre of the distribution and the width at half maximum. For simplicity we will consider the average frequency of oscillations ω 0 = 2 √ η 0 /C as a single parameter, distributed with a width at half maximum ω = 2 √ /C. In the coupled network, and if the frequencies of the individual neurons are similar enough, then one may expect some degree of phase locking (ranging from synchrony to asynchrony), itself controlled in part by the time-to-peak, 1/α, of the synaptic filter.
As shown in Fig. 5, for a model with predominantly inhibitory connections, we see patterns of coherent spiking. The degree of coherence is mainly controlled by the degree of heterogeneity of the constant current drives η i .  3) and (4). The top plot of the zoom shows the shape of the synaptic filter for the case that s(t) is an α-function: s(t) = α 2 te −αt , where 1/α is the time-to-peak. I i is the total synaptic current that enters the cell body and v i is the voltage of the cell which oscillates as shown in the middle plot. The corresponding output spike train is given by a sequence of Dirac-delta functions δ i = m∈Z δ(t − T m i ), as illustrated in the bottom plot In this figure we also track the evolution of two macroscopic order parameters. These are respectively the average membrane potential, given by Eq. (6), and the instantaneous mean firing rate r: For large N both the order parameters V and r show a smooth temporal variation. In the case of synchrony we would expect these mean field signals to show a periodic temporal variation, essentially following a trajectory reminiscent of a single QIF neuron receiving periodic drive, whilst for an asynchronously firing population these mean field signals would be constant in time (modulo finite size fluctuations). To quantify the degree of coherence (or phaselocking) within an oscillatory population it is convenient to use a Kuramoto order parameter. First though it is necessary to recast the model in terms of phase variables.
Given the well known link between the QIF neuron and the θ -neuron it is natural to introduce the phase variable θ i ∈ [−π, π) according to v i = tan(θ i /2) (so that cos θ i = (1 − v 2 i )/(1 + v 2 i ) and sin θ i = 2v i /(1 + v 2 i )). In this case we arrive at the θ -neuron network Here P (θ) = δ(θ − π) and is periodically extended such that P (θ) = P (θ + 2π), and we have used the result that The network defined by Eqs. (9) and (10) describes a set of N phase variables interacting via spike triggered currents every time that θ j passes through π. We will only consider the case that θ j increases through π (so that spikes are only generated on the upswing of the corresponding voltage variable). In the absence of synaptic coupling an isolated θ -neuron supports a pair of equilibria θ ± , with θ + < 0 and θ − > 0 for η i < 0, and no equilibria for η i > 0. In the former case the equilibria at θ + is stable and the one at θ − unstable. In neurophysiological terms, the unstable fixed point at θ − is a threshold for the neuron model. Any initial conditions with θ ∈ (θ + , θ − ) will be attracted to the stable equilibrium, while initial data with θ > θ − will make a large excursion around the circle before returning to the rest state. For η i > 0 the θ -neuron oscillates with frequency 2 √ η i /C. When η i = 0 the θ -neuron is poised at a saddle-node on an invariant circle (SNIC) bifurcation.
As well as naturally providing a phase variable the θneuron network is more straight forward to simulate as the model has continuous trajectories on an N-torus (and there is no need to handle the discontinuous reset conditions). The Kuramoto order parameter is then defined as Here R provides a measure of the degree of coherence within the network and is the average phase. If a population is perfectly synchronised then R = 1 and similarly if the system is perfectly asyncronous then R = 0. In Fig. 6 we show a sequence of snapshots of the Kuramoto order parameter for the dynamics shown in Fig. 5, as well as the time evolution of the degree of coherence. Much as the order parameters (r, V ) vary smoothly with time (for large N) then so does the pair (R, ). Indeed in the large N limit Luke et al. (2013), making use of the Ott-Antonsen (OA) ansatz (Ott and Antonsen 2008), have shown that a globally pulse-coupled network of θ -neurons has an exact mean field description. The OA ansatz allows the reduction of globally coupled phase oscillators in the infinite size limit to an explicit finite set of nonlinear ODEs. These describe the macroscopic evolution of the system, in terms of the Kuramoto order parameter for synchronisation, as long as the distribution of phases is at most single peaked. In essence the ansatz is well suited to describing systems which dynamically evolve between an incoherent asynchronous state and a partially synchronised state, which is often the case in systems with interactions that are prescribed by harmonic functions, such as found in Eq. (9). To illustrate the type of network evolution that can be generated with different values of synchrony, see Fig. 7. Here we show some plots of the phase distribution for different values of the network coherence as well as the average network current that would be produced. Interestingly, Montbrió et al. have recently shown how to move between order parameters for the phase and voltage descriptions with the use of a conformal transformation (Montbrió et al. 2015). If we introduce the complex order parameter for the voltage description as then the following transformation allows one to switch perspectives and obtain the order parameter for the voltage description in terms of the Kuramoto order parameter as where Z denotes the complex conjugate of Z. Importantly the OA ansatz can also be used to obtain a mean field model in the presence of non-pulsatile synaptic, extending the approaches in Luke et al. (2013) and Montbrió et al. (2015). In Appendix B we show that this yields the mean field model described by the fourth order ODE system: Qg = kH (Z). and the average synaptic current is very spiky. In the regime where R 0.5 the phases are more distributed. Although a dominant phase can be clearly identified (by the peak value) not all neurons have this phase. The OA ansatz gives the shape of the distribution in the form F (θ) = (2π) −1 (1−|Z| 2 )/(|e iθ −Z| 2 ). This spread in the phase distribution acts to smooth out the spikes in the average synaptic current to create a smooth oscillatory signal. When the population of neurons is completely asynchronous there is no dominant phase and every phase is equally probable so that F (θ) = 1/(2π). In this case all of the neurons fire at different times with their phases uniformly distributed to yield a constant synaptic current. Note that the peak in the distribution of phases move as the system evolves with time with a velocity˙ Here H (Z) is a global state dependent drive to the population given by A plot of this scalar function of a complex variable is shown in Fig. 8. It is illuminating to express H as function of W using (13) from which we find Hence we may interpret H (Z) as the firing rate of the population that drives the global synaptic current. Figure 8 shows H as a function of Z. As expected H takes its highest value when Z e iπ , corresponding to high synchrony where all of the neurons fire and reset at the same time. Figure 9 shows results for a simulation of 500 theta neurons (red) and a simulation of the reduced mean field model (blue). It is strikingly clear that the two simulations agree very well. If the size of the population in the large scale simulations is reduced then one can begin to see finite size fluctuations, as expected. The macroscopic order parameters (r, V ) in the reduced mean field model are plotted in Fig. 10. Unsurprisingly they behave similarly to the corresponding order parameters for the large scale simulations plotted in Fig. 5. Likewise, the mean field representation of (R, ), plotted in Fig. 11, agree extremely well with those  Fig. 6. For a further discussion of the bifurcation structure of this model see Coombes and Byrne (2017).

A mechanistic interpretation of movement induced changes in the beta rhythm
In Section 2 we demonstrated how an externally cued thumb movement caused a ∼ 0.5 s decrease in beta band power followed by a ∼ 2 − 4 s increase in beta band power, typifying MRBD and PMBR, respectively. The median nerve stimulation lasts ∼ 50 ms, however the evoked response lasts significantly longer. Upon examining the time-frequency spectrograms in Fig. 1 we observed an increase in low frequency activity at t = 0, which appears to last for ∼ 0.3 − 0.4 s, corresponding to the transduced median nerve stimulation and corresponding movement. We base the design of the external drive on this transduced signal.
We model the transduced signal as a temporally filtered drive A = A(t) that is received by every neuron in the Fig. 10 Mean field reduction for a QIF network: Time series for the mean field variable W = πCr + iV , where r is the population firing rate and V is the average voltage. Comparing these plots to the corresponding plots for a 500 neuron simulation in Fig. 5 it is clear to see that they agree well. The finite size fluctuation for V are quite apparent when comparing the results for the large scale simulation to those of the reduced mean field model. However, the overall behaviour is similar. Parameters as in Fig. 5 model. In this case the dynamics of Z obey (14) under the replacement η 0 → η 0 + A, with Q D A(t) = (t), where Q D is the differential operator obtained from Q in Eq. 5 under the replacement α → α D , and (t) is a rectangular pulse, (t) = (t) (τ − t), where can be interpreted as the strength of the drive. Note that the pulse is not applied until after transients have dropped off. As the evoked response in the experimental data last for ∼ 0.3−0.4 s we set τ = 0.4 s. Fig. 9 Validity of reduction: Comparison between the reduced mean field network (blue) and simulation a network of 500 θ -neurons (red). Phase plane for the Kuramoto order parameter Z = Re i is shown on the left and the phase plane for the synaptic conductance g and its derivative g is shown on the right. Parameter values as in   Figure 12 shows the phase plane for the Kuramoto order parameter Z = Re i , as well as a time series for the within population synchrony R, in response to the drive described above. The colours correspond to the different time periods; before drive (blue), during drive (red), after drive (green). The system oscillates in partial synchrony with R oscillating between ∼ 0.05 − 0.6 in the absence of drive. Once the drive is switched on the amplitude of these oscillations decreases and hence the power is also reduced, corresponding to MRBD. Note that the frequency also increases during this period. After the drive is switched off the level of coherence is increased as Z is drawn towards the edge of the unit disk before spiralling back to the original limit cycle, corresponding to PMBR. Importantly the system does not rebound until t 0.5 s as seen in the real data. It should be noted that the stimulus corresponds to ∼ 80% of this time to rebound, however as the evoked response is present in ∼ 60 − 80% of the ∼ 0.5 s of MRBD in the real data we believe that this is a good fit.
The corresponding response of the synaptic current is shown in Fig. 13. The time series (left) shows that when the drive is switched on the synaptic current is reduced, however the neurons are now also receiving a strong excitatory current in the form of the drive. There is a large increase in the amplitude of the oscillations of the synaptic current at Fig. 12 Response of Z to drive Top: Phase plane for Z, demonstrating the behaviour of the system in response to the drive A(t). The blue curve represents the system before the pulse arrives, as it settles to its non-perturbed dynamics (t < 0), the red curve demonstrates how the system behaves when the pulse is switched on (0 < t < τ) and the green how the system reacts once the drive is switched off (t > τ s). Bottom: Time series of R showing the change in the level of coherence, before, during and after the drive is switched on. The amplitude of the oscillations in R appear significantly reduced while the drive is switched on. Parameter values as in Fig. 5, apart from ω 0 which was increased to 0.279 Hz, as this value gave a stronger PMBR, τ = 0.4, 1/α D = 5.6 ms and = 15 mA ∼ 0.5 s, corresponding to PMBR, which can also be seen in the time-frequency spectrogram (right). The initial increase in amplitude is very large, however the percentage increase between 0.5 − 1.5 s is relatively small. The synaptic current appears to have fully settled back to its pre-drive behaviour by t 1.5 s, indicating a PMBR of roughly 1 s, which is not as long as the PMBR seen in our experimental data. An increase in power can be seen at around 26 Hz at t 0 s, corresponding to the increase in frequency during the drive on period. This high beta activity can be interpreted as the processing of the motor input.
Interestingly, we see a direct correlation between synchrony and the synaptic current. The time series in Fig. 12 (bottom) shows a peak in synchrony at t 0.5 s, just as the time series in Fig. 13 (left) shows a sharp increase in the amplitude of the synaptic current at t 0.5 s. This increase in amplitude can also be seen in the spectrogram (right). The strength of the drive dictates the extent of drive (red), after drive (green). Both figures clearly demonstrate the rebound of the system, there is an increase in amplitude (and hence power) ∼ 0.5 s after the drive was initially switched on. Parameters as in Fig. 12 the rebound. However it also prescribes the frequency of the oscillations amid the period when the drive is switched on. Therefore it is important to find the balance, where we have a prominent PMBR but also a physically realistic frequency during the interval when the stimulus is switched on. Although the increase in power at a higher frequency cannot be seen in our experimental data (Fig. 1), these timefrequency spectrograms were calculated for a small area of the motor cortex, it can be seen in the results obtained in Robson et al. (2015) (Fig. 2). It is widely believed that an increase in high beta and gamma activity is present in motor preparation and execution, in a more frontal region of the motor cortex. The parameters were chosen such that the system oscillated at beta frequency and a significant MRBD and PMBR could be observed. The model is robust and can reproduce MRBD and PMBR for a wide region of parameter space.

Discussion
We have presented a mechanistic model that exhibits both MRBD and PMBR. This low dimensional model is derived from a corresponding high dimensional spiking network model and maintains a faithful representation of synaptic currents. In the reduced model these currents are driven by a firing rate that is itself a function of the complex Kuramoto order parameter. This makes a significant departure from the usual phenomenological neural mass description of neuronal population dynamics for which the firing rate is usually only a function of synaptic activity or mean membrane potential. Importantly the transient response of the reduced model is sufficiently rich to capture the emergent time scales of both MRBD and PMBR, when it is stimulated whilst operating in the beta frequency range. Although the length of PMBR observed in the reduced model was shorter than that seen in the experimental data, it is still within the documented 1 − 10 s range. Given that model responses are linked to changes in within-population coherence, this gives further support to the notion that beta band amplitude changes, and in particular those in MRBD and PMBR, are in fact due to changes in synchrony. More generally the model parameters can be altered so that the population oscillates at other frequencies, and hence, used to explain other ERD/ERS phenomena in the brain.
One natural extension of the model is to small networks, with coupling through the mean-field variables of a set of local populations, to describe systems with a mixture of excitation and inhibition. This would lead to a richer set of structures within the phase space for the network and provide further mechanisms for controlling emergent time-scales (say as orbits can be made to approach saddle structures, leading to a slow down in dynamics), which may help lengthen the PMBR, see Coombes and Byrne (2017) for a discussion of the bifurcation structure for a two population model. An interesting study would be to couple two identical populations, corresponding to left and right motor cortex areas, and drive one of the populations to inspect the bilateral response as seen in experimental data (Robson et al. 2015;Liddle et al. 2016). It is also possible that noise may play a constructive role, and the OA ansatz for a meanfield reduction can also be performed in this case (Lai and Porter 2905). The introduction of noise may also result in a lengthening of the PMBR.
Perhaps of more interest in using these next generation neural models to replace existing neural mass descriptions, is the development of reductive techniques to handle other nonlinear integrate-and-fire models, such as those of Izhikevich type (Izhikevich 2003). However, this is a substantial mathematical challenge since the OA reduction technique that we have employed here will break down. However, it is possible to extend the work presented here to include gap junction coupling (Laing 2015). There is now little doubt that gap junctions play a substantial role in the generation of neural rhythms, both functional (Hormuzdi et al. 2004;Bennet and Zukin 2004) and pathological (Velazquez and Carlen 2000;Dudek 2002). It is also interesting to consider the spatially extended version of this model, we report on this elsewhere .
Another possible extension would be to include a variety of different synaptic receptor. We have assumed that PMBR and MRBD are mediated by the same type of synaptic receptor. However, Hall et al. (2011) suggest that MRBD is a GABA-A mediated process, whilst PMBR appears to be generated by a non-GABA-A receptor mediated process. A further model that distinguishes between receptors, may offer important insights into motor processes, and can be readily accomplished within the framework that we have presented here.

Compliance with Ethical Standards
Conflict of interests The authors declare that they have no conflict of interest.
Ethical approval All procedures performed in studies involving human participants were in accordance with the ethical standards of the University of Nottingham Medical School Ethics Committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

A.1 Data collection
In order to demonstrate directly the robustness of beta band modulation in sensorimotor cortex, a somatosensory paradigm was used. Two subjects took part in the study, which was approved by the University of Nottingham Medical School Ethics Committee. The paradigm comprised electrical stimulation of the subject's left median nerve, which was achieved by applying a series of 500 μs duration current pulses to two gold electrodes placed on the subject's wrist. The current was delivered using a Digitimer DS7A constant current stimulator, and the amplitude was increased slowly until a visible movement of the thumb was observed. Each experimental run comprised a total of 80 pulses delivered with an inter-stimulus-interval (ISI) of 10 s. A single experimental run lasted approximately 13 minutes. Each subject performed the study twice to assess robustness.
MEG data were captured using the third order synthetic gradiometer configuration of a 275-channel CTF wholehead MEG system (MISL, Port Coquitlam, Canada). The subject was positioned supine, with their head in the MEG helmet, whilst data were recorded at a 600 Hz sampling rate. Three localisation coils were attached to the head as fiducial markers (nasion, left preauricular and right preauricular) prior to the recording. Energising these coils at the start and end of data acquisition enabled localisation of the fiducial markers relative to the MEG sensor geometry as well as determination of total head movement. In order to coregister brain anatomy to the MEG sensor array, prior to the MEG recording each subject's head shape was digitised relative to the fiducial markers using a 3D digitiser (Polhemus IsoTrack, Colchester, VT, USA). Volumetric anatomical MR images were acquired using a 3 T MR system (Phillips Achieva, Best, Netherlands) running an MPRAGE sequence (1 mm 3 resolution). Following data acquisition, the head surface was extracted from the anatomical MR image and coregistered (via surface matching) to the digitised head shape for each subject. This allowed complete coregistration of the MEG sensor array geometry to the bain anatomy, thus facilitating subsequent forward and inverse calculations.

A.2 Data analysis
MEG data were inspected visually and any trials containing excessive interference removed. Data were then analysed using synthetic aperture magnetometry (SAM) (Vrba and Robinson 2001), a beamforming variant (van Drongelen et al. 1996;Gross et al. 2001;Hillebrand et al. 2005;Robinson and Vrba 1998;van Veen et al. 1997) that has been applied successfully in many studies to localise neural oscillatory amplitude changes. Data were first filtered to the beta (13-30 Hz) band. Following this the SAM beamformer was applied and oscillatory amplitude was contrasted in active and control time windows in order to delineate the spatial signatures of beta amplitude change. The forward model was based upon a multiple local sphere head model and the forward calculation by Sarvas (Huang et al. 1999;Sarvas 1987). Covariance was computed within the prescribed windows (Brookes et al. 2008). Pseudo-tstatistical images were generated.

Appendix B: Mean field reduction
In the limit N → ∞ the state of the network at time t can be described by a continuous probability distribution function ρ (η, θ, t), which satisfies the continuity equation (arising from the conservation of oscillators): ∂ρ ∂t where c is a given realisation ofθ , The global drive to the network, given by the right hand side of Eq. (10), can be constructed as where we have used the result that 2πP (θ) = m∈Z e im(θ−π) . The formula for c above may be written conveniently in terms of e ±iθ as where f = ((η−1)+v syn g +ig)/2 and h = (η+1)+v syn g, and f denotes the complex conjugate of f . The OA ansatz assumes that ρ(η, θ, t) has the product structure ρ(η, θ, t) = L(η)F (η, θ, t) where L(η) is taken to be the Lorentzian given by Eq. (7). Since F (η, θ, t) should be 2π periodic in θ it can be written as a Fourier series: where cc denotes complex conjugate. The insight in Ott and Antonsen (2008) was to restrict the Fourier coefficients such that F n (η, t) = α(η, t) n , where |α(η, t)| ≤ 1 to avoid divergence of the series. There is also a further requirement that α(η, t) can be analytically continued from real η into the complex η-plane, and that this continuation has no singularities in the lower half η-plane, and that |α(η, t)| → 0 as Im η → −∞. If we now substitute (22) into the continuity Eq. (18), use the OA ansatz, and balance terms in e iθ we obtain an evolution equation for α(η, t) as It is now convenient to remember the Kuramoto order parameter (11), which in the large N limit is given by We note that |Z| ≤ 1. Using the OA ansatz (and using the orthogonality properties of e iθ , namely 2π 0 e ipθ e iqθ dθ = 2πδ p+q,0 ) we then find that By noting that the Lorentzian (7) has simple poles at η ± = η 0 ±i the integral in Eq. (27) may be performed by choosing a large semi-circle contour in the lower half η-plane. This yields Z(t) = α(η − , t), giving Z(t) = α(η + , t). Hence, the dynamics for g given by Eq. (21) can be written as where with |Z| < 1. The dynamics for Z is obtained from Eq. (25) as