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