many simple rhythmic behaviors are controlled by specialized central pattern generator ( cpg ) networks in the spinal cord or brainstem with intrinsic oscillatory characteristics ( grillner , 2006 ; kiehn , 2006 ) .
however , goal - directed upper - limb movements under cortical control can also exhibit rhythmicity . when tracking a moving target , trajectories are composed of multiple submovements ( craik , 1947 ) at a frequency ( usually one to four per second ) that is largely independent of movement speed ( miall et  al .
, 1986 ; roitman et  al . , 2004 ; pasalar et  al . , 2005 ; selen et  al . , 2006 ) .
an oscillation at around 3  hz has also been reported in the kinematics of finger - tracking movements ( mcauley et  al . , 1999 ) .
the motor cortical electroencephalogram ( eeg ) is phase locked to submovements ( dipietro et  al . , 2011 ) , and
coherence spectra between the magnetoencephalogram ( meg ) and movement speed show peaks around 3  hz during visuomotor tracking ( jerbi et  al . , 2007 ) .
brain - machine interface ( bmi ) studies have found low - frequency bands to be particularly informative for decoding direction from local field potentials ( lfps ; rickert et  al .
2011 ) , megs ( waldert et  al . , 2008 ) , and eegs ( waldert et  al . , 2008 ; bradberry et  al . ,
it has been argued that submovements reflect intermittent corrections driven by visual feedback of errors ( craik , 1947 ; miall et  al .
, 1986 , 1993 ) and that their frequency should therefore be  determined by extrinsic factors such as feedback loop delays . in support of this
submovements are locked to eye movements ( mcauley et  al . , 1999 ) and often disappear in the absence of vision ( miall et  al . , 1993 ;
mcauley et  al . , 1999 ; but see doeringer and hogan , 1998 ) , whereas the introduction of artificial feedback delays alters their frequency ( miall et  al . , 1986 ; miall and jackson , 2006 ) .
nevertheless , submovements are not restricted to tracking tasks , and a natural rhythmicity is observed across diverse upper - limb behaviors ( kunesch et  al . ,
1989 ) including self - paced isometric drawing ( massey et  al . , 1992 ) and finger tapping ( schner and kelso , 1988 ) .
moreover , low - frequency cortical oscillations have long been associated with slow - wave sleep , when large k complex potentials signifying transitions from down to up states of the cortex ( colrain , 2005 ; cash et  al . , 2009 )
are accompanied by bursts of activity in the delta ( 14  hz)-frequency range ( amzica and steriade , 1997 ) .
at least two mechanisms contribute to these delta oscillations : intrinsic currents that cause bursting patterns in thalamic relay cells ( amzica et  al .
, 1992 ; destexhe and sejnowski , 2003 ) and a second , purely cortical circuit ( amzica and steriade , 1998 ; carracedo et  al . ,
therefore , it remains possible that oscillatory properties of cortical ( and perhaps thalamic ) circuits contribute to low - frequency rhythms in movement , functioning much like a cpg .
recently , this  intrinsic hypothesis  has been proposed to explain the complex , multiphasic profiles of motor cortical firing rates observed during reaching ( churchland et  al . , 2012 ; shenoy et  al . ,
the high - dimensional neural state was projected onto a plane revealing low - frequency cycles , even though movements in this case were not overtly rhythmic . it was proposed that this dynamical structure reflects  an engine of movement 
( churchland et  al . , 2012 ) , and could be reproduced by a recurrent neural network model trained to generate muscle patterns given static initial inputs and no sensory feedback ( shenoy et  al . , 2013 ) .
the intrinsic hypothesis suggests that low - frequency cortical dynamics may be preserved across different movements and  resemble spontaneous delta oscillations during sleep .
therefore , we compared motor cortex activity in monkeys during an isometric movement task , while retrieving food from a klver board , during natural sleep , and under ketamine sedation .
we found clear evidence for a common lfp correlation structure that could be explained by a single model of 3  hz oscillatory dynamics underlying all behavioral states .
our results thereby unify two previously unrelated phenomena : low - frequency structure in movement kinematics and low - frequency oscillations during slow - wave sleep , providing a new insight into how the dynamics of cortical networks influence complex upper - limb behaviors .
three monkeys controlled the 2d position of a cursor with isometric wrist torque to acquire targets in a center - out fashion ( figure  1a ) . similar to isometric trajectories in humans
( massey et  al . , 1992 ) , the movements made by the monkeys were often composed of multiple , regular submovements . the representative single trial shown in figure  1b shows five submovements between the go cue and successful acquisition of the target , appearing as peaks in the radial speed of the cursor .
the distribution of intersubmovement intervals ( figure  2a ) , as well as their autocorrelation structure ( figure  2b ) , revealed a tendency for submovements to occur rhythmically at a frequency of around 3  hz in all three animals .
movement intermittency was also evident in the electromyogram ( emg ) from wrist muscles involved in the task , with a peak in coherence between radial position and rectified emg at 3  hz ( figure  2c ) .
lfps from multiple electrodes in primary motor cortex ( m1 ) were low - pass filtered and mean referenced , revealing slow oscillations with a phase that varied across electrodes ( figure  1c ) .
principal component analysis ( pca ) reduced the lfp to two principal components ( lfp - pcs ; figure  1d ) capturing orthogonal projections of the dominant oscillatory mode .
there was a striking correlation between lfp - pcs and submovements , evident even in single trials ( figure  1d , bottom ) .
coherence analysis over the whole session confirmed strong correlation between lfp - pcs and cursor speed at the 3  hz frequency of submovements ( figure  2d ) .
the relationship between lfp oscillations and submovements was visualized by projecting the lfp trajectory over time onto the pc plane ( movie s1 available online ) . due to the 90 phase difference between components , this trajectory was cyclic with constant direction of rotation .
each submovement was associated with a single lfp cycle , and the peak cursor speed occurred at a similar phase within each cycle ( figure  2e ) . to examine how lfp - pc trajectories were related to submovement kinematics , we binned submovements into nine equal groups according to peak cursor speed .
average lfp - pc trajectories for each group ( figures 2f and 2i ; movie s2 ) had constant direction and frequency of rotation ( figures 2 g and 2j ) but an areal velocity ( area swept out per unit time about the origin ; see equation  2 in experimental procedures ) that increased with cursor speed ( figures 2h and 2k ) . when binned according to the direction of submovements , average lfp trajectories in  the pc plane appeared similar ( figure  2l ) . however , plotting the trajectories in 3d pc space revealed a subtle variation in the axis of rotation for different submovement directions ( figure  2 m ; movie s3 ) . in three dimensions ,
we hypothesized that its magnitude should encode information about submovement speed , whereas its orientation might be informative of submovement direction .
we tested this directly by decoding speed ( or direction ) from the magnitude ( or orientation ) of areal velocity vectors associated with individual submovements ( see experimental procedures ) .
decoding performance for 13 data sets across three monkeys is summarized in table s1 , using a coefficient of determination ( cod ) between zero ( chance decoding ) and one ( perfect decoding ) . in every case
, we obtained significant decoding of the speed from areal velocity magnitude , with mean ( sd ) cod  = 0.30   0.13 .
decoding submovement direction from areal velocity orientation was significant in 12/13 sessions , with mean cod  = 0.15   0.07 . as might be expected , direction decoding was better for faster submovements ( figure  2n ) , with a cod of 0.26   0.17 for the fastest submovements .
although statistically significant , this nevertheless corresponds to an average decoding error of approximately 75 , only slightly better than chance ( 90 average decoding error ) . in summary ,
the rhythmic structure of submovements is reflected in low - frequency m1 lfps and can be revealed using pca to reduce the dimensionality of the multichannel data . in the space defined by the first three pcs , each submovement
is associated with a cyclic lfp trajectory , and the peak of the submovement occurs at a consistent phase of the cycle .
the areal velocity of the trajectory is proportional to the speed of the submovement , whereas the axis of rotation provides statistically significant ( albeit modest ) information about the direction of movement . on separate days , we recorded from the same electrodes during ketamine sedation ( figure  3a ) .
m1 lfps exhibited typical signatures of slow - wave sleep including spindles and large k complex potentials , thought to reflect transitions between cortical down and up states .
consistent with this , most neurons had low firing rates prior to each k complex , and fired particularly strongly during its rising phase and peak .
each k complex was associated with large - amplitude delta activity in the low - pass - filtered , mean - subtracted lfp , comprising either a single cycle or an extended burst of two or more cycles of low - frequency oscillation .
figure  4 compares cortical activity aligned to the peak speed of  submovements ( during task performance ) and aligned to the peak of k complexes ( under ketamine sedation ) .
both events were associated with phasic bursts of neural activity ( figure  4b ) and phase - locked low - frequency lfp oscillations in m1 ( figure  4c ) , which were an order of magnitude larger in the sedated state .
note that the same color scheme is used throughout figures 1 and 4c to represent the lfp phase relative to submovements during task performance but each electrode shows a similar phase relative to the k complex under sedation . as a result , average lfp - pc trajectories followed similar rotational cycles aligned to both submovements and k complexes ( figure  4d ) .
we calculated the circular - circular correlation coefficient ( cc ) across electrodes of phase relative to submovements ( during task performance ) against phase relative to k complexes ( under sedation ) .
for the sample pair of movement and sedation sessions shown in figure  4c , these phases were highly correlated ( n  = 10 , cc  = 0.81 , p  = 0.025 ; figure  4e ) , and across all the data sets this correlation was significant ( p  < 0.05 ) in 11/13 pairs of sessions in three monkeys ( mean   sd ; cc  = 0.75   0.25 ; table s1 ) .
pooling all the sessions for each animal also yielded significant correlation ( monkey d : n  = 52 , cc  = 0.63 , p  = 2   10 ; monkey r : n  = 36 , cc  = 0.49 , p  = 0.002 ; monkey s : n  = 45 , cc  = 0.81 , p  = 4   10 ; figure  4f ) .
moreover , advancing several microwires from the most superficial depth down through the gray matter revealed that both submovement- and k complex - related potentials underwent polarity reversals at  the same depth , indicating a common cortical source ( figure  s1 ) .
in summary , the same patterns of cortical activity seen during isometric movements ( and related to submovement kinematics ) also arose endogenously in the absence of behavior , suggesting that intrinsic circuits rather than extrinsic sensorimotor feedback loops impose this dynamical structure on low - frequency cortical activity .
because isometric torque tracking and ketamine sedation are somewhat unnatural experimental conditions , we proceeded to examine whether the same low - frequency dynamics were also present in brain activity during more naturalistic conditions including sleep and unrestrained reach - to - grasp .
two of our subjects often fell asleep at the end of recording sessions , providing an opportunity to examine low - frequency activity during natural sleep .
as the eyes closed , large - amplitude slow waves were observed in the lfp ( figure  3b ) , of a comparable amplitude to sedation recordings and approximately an order of magnitude greater than that seen in the awake state . however , cortical activity appeared disorganized and lacked clear up-/downstate transitions or k complexes , which is consistent with stage 1 sleep .
in addition , we collected data when the same animals retrieved food from wells in a klver board with the arm unrestrained .
in general , lfps showed less rhythmicity than during isometric tracking , but there were nevertheless periods of pronounced low - frequency oscillation in m1 with a phase that varied systematically across electrodes ( figure  3c ) .
we did not measure kinematics during these complex whole - limb movements , but characterized behavior instead using emgs recorded from multiple hand and wrist muscles .
when projected onto the  pc plane determined from isometric task recordings , the first two lfp - pcs were coherent with rectified emg in the delta band ( figure  3d ) , suggesting a consistent relationship between lfp and muscle activity even during unrestrained reach - to - grasp .
figure  5a shows power spectra for a representative m1 lfp under the four behavioral states : isometric task performance , free reaching , natural sleep , and ketamine sedation .
awake isometric and naturalistic reaching behaviors are characterized by a peak in the beta band around 20  hz , whereas sleep and sedation recordings show increased power at low frequencies .
a clear peak in the delta band is seen during sedation ( but not sleep ) , whereas awake behaviors are associated with a broad distribution of power at low frequencies .
the event - triggered analysis used in figure  4 was not applicable to these naturalistic recordings because the sleep data lacked k complexes , and the timing of submovements could not be accurately determined for free movements .
therefore , we examined whether a conserved dynamical structure could be found in the correlation structure between multichannel lfps , and between lfps and spiking activity , under the four different behavioral states . in all analyses
, we used the same 2d projection of the lfp data , which was determined by pca of the low - pass - filtered isometric torque data ( lfp - pcs ) . despite differences in lfp power spectra ,
the low - frequency correlation structure between lfp - pcs was preserved under all four behavioral states ( figure  5b ) .
cross - correlation of both unfiltered and low - pass - filtered lfp - pcs revealed strong , consistent peaks and troughs separated by about 150  ms , corresponding to an oscillatory cycle of around 3  hz . this is not a trivial consequence of pca decomposition , because although pcs must be uncorrelated at zero lag , there is no reason why they should be strongly correlated at any other lag .
moreover , it is not trivial that lfps recorded under other behavioral states , when projected onto the pc axes determined from the isometric data , should exhibit the same correlation structure . indeed
, this analysis revealed a subtle but systematic difference in the frequency of oscillation ( determined from the interval between cross - correlation peak and trough ) . in both animals for which we recorded in  all four conditions ,
the frequency of correlation was highest ( 3  hz ) during isometric movements and natural sleep , and slightly lower ( 2.8  hz ) during free reaching and ketamine sedation ( figure  5c ) .
although small , this difference was individually significant in both animals ( one - factor anova ; monkey d : f3,16  = 11.4 , p  = 3   10 ; monkey s : f3,10  = 8.0 , p  = 0.005 ) . another way to visualize
the similarity in correlation structure is to plot the lfp trajectory over time in the pc plane .
movie s4 shows real - time lfp data recorded in the different conditions , alongside its pc projection ( note that the awake data are expanded 2-fold to compensate for the increase in slow - wave amplitude during sleep and sedation ) . in all cases
the lfp trajectory rotated in the same direction , with a frequency of around 3  hz .
we quantified the extent to which these lfp trajectories could be captured by a single first - order linear dynamical equation of the form(1)x(t)=a.x(t),where the time evolution of the first two lfp - pcs , x(t ) , is determined only by a 2   2 matrix , a , with a trace equal to zero .
solutions of equation  1 are closed elliptical trajectories with constant frequency and direction of rotation ( figure  s2c ) , similar to the real data ( figure  s2d ) .
the three free parameters of a ( the fourth is fixed by the trace constraint ) effectively determine the frequency , orientation , and eccentricity of trajectories .
therefore , the extent to which a single matrix a can describe the time evolution of lfp - pc trajectories provides a measure of the consistency of the underlying dynamics .
we quantified the fit over sessions of isometric task performance using the coefficient of determination ( equation  10 in experimental procedures ) and obtained an average ( sd ) cod of 0.20   0.08 ( n  = 13 sessions in three animals ; table s3 ) .
the quality of this fit is not a trivial consequence of the orthogonality of pcs , because most orthogonal signals can not be described by equation  1 .
when the same analysis steps ( low - pass filtering , mean referencing , pca , and model fitting ) were applied to equivalent lengths of white noise , the 95% percentile of the distribution of the resultant cod was only 0.0013 .
equivalent analysis of beta - band lfp data ( filtered between 10 and 30  hz ; figure  s2b ) yielded an average cod of only 0.06   0.05 , significantly worse than the low - frequency fit ( n  = 13 , t  = 6.8 , p  = 2   10 , paired t test ; figure  5d ) .
this was not due to an absence of signal at this frequency , because beta - frequency oscillation was evident in the raw signal ( figure  s2a ) and power spectrum ( figure  5a ) .
rather , the oscillation at this frequency comprised predominantly a single phase leading to a high proportion of variance in the first pc ( figure  s2e ) , whereas the second pc had no consistent phase relationship .
therefore , trajectories in the pc plane lacked rotational structure ( figure  s2d ) , and hence could not be described by first - order linear dynamics .
next , we tested how well the model that best described the isometric task data could explain lfps recorded under ketamine sedation .
we applied the best - fit parameters obtained from the task recordings to predict the time derivative of the sedation data using equation  1 and achieved a comparable cod of 0.20   0.10 .
this was significantly better than the generalization of the beta - band model , which failed to fit these frequencies in the sedation data ( mean cod  = 0.05   0.05 , n  = 13 , t  = 8.5 , p  = 2   10 , paired t test ) . for comparison ,
the model with parameters best fit to the sedated state explained the delta - band data only marginally better , with a cod of 0.26   0.13 ( figure  5d ; table s3 ) , whereas the best fit to the beta - band activity remained poor ( cod  = 0.02   0.03 ) .
similar results were obtained for the generalization of the isometric task model to data recorded during natural sleep ( figure  5e ; table s4 ) and free reaching ( figure  5f ; table s5 ) . in both cases ,
the model parameters that best fit the task data captured a significantly higher proportion of lfp dynamics in the delta band compared with the beta band , and the quality of the fit was only marginally improved by fitting model parameters to the corresponding behavioral state .
these analyses confirm the consistent correlation structure in  multichannel m1 lfp activity under all four behavioral states , albeit with a minor ( 10% ) reduction in frequency during free  reaching and sedation .
next , we examined whether similar slow lfp oscillations were also observed in ventral premotor cortex ( pmv ) , and how they related to neural activity in each area .
figure  6 compares m1 and pmv activity during a single trial of isometric task performance , and figures 7a7c show average data for an entire session aligned to the end of each successful trial .
firing rates in m1 ( figures 6b and 7b ) were highest during the rising torque phase , as the animal made multiple submovements to acquire peripheral targets .
by contrast , pmv firing rates were highest after the end of the trial , as the animal took a food reward with the ipsilateral limb .
this is consistent with greater bilateral tuning of premotor neurons ( hoshi and tanji , 2006 ) , as well as with a strong preference for object - grasping movements within the bank of the arcuate sulcus ( umilta et  al . , 2007 ) .
these distinct periods of high neuronal activity were each associated with low - frequency lfp activity within the same cortical area ( figure  6c ) . in both m1 and pmv , the low - pass - filtered lfp could be decomposed into two orthogonal components ( figure  6d ) .
submovements during the trial were phase locked to the m1 cycle ( figure  6e ) , but had no consistent relationship to the pmv lfp .
the areal velocity of the lfp - pc trajectory in m1 and pmv was maximal during periods of high neural activity in the same cortical area ( figures 6f and 7c ) .
lfp - pcs within each area exhibited a similar low - frequency correlation structure ( figure  6 g ) , indicating a consistent phase lag throughout the recording .
however , the oscillations in each area were largely independent of each other during the isometric task , occurring at different phases of the task . as a result ,
correlations between lfp - pcs across areas were weaker than within areas ( figure  6h ) .
by contrast , neurons in m1 and pmv were coactive during free reaching - to - grasp with the contralateral limb ( figure  s3 ) . in this case , slow oscillations in both areas were phase locked , leading to robust correlation between lfp - pcs across areas .
finally , neural activity in both m1 and pmv under sedation was synchronized to k complexes ( figures 7d and 7e ; figure  s4 ) .
each k complex was also associated with synchronous bursts of low - frequency oscillation , resulting in peaks of lfp - pc areal velocity in both areas ( figure  7f ) .
last , we examined how spiking activity was related to the phase of the low - frequency oscillation in each area .
figure  8a and movie s5 show sample spike - triggered average trajectories of m1 and pmv lfp - pcs for the same set of neurons recorded during the four behavioral states . during isometric task performance
, neurons showed greater locking to lfps within the same cortical area , as expected from the dissociation of activity patterns during different task phases ( figures 7b and 7c ) .
however , during free reaching , sleep , and sedation , spike activity in both m1 and pmv was associated with cyclical lfp - pc trajectories within and across areas , consistent with the synchronization of low - frequency rhythms under these conditions . to assess phase locking to m1 lfp cycles across different data sets , we rotated the m1 lfp - pc plane such that the peak speed during isometric submovements occurred at a phase of zero . across the population ,
m1 neurons were significantly phase locked to m1 lfp oscillations in all three animals ( figure  8b ) , with an average preferred phase that preceded peak speed , consistent with these neurons having a causal role in movement . during task performance , pmv neurons in two out of the three animals ( monkeys d and r ) did not show consistent locking to m1-lfp cycles , whereas in the third animal ( monkey s ) the distribution of preferred phases was significantly nonuniform but nevertheless broad relative to m1 neurons .
however , during free reaching , natural sleep and ketamine sedation , neurons in both m1 and pmv became synchronized and fired at a similar preferred phase of the m1 cycle in all three animals . in summary , the relationship between spiking activity and lfps suggests that each cortical area is governed by its own intrinsic dynamics , allowing distinct slow oscillations to emerge in m1 and pmv when those areas are individually active during different phases of the isometric task .
however , during free behavior involving coactivation of m1 and pmv , as well as during sleep and sedation , the slow oscillations become coupled across cortical areas .
we have described a common 3  hz correlation structure in lfp recordings during an isometric movement task , free reaching , natural sleep , and ketamine sedation .
individual lfps exhibited oscillatory activity , albeit of lower amplitude in the awake state ,  with a consistent distribution of phase across electrodes relative  to submovements and k complexes .
because this phase  distribution was preserved across behavioral states , the multielectrode lfp could be decomposed into two orthogonal components that evolved according to the same underlying dynamics during all behavioral conditions . during isometric task performance , m1 neurons fired at a consistent phase of the cortical cycle , and this modulation of descending drive led to a 3  hz submovement structure in muscle activity and movement kinematics .
a similar cycle was evident during free reaching movements , which also comprise submovements ( milner and ijaz , 1990 ; roitman et  al . , 2004 ) , although peripheral interactions with limb biomechanics and afferent feedback may lower the frequency and obscure the clear rhythmicity seen in isometric tasks . interestingly ,
the isometric task also revealed a dissociation of neural activity within m1 and pmv during different task phases , each associated with a distinct low - frequency oscillation .
however , during sleep and under sedation , these rhythms became coupled across different areas , which may explain the increased slow - wave amplitude seen in these brain states .
it has previously been thought that the frequency of submovements during visuomotor tracking was determined by feedback loop delays , because their rhythmicity is disrupted under conditions of absent or delayed visual feedback ( miall et  al .
, 1986 , 1993 ; mcauley et  al . , 1999 ; miall and jackson , 2006 ) .
however , our finding of a common oscillatory structure in the cortical lfp that is ( 1 ) coherent with movement speed and ( 2 ) present during sleep and sedation reveals an intrinsic periodicity in motor circuitry at the submovement frequency .
arm stiffness ( selen et  al . , 2006 ) , or learning novel visuomotor mappings ( sailer et  al . , 2005 )
one possibility is that the intrinsic dynamics are tuned appropriately for visuomotor control such that different phases of the cycle are associated with the various computations involved in planning and generating the next submovement based on feedback from the previous one .
in fact , adaptation to delayed visual feedback is extremely limited ( miall and jackson , 2006 ) , suggesting that the motor system may in fact be tuned only to a narrow range of naturally occurring loop delays . indeed , kunesch et  al .
( 1989 ) concluded that the temporal characteristics of manipulative hand movements requiring tactile feedback were determined not by ( shorter ) sensorimotor loop delays but instead by central neural mechanisms responsible for interpreting sensory inputs .
this would be consistent with a common intrinsic oscillator shaping the structure of feedback - controlled movements , irrespective of the feedback modality .
finally , it is interesting to note that during verbal articulation there is coherence between cortical signals and mouth emg at a frequency of 23  hz , which reflects the spontaneous rhythmicity of speech ( ruspantini et  al . , 2012 ) .
care must be taken when inferring neural substrates of lfp activity , because synaptic and intrinsic currents from multiple neuronal populations contribute to the extracellular field ( buzski et  al . , 2012 ) .
moreover , rotation in the pc plane does not require underlying oscillatory sources that are orthogonal , because any consistent phase difference , or a traveling wave appearing with a different phase on each electrode , could equally be decomposed into orthogonal components ( rubino et  al .
2009 ; nauhaus et  al . , 2009 ; ray and maunsell , 2011 ) .
importantly , the distribution of preferred phase for neural firing was narrow compared to the lfp ( figure  s5a ) .
moreover , during sleep and sedation , this phase was common to neurons in both m1 and pmv ( figure  8b ) .
this appears incompatible with a traveling wave , which would cause neurons at different locations to fire at different preferred phases of the global cycle .
( 2012 ) reported complex , multiphasic patterns of cortical activity that could be projected onto a plane using the jpca method to reveal consistent cycles with notable similarity to the lfp trajectories we describe here .
however , it is not clear from that study whether all phases of the cycle were represented equally , because the jpca method is again based on orthogonal projections of the neural activity .
( 2013 ) found that movement - related potentials were composed of multiple components with amplitude and latency that varied systematically across the cortical surface , even though recorded neurons tended to be maximally active around movement onset .
a parsimonious explanation of the consistent correlation structure we describe is that the multichannel lfp comprises a mixture of at least two underlying sources with a fixed time / phase delay ( figures s5b and s5c ) . if one source reflects ( relatively ) synchronous neural activity occurring around submovements and k complexes ,
one possibility is that neural activity at other phases is undersampled in our recordings , either because the neurons are located in a different area of cortex or a subcortical structure or have smaller soma size ( for example , inhibitory interneurons ) .
an alternative explanation is that the field potential associated with synchronous neural activity may be composed of multiple sources with different time courses .
these sources are cortical , because submovement- and k complex - related lfp oscillations underwent polarity reversal within the gray matter , and we speculate that they may reflect excitatory and inhibitory synaptic potentials contributing to the generation of low - frequency rhythms .
delta oscillations can arise in the thalamus due to low - threshold calcium currents active in the hyperpolarized state ( amzica et  al .
however , in the awake state , thalamic neurons are depolarized and generally fire in a tonic mode ( steriade and llins , 1988 ) , suggesting that the low - frequency dynamics we observe during behavior may relate to a cortical delta rhythm that has recently been characterized in slice preparations .
this rhythm originates from intrinsic bursting cells in layer v that activate a source of gabab - mediated inhibition ( carracedo et  al . ,
the slow kinetics of this g protein - coupled receptor lead to sustained hyperpolarizing currents that can be delayed by several hundred milliseconds relative to inhibitory cell activity .
these slow currents are observed in the lfp ( dine et  al . , 2014 ) , and might be expected to contribute a low - frequency component with a substantial phase lag relative to ionotropic currents .
nevertheless , occasional bursting has been reported in the thalamus in the awake state ( guido and weyand , 1995 ) , and the relative contributions of cortical and corticothalamic mechanisms in generating delta activity in  vivo during behavior and sleep remain an important area for further investigation .
low - frequency lfps have several practical advantages for bmis ( rickert et  al . , 2005 ; bansal et  al . , 2011 ; hwang and andersen , 2013 ) , but our understanding of how these signals arise and how best to extract information from them is limited .
we found that the areal velocity swept out by lfp trajectories was proportional to movement speed , and suggest that this may prove a useful feature to examine for bmi applications , as it is robust to sources of synchronous noise ( because correlated signals lead to radial trajectories ) . in 3d pc space , there was a slight variation in the axis of rotation for different directions of movement . in effect
, the first two pcs captured the lfp trajectory that was common across all submovements , whereas the third component reflected more subtle variation in the neuronal sources associated with different directions ( waldert et  al . , 2009 ) .
these observations suggest that understanding the lawful dynamics that generate low - frequency behaviors may inform and constrain the search for more sophisticated approaches to decoding kinematics from lfps . by examining the dynamics of motor cortex activity
, we can unite two previously distinct phenomena : the rhythmicity of submovements during isometric tracking and delta oscillations during sleep and under sedation . in both cases , cortical neurons fire at distinct phases of the same underlying 3  hz lfp cycle , and
we suggest that this intrinsic rhythmicity reflects an underlying organization of motor cortical circuits engaged in feedback control of movement .
experiments were approved by the local ethics committee and performed under appropriate uk home office licenses in accordance with the animals ( scientific procedures ) act 1986 .
three purpose - bred female rhesus macaques ( monkey r : 5 years old , 5  kg ; monkey d : 6 years old , 6.5  kg ; monkey s : 5 years old , 5.4  kg ) were trained to control a cursor by generating isometric flexion - extension ( vertical ) and radial - ulnar ( horizontal ) torque with the left wrist restrained in pronated posture to move to eight peripheral targets presented in a pseudorandomized center - out sequence on a computer monitor .
wrist torque was measured using a six - axis force / torque transducer ( nano25 ; ati industrial automation ) .
cursor position was expressed as percentage of the distance to screen edge , with 100% corresponding to a torque of 0.67 nm .
targets were centered at 70% of the distance to screen edge and had a diameter of 25% . after training
, we implanted emg electrodes onto forearm and hand muscles , tunneled subcutaneously to a connector on the head . in a separate surgery , two custom arrays of 12 moveable 50  m diameter tungsten microwires ( impedance 200 k at 1 khz ) were implanted into the right m1 and pmv ( jackson and fetz , 2007 ) .
head - free recordings were made using unity - gain headstages followed by wide - band amplification and sampling at 24.4 ksp / s ( sp , sample ) ( system 3 ; tucker - davis technologies ) .
lfps were digitally low - pass filtered at 300  hz and recorded at 488 sp / s .
emgs were amplified ( 1,000 ) and band - pass filtered between 10 and 5,000  hz ( model 1700 ; am systems ) before sampling at 12.2 ksp / s .
kinematic decoding was performed on 13 sessions ( monkey r : 4 ; monkey d : 6 ; monkey s : 3 ) .
sedation data ( at least 5  min per session ) were recorded on separate days after induction with ketamine ( 10  mg / kg , intramuscularly ; i.m . ) and medetomidine ( 0.02  mg / kg , i.m . ) .
we report 13 pairs of movement and sedation sessions ( separated by no more than 3  days ) for all animals ( monkey d : 5 ; monkey r : 4 ; monkey s : 4 ) . in monkeys
d and s , we collected natural sleep data at the end of behavioral sessions , and report four and five pairs of sessions , respectively . in monkeys d and s , we also collected data while animals retrieved food from a klver board , and report six and four such sessions . on average ,
remaining lfps recorded during isometric task performance were separated by area ( m1 and pmv ) and processed by low - pass filtering ( 5  hz , four - pole , zero - phase digital butterworth filter ) , mean referencing ( i.e. , subtraction of the mean lfp across electrodes within the same cortical area ) , and dimensionality reduction using standard pca .
the pc plane was oriented such that the predominant rotational structure during task performance was in the anticlockwise direction . in all analyses , lfps recorded during other behavioral states ( free reaching , sleep , and sedation ) were always projected into the same pc space obtained from the corresponding isometric task data set .
cursor speed was calculated as the derivative of the magnitude of the 2d torque vector , that is , the radial component of velocity with a positive sign for movements away from the center of the screen .
k complexes were identified from a single surface - referenced lfp channel as the peak of a positive deflection that exceeded 250  v .
coherence spectra were calculated between unfiltered cursor position and speed ; rectified emg and lfp - pcs used a 2,048-point rectangular window with no overlap . although the pc axes were determined from low - pass - filtered data , we used unfiltered lfp projections for coherence spectra so as to include frequencies above the filter cutoff . online , semisupervised spike classification used principal component feature extraction and k means clustering ( spikepac ; tucker - davis technologies ) .
firing - rate profiles for each neuron were calculated offline by binning spike events ( into 488 bins / s ) , low - pass filtering at 5  hz , and normalizing to zero mean and unity standard deviation across the entire recording .
submovements were binned into nine groups of equal size according to increasing peak cursor speed , or alternatively into six groups according to  submovement direction .
2d lfp - pc trajectories from 200  ms before to 200  ms after the midpoint of each submovement were averaged within each group and quantified using areal velocity ( areal ) and frequency of rotation ( f):(2)vareal(t)=12x(t)x(t)(3)|vareal(t)|=12|x(t)|2.2.f(t).here , x(t ) is a 2d or 3d vector of lfp - pcs at time t , x(t ) is its derivative with respect to time , and  denotes the vector cross - product . for each submovement group ,
we used the average 3d areal velocity vector ( equation  2 ) from 200  ms before to 200  ms after the time of peak speed to decode the kinematics of individual submovements with leave - one - out cross - validation , as follows.(1)we parameterized the relationship between areal velocity vi ( for submovement i ) and the speed and direction of that submovement ( si,i ) assuming ( 1 ) for a given direction of submovement , the areal velocity magnitude increased linearly with cursor speed , and ( 2 ) for a given speed of submovement , the orientation of areal velocity vector varied with submovement direction .
specifically,(4)vi = si.b(i),with the direction - dependent component composed of a fourier series,(5)b(i)=b0+b1cos(i)+b2sin(i)+b3cos(2i)+b4sin(2i).we included only terms up to 2i to prevent overfitting .
the 15 free parameters were obtained by linear regression over the entire data set , excluding one submovement that was used for cross - validation.(2)we assessed how well the model predicted the speed ( or direction ) of the excluded submovement from the magnitude ( or orientation ) of the areal velocity vector associated with that submovement .
the decoded submovement speed was proportional to the magnitude of the areal velocity vector,(6)decodedspeed(sj)=|vj||b0|,whereas the decoded submovement direction was that which minimized the angular deviation between actual and predicted areal velocity vectors , calculated by the vector dot product,(7)decodeddirection(j)=argmax(vj.b()|b()|).steps 1 and 2 were repeated with a different submovement excluded until the speed and direction of all n submovements had been estimated.(3)we quantified decoding performance using coefficients of determination,(8)cod(speed)=1j(sjsj)2jsj2(9)cod(direction)=1njcos(jj).(4)finally , we determined significance thresholds ( p  < 0.05 ) for cod values by repeating the entire procedure for 1,000 surrogate data sets in which either the speed or direction was shuffled across submovements .
we parameterized the relationship between areal velocity vi ( for submovement i ) and the speed and direction of that submovement ( si,i ) assuming ( 1 ) for a given direction of submovement , the areal velocity magnitude increased linearly with cursor speed , and ( 2 ) for a given speed of submovement , the orientation of areal velocity vector varied with submovement direction .
specifically,(4)vi = si.b(i),with the direction - dependent component composed of a fourier series,(5)b(i)=b0+b1cos(i)+b2sin(i)+b3cos(2i)+b4sin(2i).we included only terms up to 2i to prevent overfitting .
the 15 free parameters were obtained by linear regression over the entire data set , excluding one submovement that was used for cross - validation .
we assessed how well the model predicted the speed ( or direction ) of the excluded submovement from the magnitude ( or orientation ) of the areal velocity vector associated with that submovement .
the decoded submovement speed was proportional to the magnitude of the areal velocity vector,(6)decodedspeed(sj)=|vj||b0|,whereas the decoded submovement direction was that which minimized the angular deviation between actual and predicted areal velocity vectors , calculated by the vector dot product,(7)decodeddirection(j)=argmax(vj.b()|b()|).steps 1 and 2 were repeated with a different submovement excluded until the speed and direction of all n submovements had been estimated .
finally , we determined significance thresholds ( p  < 0.05 ) for cod values by repeating the entire procedure for 1,000 surrogate data sets in which either the speed or direction was shuffled across submovements .
the phase of each lfp relative to submovements and k complexes was determined at the time of the event from a hilbert transform of the event - triggered average .
the correlation between lfp phase relative to each event was tested using the circular - circular correlation coefficient available in the circstat toolbox ( berens , 2009 ) .
the simplest linear system with oscillatory dynamics is a 2d state vector x(t ) that evolves according to equation  1 . under the conditions trace ( a )  = 0 and det ( a ) > 0 , this system exhibits stable periodic solutions .
we regressed the time derivative , x , against 2d lfp - pcs , x , to find the three free parameters of a ( the fourth is fixed by the trace constraint ) . because the lfp - pcs during the isometric task are orthogonal
however , we did not impose this constraint during model fitting , because it will not necessarily hold when data from other conditions are projected onto the same pc axes .
we measured the quality of fit to the lfp - pc derivative x(t) in each case with a vector coefficient of determination defined by(10)cod=1|x(t)x(t)|2dt|x(t)|2dt , where integration was performed over the entire recording .
the same analysis was applied separately to lfp data that had been filtered into delta- ( 05  hz ) and beta- ( 1030  hz ) frequency bands using four - pole , zero - phase digital butterworth filters before mean referencing and pca .
m1 and pmv lfps were averaged from 200  ms before to 200  ms after each spike .
the lfp averages were then projected onto the pc axes determined from the isometric task data and normalized by the standard deviation of the firing rate . to allow comparison across different data sets ,
the m1 lfp - pc plane was rotated such that the average submovement - triggered trajectory had zero phase at the moment of peak cursor speed , and the phase of the spike firing was measured relative to this .
significant phase locking across the population of neurons was assessed using the rayleigh test for circular uniformity in the circstat toolbox .
 summaryupper - limb movements are often composed of regular submovements , and neural correlates of submovement frequencies between 1 and 4  hz have been found in the motor cortex . 
 the temporal profile of movements is usually assumed to be determined by extrinsic factors such as limb biomechanics and feedback delays , but another possibility is that an intrinsic rhythmicity contributes to low frequencies in behavior . 
 we used multielectrode recordings in monkeys performing an isometric movement task to reveal cyclic activity in primary motor cortex locked to submovements , and a distinct oscillation in premotor cortex . during ketamine sedation and natural sleep , cortical activity traversed similar cycles and became synchronized across areas . because the same cortical dynamics are coupled to submovements and also observed in the absence of behavior , we conclude that the motor networks controlling the upper limb exhibit an intrinsic periodicity at submovement frequencies that is reflected in the speed profile of movements . 
Introduction
Results
Discussion
Experimental Procedures