A Recurrent Neural Network that Produces EMG from Rhythmic Dynamics -‐ David Sussillo*, Mark Churchland^*, Matt Kaufman#* & Krishna Shenoy* *
-‐ Stanford University, ^ -‐ Columbia University, # -‐ Cold Spring Harbor
It remains an open question how the firing rates of neurons in motor cortex (M1) lead to the EMG activity that ultimately drives movement. Recently, Churchland et al. 1 reported that neural responses in monkey M1 exhibit a prominent quasi-‐rhythmic pattern during reaching, even though the reaches themselves are not rhythmic. They argued that M1 could be understood as “an engine of movement that uses lawful dynamics”, i.e., that M1 could be viewed as a dynamical system. A major question posed by their work is finding a concise set of equations for a dynamical system that uses rhythmic patterns to drive EMG. We approached this problem by training a nonlinear recurrent neural network (RNN)2 to generate the recorded EMG during the same reach tasks used in 1. Because feedback connections endow the system with the ability to change dynamically in time, RNNs are a natural class of models to use when studying cortical circuits. We trained the RNN to simultaneously generate the EMG activity recorded from three muscles (deltoid, pectoral, and biceps) for 27 ‘conditions’ (reach types). The network was provided with condition-‐specific static inputs as an initial condition, derived from the actual preparatory activity of recorded neurons (panel A and B). The RNN architecture consisted of a simulated M1 circuit (sM1, 150 neurons), which provided input to three separate spinal cord circuits (sSC1-‐3, 25 neurons each performing nonlinear filtering of sM1 drive input). There were only two constraints on the system during optimization: 1) successfully generate the EMG and 2) using regularization techniques, do so as simply as possible. After training the RNN, it generated EMG with normalized RMS of 0.04 (panel B). We examined the network dynamics and uncovered a remarkably simple system that showed similarities to M1 on the individual neuron level (panels C and D). Further, the sM1 circuit exhibited oscillatory dynamics as a major component of the network activity. These dynamics, in turn, drove the spinal circuits to generate the EMG. The dimensionality of sM1 activity during simulation of the plan and movement required 15 principal components (PCs) to capture 99% of the variance, in reasonable agreement with M1 data. The spinal cord circuits required 3-‐5 PC dimensions. We investigated the nature of the sM1 population dynamics by applying a recently-‐ developed technique for identifying dimensions containing dynamical structure, jPCA1. The dynamics in the 1st jPC plane were strongly oscillatory and explained 23% of the variance of the network activity (panels E for monkey J, panel F for the RNN that generated EMG of monkey J). These rotations were produced by dynamics in the RNN whose linear approximation – around a local fixed point – contained strongly oscillatory structure reflected by eigenvalues with a large imaginary component. In addition to the rotational dynamics, we found a strong component of the neural trajectory, roughly orthogonal to the jPC plane (80 degrees), which carried the trajectories into the rotation. This component was similar across all conditions, and is thus captured by the ‘cross-‐condition mean’. Panel G shows a cartoon from 1, illustrating the idea, while panel H shows data from the sM1 circuit visualized in the space spanned by the jPC plane and the cross-‐condition mean. In summary, these simulations provide an existence proof that a dynamical system, when appropriately seeded, can generate the EMG of multiple muscles. Crucially, the dynamics are simple and consist primarily of (1) rotational dynamics and (2) a cross-‐ condition mean that brings the trajectories near the region in phase space where the rotations occur. We emphasize that neither the similarities of the RNN units to M1 neurons, nor the oscillatory patterns were built into the system. 1. Churchland, M. M. et al. Nature (2012), 2. Sussillo, D. & Abbott, L. F. Neuron 63 (2009).
B
...
A
Plan Input
sM1
sSC1-3
EMG Output
2 (a.u.)
projection onto jPC
projection onto jPC 2 (a.u.)
2
D A - Network architecture. C Condition-dependent prepatory activity provides the initial conditions (ICs) for the RNN to dynamically generate the EMG output. The RNN has four parts, a simulated M1 and 3 simulated spinal cord circuits, one for each muscle. B - 4 out of the 27 example triplets: input (black), RNN output (blue), and target SUPPLEMENTARY INFOR EMG (orange). C - 4 example PSTHs from M1 from monkeys performing a reach task. The 27 conditions are color coded from green through black to red, based !"##$%&%'()*+,-./"*%,001 on the level of plan-period !"#$!"%''()*+',#!+-.#"#+#/&.+)#0"%)1' activity. D - 4 example PSTHs from "#'!%)'#+-/'+),#"#2+34+"#5%6()*+,.#+ the RNN sM1 circuit chosen to &%)2(,(%)+5#/)7++8.#+*%/9+%:+,.('+!"# highlight the similarity of dynam!"%''()*+',#!+-/'+,%+:%&0'+/99+:0",.#"+ ics between neurons in M1 and /)/94'#'+%)+2(5#)'(%)'+-.#"# sM1. E - Projection of monkey J3 ',"%)*94+/&"%''+&%)2(,(%)'7++8.#+!/)#9'+3 ARTICLERESEARCH data onto first jPC plane, from 1 . '.%-+.%-+,.#+&"%''$&%)2(,(%)+5#/)+-/' F - Projection of sM1 data onto -400 TARGET 400 200 GO -200 MOVE -400 TARGET 400 200 GO -200 600 600 "#5%6#2+;!"#$%"&<=+,.#+&%)'#>0#) first jPC plane. This jPC plane b ofMonkey such patterns could have appeared by accident or for trivial reasons. Monkey B A Monkey J3 '03,"/&,()*+,.#+&"%''$&%)2(,(%) explained 23% the variance of c %"&<=+/)2+!%''(39#+"#9/,(%)'.(!'+3#,-##) sM1 activity. G - Cartoon provided E F &"%''$&%)2(,(%)+5#/)+/)2+,.#+"%,/,(%)/9+ in 1 to provide an intuition for how To address this possibility, multiple ‘shuffle’ controls demonstrate M1 might organize the cross!/,,#")'+;'"!!"($%"&<7++)1+?("()*+"/,#+/'+ that jPCA does not find rotational structure when such structure is condition mean with respect to :0)&,(%)+%:+,(5#+:%"+%)#+#@/5!9#+)#0"%) not present (Supplementary Figs 2, 3 and Supplementary Movie 4). the rotations in the jPC plane. The ;5%)A#4+B$/""/4+2/,/'#,C+DEF+&%)2(,(%)' Similarly, rotations in the walking monkey were not erroneously cross-condition mean takes G/&.+*"##)H"#2+,"/+!9%,'+,.#+/6#"/* found when the monkey was stationary (Supplementary Movie 1). system trajectories from all ICs "/,#+:%"+%)#+&%)2(,(%)=+'./2#2 The fact that a single plane (two dimensions)2captures an average of together to the oscillatory region 28% of the total data variance is9#6#9+%:+!"#!/"/,%"4+/&,(6(,47++ notable, given the high dimensionality and back again. H - A phase space !9%,'+,.#+&"%''$&%)2(,(%)+5#/)+;,.#+5#/ the dimensions defined by PC of the data itself14. As a comparison, diagram of sM1 activity during ,.#+%,.#"+,"/'<7++I#:%"#+'03'#>0#),+/)/ and PC3 (which by definition capture the second- and third-most data prepatory and movement phases variance possible) together capture 29% of the total variance. Thus, ;,.#+/!!9(&/,(%)+%:+JKL+/)2+MJKL e Monkey J-array f Monkey N-array (a.u.) N across all 27 Monkey conditions. The projection onto jPC projection jPC (a.u.) the jPCA projection simplyonto captures &%)2(,(%)+5#/)+-/'+'03,"/&,#27 1 1 response patterns that were visualized subspace is spanned by always present in the top PCs, but werefrequencies difficult to see because they gonal planes that captured rotational structure at different ()2#!#)2#),94+:%"+#/&.+)#0"%)7 the two jPC vectors and an were axis aligned. In fact, there wereoftypically two or three ortho(Supplementary Fig. 4). not Together these captured 50–70% the total ,.#+'/5#+#@/5!9#+)#0"%)+/:,#"+5#/)+ G H additional vector that captures data variance. Thus, rotations are a dominant feature of the popu'03,"/&,(%)7++8.#+&"%''$&%)2(,(%) the variance of the first principal lation response. This was true for primary motor;4#99%-<+('+)%-+N#"%+/,+/99+,(5#'7 cortex and dorsal component of the cross-condition premotor cortex independently (Supplementary L!!9(&/,(%)+%:+MJKL+,%+,.#+5%)A#4+B Fig. 5). mean (axis in red). Colored circles 2/,/'#,=+-(,.+)"$*+'!%,-!.")$"/$!01$ show the ICs provided by the Rotations, kinematics and EMG -")2.!.")$(1,)+;.313=+,.#+!"#$!"%''()*+ prepatory input to the RNN. Projection onto jPC Traditional views posit that motor cortex neurons are tuned for Projection onto jPC 1 Projection onto jPC 1 1 ,4'+-/'+)"!+/!!9(#2<7++8.('+!/)#9 During movement, the sM1 (a.u.) (a.u.) (a.u.) movement parameters such as direction. This perspective does not CCM &%),"/',#2+-(,.+,./,+()+:(*0"#+O dynamics move towards the jPC naturally account for the data in Fig. 3. 1We simulated neural populaProjections of the neural population response. a, Projection for ,#@,+;:%"+-.(&.+,.#+&"%''$&%)2(,(%) plane, yielding the rotational tions that were directionally tuned for velocity with an additional 74 neurons; 28 straight-reach conditions). Each trace (one jPC1 '03,"/&,#2<7++P.#)+,.#+&"%''$ 27 dynamics that ultimately drive the . Simulated preparatory activity non-directional sensitivity to speed lots the first 200 ms of movement-related activity away from the 28 )%,+'03,"/&,#2=+,.#+!"%M#&,(%)+%),%+,.#+:( simulated cord circuits to was tuned for reach direction and distance. We simulated one ‘velostate (circles). Traces are coloured on thespinal basis of the preparatoryjPC2 MJKL+!9/)#+&/!,0"#'+"#'!%)'#+!/,,#")'+,. produce their respective EMG. . a.u., arbitrary units.b, Projection for monkey A (64 city model’ data set per recorded data set, based upon the recorded 6/"4+%)94+-#/A94+/&"%''+&%)2(,(%)'7++8. straight-reach conditions)., Monkey J, data set 3 (55 neurons; 27 velocities and endpoints. Firing rates, trial counts, neuron counts and d, Monkey N (118 neurons; 27 straight- spiking )%,+'0"!"('()*Q+5/)4+)#0"%)'+2('!9/4+',"%)*+"#'!%)'#+:#/,0"#'+,./,+/"#+'(5(9/"+/&"%''+&%)2(,(%)'+;,.#+-#99 noise were matched to the recorded data. For velocity-model