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

Journal of Computational Neuroscience, Jul 2017

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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:


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

A mean field model for movement induced changes in the beta rhythm A´ine Byrne 0 1 2 Matthew J Brookes 0 1 2 Stephen Coombes 0 1 2 Action Editor: Maxim Bazhenov 0 1 2 0 Sir Peter Mansfield Imaging Centre, School of Physics and Astronomy, University of Nottingham, University Park , Nottingham, NG7 2RD , UK 1 Centre for Mathematical Medicine and Biology, School of Mathematical Sciences, University of Nottingham, University Park , Nottingham, NG7 2RD , UK 2 A ́ ine Byrne 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 stimulusrelated 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. Post-movement beta rebound; Movement related beta decrease; Neural mass; Synchrony; Power spectra; Magnetoencephalography; MEG; Electroencephalography; EEG; Mean field 1 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 (30 − 200 Hz). Moreover, it has long been known, since the early works of Jasper and 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 (Stanca´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 (Brookes et al. 2011; 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 106−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 (Pazo´ and Montbrio´ 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. 2 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 Fig. 1 Robust beta rebound for median nerve stimulation: Timefrequency spectrograms showing the percentage change from baseline of the activity in the motor cortex, for 2 participants on two separate days. The top row shows the results for participant 1 and the bottom row represents data from participant 2. Each participant displayed a clear difference in their PMBR. However, the strength and timing of both the MRBD and PMBR are consistent between subjects and trials. PMBR for each participant appears to be reproducible 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. (2011, 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 Fig. 2 The beta rebound and its disruption in patients with schizophrenia: (a) Time-frequency spectrograms showing changes in the amplitude of neural oscillations, in contralateral sensorimotor cortex, when subjects execute a 2s finger movement. Note that, in the beta band the loss in oscillatory power during movement is accompanied by 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 an increase in power on movement cessation. (b) Equivalent timefrequency spectrogram in patients with schizophrenia. Note the significant reduction in PMBR (Figure reproduced with permission from Robson et al. 2015.) 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 Fig. 3 Task induced beta band decrease and rebound phenomena in other cortical regions: (a) Time-frequency dynamics in a network of brain areas including bilateral insula. The task involved visual stimuli that were relevant and irrelevant to the task. Note the significant reduction and rebound in beta oscillations in the relevant condition. (Reproduced with permission from Liddle et al. 2016.) (b) Timecourse showing the envelope of beta oscillations in primary visual cortex during passive viewing of a visual grating. Visual stimulation occurred in the 0 s to 4 s window. Note again the task induced power loss and post stimulus rebound (Reproduced with the authors’ permission from Stevenson et al. 2011.) 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 (Stevenson et al. 2011). 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. 3 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 (2) (3) 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 Montbrio´ 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 vi , for i = 1, . . . , N , evolving according according to the following set of coupled ordinary differential equations (ODEs): d C dt vi = ηi + κvi2 + Ii , 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 Ii is the synaptic current, i = 1, . . . , N . (1) Ii = g(t )(vsyn − vi ), 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 vsyn. If the synaptic current is positive (negative) we say that the synapse is excitatory (inhibitory). The QIF network has discontinuous trajectories since whenever vi reaches a threshold value vth it is reset to the value vreset. This firing condition is also used to define an implicit set of firing times according to vi (Tim) = vth, where Tim 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 − Tim), 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 k N g(t ) = N 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(t ) = (1/α1 − 1/α2)−1[α1e−α1t − α2e−α2t ] (t ). 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 k Qg(t ) = N N j =1 m∈Z 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 ) = α2t e−α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 Q = 1 d 1 + α dt 2 . Thus Eqs. (1)–(5) define our spiking model. This is further illustrated in Fig. 4 for an all-to-all coupled neural (4) (5) Fig. 4 Neural network: The diagram on the right shows an all-to-all coupled network. The zoom on the left shows each of the components of Eqs. (1), (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 ) = α2t e−αt , where 1/α is the time-to-peak. Ii is the total synaptic current that enters the cell body and vi 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 as illustrated in thme∈bZotδt(otm−pTloimt), functions δi = 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 g(t )(vsyn − V (t )), where 1 V (t ) = N N which describes the average membrane potential. For a single uncoupled neuron (k = 0) it is a simple matter to show that, when vth → ∞ and vreset → −∞, the frequency of oscillation is given by 2√ηi /C. For further simplicity we shall work with these choices for the values of vth and vreset. From now on we shall choose ηi to be random variables drawn from a Lorentzian distribution L(η) with 1 L(η) = π (η − η0)2 + 2 , 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 2pa√ram/eCte.r,Indistthreibucoteudpwleidthnaetwwiodrtkh,aatnhdalfifmtahxeimfruemquenωcie=s 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 . (6) (7) 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: Given the well known link between the QIF neuron and the θ -neuron it is natural to introduce the phase variable θi ∈ [−π, π ) according to vi = tan(θi /2) (so that cos θi = (1 − vi2)/(1 + vi2) and sin θi = 2vi /(1 + vi2)). In this case we arrive at the θ -neuron network d C dt θi = (1−cos θi )+(1+cos θi )(ηi +gvsyn)−g sin θi , (9) (10) N P (θj ). j=1 Here P (θ ) = δ(θ − π ) and is periodically extended such that P (θ ) = P (θ + 2π ), and we have used the result that δ(t − T m j ) = δ(θj (t ) − π )|θ˙j (Tjm)|. 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 1 N Z(t ) = N j=1 eiθj (t) ≡ R(t )ei (t). (11) 1 N r(t ) = N 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. (8) 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 middle plot demonstrates that the length from the centre of the disk to the black dot represents R, the population synchrony, and , the average phase, is represented by the angular position of the black dot. The top right plot shows the system at a later point in time. The lower plots show R and as a function of time. One observes that both the population synchrony and average phase continuously oscillate. Parameter values as in Fig. 5 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, Montbrio´ et al. have recently shown how to move between order parameters for the phase and voltage descriptions with the use of a conformal transformation (Montbrio´ et al. 2015) . If we introduce the complex order parameter for the voltage description as W = π Cr + iV , W = 11 −+ ZZ , 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 Montbrio´ et al. (2015) . In Appendix B we show that this yields the mean field model described by the fourth order ODE system: 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 (13) (14) Fig. 7 Distribution of phases: Figure illustrating the distribution of phases F (θ ) in the large N limit and the average synaptic current for different values of the population coherence R. For simplicity we have fixed the choice of time so that (t ) = 0. When the population is completely synchronous (R = 1) all of the neurons have the same phase and as a result all of the neurons fire together so that F (θ ) = δ(θ ) 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 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 H (W ) = 1 Cπ W + W 2 = r. Hence we may interpret H (Z) as the firing rate of the population that drives the global synaptic current. Figure 8 shows (16) (17) phase. The OA ansatz gives the shape of the distribution in the form F (θ ) = (2π )−1(1−|Z|2)/(|eiθ −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 ˙ H as a function of Z. As expected H takes its highest value when Z eiπ , 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 shown in Fig. 6. For a further discussion of the bifurcation structure of this model see Coombes and Byrne (2017) . 4 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. 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 = Rei 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 Fig. 5 0.5 1 0 -0.5 -1 -0.5 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 QD A(t ) = (t ), where QD 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. Figure 12 shows the phase plane for the Kuramoto order parameter Z = Rei , 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 Fig. 13 Response of the synaptic current to drive: Time series and spectrogram of the synaptic current showing the response of the system to the external drive (t). The colours in the time series (left) correspond to the different time periods; before drive (blue), during 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. 5 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 (Byrne et al. 2017) . 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. Acknowledgements SC was supported by the European Commission through the FP7 Marie Curie Initial Training Network 289146, NETT: Neural Engineering Transformative Technologies. MJB was funded by a Medical Research Council (MRC) New Investigator Research Grant (MR/M006301/1). We also acknowledge Medical Research Council Partnership Grant (MR/K005464/1). Compliance with Ethical Standards Conflict of interests of interest. The authors declare that they have no conflict 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. Appendix A: Experimental method 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 mm3 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): 1 C −g sin θ ] . ∂ρ ∂t + ∂θ ∂ρc = 0, where c is a given realisation of θ˙, c = (1 − cos θ ) + (1 + cos θ )(η + gvsyn) The global drive to the network, given by the right hand side of Eq. (10), can be constructed as (18) 2π δp+q,0) we then find that Z(t ) = ∞ −∞ dηL(η)α(η, t ). It is now convenient to remember the Kuramoto order parameter (11), which in the large N limit is given by Z(t ) = 0 2π dθ ∞ −∞ dηρ(η, θ , t )eiθ . orthogonality properties of eiθ , namely We note that |Z| ≤ 1. Using the OA ansatz (and using the 2π eipθ eiqθ dθ = 0 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 dρ(η, θ , t )P (θ ). (20) dρ(η, θ , t )eim(θ−π), (21) lim N→∞ N Hence, Qg = where f = ((η −1)+vsyng +ig)/2 and h = (η +1)+vsyng, 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: 1 F (η, θ , t ) = 2π 1 + Fn(η, t )einθ + cc , ∞ n=1 where cc denotes complex conjugate. The insight in Ott and Antonsen (2008) was to restrict the Fourier coefficients such that Fn(η, 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 eiθ we obtain an evolution equation for α(η, t ) as ∂ 1 ∂t α + C iα2f + iαh + if (26) (27) (28) (29) (30) (31) with |Z| < 1. The dynamics for Z is obtained from Eq. (25) as − . (−1)mZm + cc C d dt Z = −i (Z − 1)2 2 + (Z + 1)2 2 − Z2 − 1 2 (19) (22) (23) (24) (25) Alegre , M. , Alvarez-Gerriko , I. , Valencia , M. , Iriarte , J. , & Artieda , J. ( 2008 ). Oscillatory changes related to the forced termination of a movement . Clinical Neurophysiology: Official Journal of the International Federation of Clinical Neurophysiology , 119 ( 2 ), 290 - 300 . Ashwin , P. , Coombes , S. , & Nicks , R. ( 2016 ). Mathematical frameworks for oscillatory network dynamics in neuroscience . Journal of Mathematical Neuroscience , 6 ( 2 ), 1 - 92 . Bennet , M.V.L. , & Zukin , R.S. ( 2004 ). Electrical coupling and neuronal synchronization in the mammalian brain . Neuron , 41 , 495 - 511 . Berger , H. ( 1929 ). U¨ ber das Elektroenzenkephalogram des Menschen . Archiv fu¨r Psychiatrie und Nervenkrankheiten , 87 , 527 - 70 . Brookes , M.J. , Vrba , J. , Robinson , S.E. , Stevenson , C.M. , Peters , A.M. , Barnes , G.R. , Hillebrand , A. , & Morris , P.G. ( 2008 ). Optimising experimental design for MEG beamformer imaging . NeuroImage , 39 ( 4 ), 1788 - 802 . Brookes , M.J. , Wood , J.R. , Stevenson , C.M. , Zumer , J.M. , White , T.P. , Liddle , P.F. , & Morris , P.G. ( 2011 ). Changes in brain network activity during working memory tasks: a magnetoencephalography study . NeuroImage , 55 ( 4 ), 1804 - 15 . Brookes , M.J. , Liddle , E.B. , Hale , J.R. , Woolrich , M.W. , Luckhoo , H. , Liddle , P.F. , & Morris , P.G. ( 2012 ). Task induced modulation of neural oscillations in electrophysiological brain networks . NeuroImage , 63 ( 4 ), 1918 - 30 . Byrne , A. , Avitabile , D. , & Coombes , S. ( 2017 ). A next generation neural field model: the evolution of synchrony within patterns and waves . Physical Review E. In prep. Donner , T.H. , & Siegel , M. ( 2011 ). A framework for local cortical oscillation patterns . Trends in Cognitive Sciences , 15 ( 5 ), 191 - 199 . Cassim , F. , Monaca , C.A.C. , Szurhaj , W. , Bourriez , J.-L. , Defebvre , L. , Derambure , P. , & Guieu , J.-D. ( 2001 ). Does post-movement beta synchronization reflect an idling motor cortex? Neuroreport , 12 ( 17 ), 3859 - 3863 . Cheyne , D.O. ( 2013 ). MEG studies of sensorimotor rhythms: a review . Experimental Neurology , 245 , 27 - 39 . Coombes , S. ( 2010 ). Large-scale neural dynamics: simple and complex . NeuroImage , 52 , 731 - 739 . Coombes , S. , & Byrne , A. ( 2017 ). Lecture notes nonlinear dynamics in computational neuroscience: from physics and biology to ICT, chapter next generation neural mass models . PoliTO Springer. Coombes , S. , beim Graben, P. , Potthast , R. , & Wright , J . (Eds.) ( 2014 ). Neural fields, theory and applications . Berlin: Springer. Dudek , F.E. ( 2002 ). Gap junctions and fast oscillations: a role in seizures and epileptogenesis? Epilepsy Currents , 2 , 133 - 136 . Ermentrout , G.B. , & Kopell , N. ( 1986 ). Parabolic bursting in an excitable system coupled with a slow oscillation . SIAM Journal on Applied Mathematics , 46 ( 2 ), 233 - 253 . Gaetz , W. , Macdonald , M. , Cheyne , D. , & Snead , O.C. ( 2010 ). Neuromagnetic imaging of movement-related cortical oscillations in children and adults: age predicts post-movement beta rebound . NeuroImage , 51 ( 2 ), 792 - 807 . Gaetz , W. , Edgar , J.C. , Wang , D.J. , & Roberts , T.P.L. ( 2011 ). Relating MEG measured motor cortical oscillations to resting γ - aminobutyric acid (GABA) concentration . NeuroImage , 55 ( 2 ), 616 - 21 . Gross , J. , Kujala , J. , Hamalainen , M. , Timmermann , L. , Schnitzler , A. , & Salmelin , R. ( 2001 ). Dynamic imaging of coherent sources Studying neural interactions in the human brain . Proceedings of the National Academy of Sciences of the United States of America , 98 ( 2 ), 694 - 9 . Hall , S.D. , Stanford , I.M. , Yamawaki , N. , McAllister , C.J. , Ro¨nnqvist, K.C., Woodhall , G.L. , & Furlong , P.L. ( 2011 ). The role of GABAergic modulation in motor function related neuronal network activity . NeuroImage , 56 , 1506 - 1510 . Hall , S.D. , Prokic , E.J. , McAllister , C. , Ronnqvist , K.C. , Williams , A.C. , Yamawaki , N. , Witton , C. , Woodhall , G.L. , & Stanford , I.M. ( 2014 ). Gaba-mediated changes in inter-hemispheric beta frequency activity in early-stage parkinson's disease . Neuroscience , 281 , 68 - 76 . Hillebrand , A. , Singh , K.D. , Holliday , I.E. , Furlong , P.L. , & Barnes , G.R. ( 2005 ). A new approach to neuroimaging with magnetoencephalography . Human Brain Mapping , 25 ( 2 ), 199 - 211 . Hipp , J.F. , Hawellek , D.J. , Corbetta , M. , Siegel , M. , & Engel , A.K. ( 2012 ). Large-scale cortical correlation structure of spontaneous oscillatory activity . Nature Neuroscience , 15 ( 6 ), 884 - 90 . Hormuzdi , S.G. , Filippov , M.A. , Mitropoulou , G. , Monyer , H. , & Bruzzone , R. ( 2004 ). Electrical synapses: a dynamic signaling system that shapes the activity of neuronal networks . Biochimica et Biophysica Acta , 1662 , 113 - 137 . Huang , M.X. , Mosher , J.C. , & Leahy , R.M. ( 1999 ). A sensor-weighted overlapping-sphere head model and exhaustive head model comparison for MEG . Physics in Medicine and Biology , 44 ( 2 ), 423 - 40 . Izhikevich , E.M. ( 2003 ). Simple model of spiking neurons . IEEE Transactions on Neural Networks , 14 , 1569 - 1572 . Jasper , H.H. , & Andrews , H.L. ( 1936 ). Human brain rhythms. I. Recording techniques and preliminary results . Journal of General Physiology , 14 , 98 - 126 . Jasper , H.H. , & Andrews , H.L. ( 1938 ). Brain potentials and voluntary muscle activity in man . Journal of Neurophysiology , 1 ( 2 ), 87 - 100 . Jasper , H.H. , & Penfield , W. ( 1949 ). Electrocorticograms in man: Effect of voluntary movement upon the electrical activity of the precentral gyrus . Archiv fu¨r Psychiatrie und Nervenkrankheiten , 183 ( 1-2 ), 163 - 174 . Jensen , O. , Goel , P. , Kopell , N. , Pohja , M. , Hari , R. , & Ermentrout , B. ( 2005 ). On the human sensorimotor-cortex beta rhythm: sources and modeling . NeuroImage , 26 ( 2 ), 347 - 355 . Jurkiewicz , M.T. , Gaetz , W.C. , Bostan , A.C. , & Cheyne , D. ( 2006 ). Post-movement beta rebound is generated in motor cortex: evidence from neuromagnetic recordings . NeuroImage , 32 ( 3 ), 1281 - 9 . Kilavik , B.r.E. , Zaepffel , M. , Brovelli , A. , MacKay , W.A. , & Riehle , A. ( 2013 ). The ups and downs of β oscillations in sensorimotor cortex . Experimental Neurology , 245 , 15 - 26 . Klinshov , V.V. , Teramae , J.N. , Nekorkin , V.I. , & Fukai , T. ( 2014 ). Dense neuron clustering explains connectivity statistics in cortical microcircuits . PloS One , 9 , e94292 . Lai , Y.M. , & Porter , M.A. ( 2905 ). Noise-induced synchronization, desynchronization, and clustering in globally coupled nonidentical oscillators . Physical Review E , 88 ( 01 ), 2013 . Laing , C.R. ( 2015 ). Exact neural fields incorporating gap junctions . SIAM Journal on Applied Dynamical Systems , page to appear. Latham , P.E. , Richmond , B.J. , Nelson , P.G. , & Nirenberg , S. ( 2000 ). Intrinsic dynamics in neuronal networks . I. Theory Journal of Neurophysiology , 83 , 808 - 827 . Liddle , E.B. , Price , D. , Palaniyappan , L. , Brookes , M.J. , Robson , S.E. , Hall , E.L. , Morris , P.G. , & Liddle , P.F. ( 2016 ). Abnormal salience signaling in schizophrenia: the role of integrative beta oscillations . Human Brain Mapping , 37 ( 4 ), 1361 - 74 . Luke , T.B. , Barreto , E. , & So , P. ( 2013 ). Complete classification of the macroscopic behaviour of a heterogeneous network of theta neurons . Neural Computation , 25 , 3207 - 3234 . Mary , A. , Bourguignon , M. , Wens , V. , Op de Beeck, M. , Leproult , R. , De Tie`ge, X. , & P. Peigneux. ( 2015 ). Aging reduces experience-induced sensorimotor plasticity. A magnetoencephalographic study . NeuroImage , 104 , 59 - 68 . Montbrio´ , E. , Pazo´, D. , & Roxin , A. ( 2015 ). Macroscopic description for networks of spiking neurons . Physica Review X , 5 ( 02 ), 1028 . Muthukumaraswamy , S.D. , Myers , J.F.M. , Wilson, S.J. , Nutt , D.J. , Hamandi , K. , Lingford-Hughes , A. , & Singh , K.D. ( 2013 ). Elevating endogenous GABA levels with GAT-1 blockade modulates evoked but not induced responses in human visual cortex . Neuropsychopharmacology: Official Publication of the American College of Neuropsychopharmacology , 38 ( 6 ), 1105 - 12 . Ott , E. , & Antonsen , T.M. ( 2008 ). Low dimensional behavior of large systems of globally coupled oscillators . Chaos , 18 , 037113 . Pazo´ , D. , & Montbrio´, E. ( 1009 ). Low-dimensional dynamics of populations of pulse-coupled oscillators . Physical Review X , 4 ( 01 ), 2014 . Pfurtscheller , G. , & Lopes da Silva , F.H. ( 1999 ). Event-related EEG/MEG synchronization and desynchronization: basic principles . Clinical Neurophysiology , 110 ( 11 ), 1842 - 1857 . Pfurtscheller , G. , & Solis-Escalante , T. ( 2009 ). Could the beta rebound in the EEG be suitable to realize a “brain switch”? Clinical Neurophysiology , 120 ( 1 ), 24 - 29 . Pfurtscheller , G. , Stanca´k, A ., & Neuper , C. ( 1996 ). Event-related synchronization (ERS) in the alpha band-an electrophysiological correlate of cortical idling: a review . International Journal of Psychophysiology , 24 ( 1-2 ), 39 - 46 . Pfurtscheller , G. , Neuper , C. , Brunner , C. , & da Silva , F.L. ( 2005 ). Beta rebound after different types of motor imagery in man . Neuroscience Letters , 378 ( 3 ), 156 - 9 . Pinotsis , D. , Robinson , P. , beim Graben, P. , & Friston , K. ( 2014 ). Neural masses and fields: Modelling the dynamics of brain activity . Frontiers in Computational Neuroscience 8 ( 149 ). Riehle , A. , & Vaadia , E. ( 2004 ). Motor cortex in voluntary movements: a distributed system for distributed functions . CRC Press. ISBN 0203503589. Robinson , S.E. , & Vrba , J. ( 1998 ). Functional neuroimaging by synthetic aperture magnetometry (SAM) . In Yoshimoto, T., Kotani , M. , Kuriki , S. , Karibe , H. , & Nakasato , N. (Eds.) Recent advances in biomagnetism (pp. 302 - 305 ): Tohoku University Press. Robson , S.E. , Brookes , M.J. , Hall , E.L. , Palaniyappan , L. , Kumar , J. , Skelton , M. , Christodoulou , N.G. , Qureshi , A. , Jan , F. , Katshu , M.Z. , Liddle , E.B. , Liddle , P.F. , & Morris , P.G. ( 2015 ). Abnormal visuomotor processing in schizophrenia . NeuroImage: Clinical. Sarvas , J. ( 1987 ). Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem . Physics in Medicine and Biology , 32 , 11 - 22 . Schnitzler , A. , Salenius , S. , Salmelin , R. , Jousma¨ki, V. , & Hari , R. ( 1997 ). Involvement of primary motor cortex in motor imagery: a neuromagnetic study . NeuroImage , 6 ( 3 ), 201 - 8 . So , P. , Luke , T.B. , & Barretto , E. ( 2014 ). Networks of theta neurons with time-varying excitability Macroscopic chaos, multistability, and final-state uncertainty . Physica D, 267 , 16 - 26 . Solis-Escalante , T. , Mu¨ller- Putz , G.R. , Pfurtscheller , G. , & Neuper , C. ( 2012 ). Cue-induced beta rebound during withholding of overt and covert foot movement . Clinical Neurophysiology: Official Journal of the International Federation of Clinical Neurophysiology , 123 ( 6 ), 1182 - 90 . Stanca´k, A ., & Pfurtscheller , G. ( 1995 ). Desynchronization and recovery of beta rhythms during brisk and slow self-paced finger movements in man . Neuroscience Letters , 196 , 21 - 24 . Stevenson , C.M. , Brookes , M.J. , & Morris , P.G. ( 2011 ). β-band correlates of the fMRI BOLD response . Human Brain Mapping , 32 ( 2 ), 182 - 97 . Timmermann , L. , & Florin , E. ( 2012 ). Parkinson's disease and pathological oscillatory activity: is the beta band the bad guy? - New lessons learned from low-frequency deep brain stimulation . Experimental Neurology , 233 ( 1 ), 123 - 5 . van Drongelen , W. , Yuchtman , M. , Van Veen , B.D. , & van Huffelen , A.C. ( 1996 ). A spatial filtering technique to detect and localize multiple sources in the brain . Brain Topography , 9 ( 1 ), 39 - 49 . van Veen , B.D. , van Drongelen, W. , Yuchtman , M. , & Suzuki , A. ( 1997 ). Localization of brain electrical activity via linearly constrained minimum variance spatial filtering . IEEE Transactions on Bio-medical Engineering , 44 ( 9 ), 867 - 80 . Velazquez , J.L.P. , & Carlen , P.L. ( 2000 ). Gap junctions, synchrony and seizures . Trends in Neurosciences , 23 , 68 - 74 . Vrba , J. , & Robinson , S.E. ( 2001 ). Signal processing in magnetoencephalography . Methods (San Diego California), 25 ( 2 ), 249 - 71 .

This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs10827-017-0655-7.pdf

Áine Byrne, Matthew J Brookes, Stephen Coombes. A mean field model for movement induced changes in the beta rhythm, Journal of Computational Neuroscience, 2017, 1-16, DOI: 10.1007/s10827-017-0655-7