Table Of ContentUsing Firing-Rate Dynamics to Train
Recurrent Networks of Spiking Model Neurons
6
1
0 BrianDePasquale,MarkM.Churchland1,L.F.Abbott1,2
2
n DepartmentofNeuroscience
a 1GrossmanCenterfortheStatisticsofMind
J
2DepartmentofPhysiologyandCellularBiophysics
8
2 ColumbiaUniversityCollegeofPhysiciansandSurgeons
NewYorkNY10032-2695USA
]
C
N
.
o
Abstract
i
b
-
q
[ Recurrent neural networks are powerful tools for understanding and modeling computa-
tion and representation by populations of neurons. Continuous-variable or “rate” model
1
v networks have been analyzed and applied extensively for these purposes. However, neu-
0
rons fire action potentials, and the discrete nature of spiking is an important feature of
2
6 neural circuit dynamics. Despite significant advances, training recurrently connected spik-
7
ing neural networks remains a challenge. We present a procedure for training recurrently
0
. connectedspikingnetworkstogeneratedynamicalpatternsautonomously,toproducecom-
1
0 plextemporaloutputsbasedonintegratingnetworkinput,andtomodelphysiologicaldata.
6 Our procedure makes use of a continuous-variable network to identify targets for train-
1
: ing the inputs to the spiking model neurons. Surprisingly, we are able to construct spiking
v
networks that duplicate tasks performed by continuous-variable networks with only a rel-
i
X
atively minor expansion in the number of neurons. Our approach provides a novel view
r
a of the significance and appropriate use of “firing rate” models, and it is a useful approach
forbuildingmodelspikingnetworksthatcanbeusedtoaddressimportantquestionsabout
representationandcomputationinneuralsystems.
Introduction
A fundamental riddle of nervous system function is the disparity between our continuous
and comparatively slow sensory percepts and motor actions and the neural representation
of those percepts and actions by brief, discrete and spatially distributed actions potentials.
Arelatedpuzzleisthereliabilitywithwhich thesesignalsarerepresenteddespitethevari-
ability of neural spiking across nominally identical performances of a behavior. A useful
approachtoaddressingtheseissuesistobuildspikingmodelnetworksthatperformrelevant
tasks, but this has proven difficult to do. Here we develop a method for constructing func-
tioningnetworksofspikingmodelneuronsthatperformavarietyoftaskswhileembodying
the variable character of neuronal activity. In this context, “task” refers to a computation
performedbyabiologicalneuralcircuit.
There have been previous successes constructing spiking networks that perform specific
tasks (see for example Seung et. al. (2000); Wang (2002); Machens, Romo & Brody
(2005); Hennequin et. al. (2014)). In addition, more general procedures have been devel-
oped (reviewed in Abbott, DePasquale & Memmesheimer (2016)) that construct spiking
networksthatduplicatesystemsoflinear(Eliasmith,2005;Boerlin&Dene`ve,2011;Boer-
lin, Machens & Dene`ve, 2013) and nonlinear (Eliasmith, 2005; Thalmeier et. al., 2015)
equations. However, most tasks of interest to neuroscientists, such as action choices based
onpresentedstimuli,arenotexpressedintermsofsystemsofdifferentialequations.
Our work uses continuous-variable network models (Sompolinsky, Crisanti & Sommers,
1988), typically called “rate” networks, as an intermediary between conventional task de-
scriptionsintermsofstimuliandresponses,andspikingnetworkconstruction.Thisresults
in a general procedure for constructing spiking networks that perform a wide variety of
tasks of interest to neuroscience (see also Thalmeier et. al. (2015); Abbott, DePasquale &
Memmesheimer (2016)). We apply this procedure to example tasks and show how con-
straints on the sparseness and sign (Dale’s law) of network connectivity can be imposed.
We also build a spiking network model that matches multiple features of data recorded
fromneuronsinmotorcortexandfromarmmusclesduringareachingtask.
Results
The focus of our work is the development of a procedure for constructing recurrently con-
nected networks of spiking model. We begin by describing the model-building procedure
andthenpresentexamplesofitsuse.
Networkarchitectureandnetworktraining
ThegeneralarchitectureweconsiderisarecurrentlyconnectednetworkofNleakyintegrate-
and-fire(LIF)modelneuronsthatreceivestask-specificinputF (t)and,followingtraining,
in
produces an approximation of a specified “target” output signal F (t) (Figure 1a). F (t)
out in
can be thought of as external sensory input or as input from another neural network, and
F (t) as the input current into a downstream neuron or as a more abstractly defined net-
out
2
a)
b)
Figure1:Networkarchitectures.a)Spikingnetwork.AnetworkofN recurrentlyconnectedleakyintegrate-
and-fire neurons (green circles) receives an input F (t) (grey circle) through synapses U, and generates an
in
output F (t) (red circle) through synapses W. Connections marked in red (recurrent connections J and
out
output connections W) are modified by training, and black connections are random and remain fixed, in-
cludingasecondsetofrecurrentconnectionswithstrengths Jf.b)Continuous-variablenetwork.Anetwork
of N˜ recurrently connected “rate” units (blue circles) receive inputs F (t) and F (t) through synapses U˜
in out
andu˜,respectively.Allconnectionsarerandomandheldfixed.ThesumofU˜F (t)andtherecurrentinput
out
determinedby J˜definestheauxiliarytargetsF (t)forthespikingnetwork.
J
work output (for example a motor output signal or a decision variable). The neurons in
the network are connected to each other by synapses with strengths denoted by the N × N
matrix J. Connections between the network and to the output have strengths given by an
N × N matrix W, where N is the number of outputs (either 1 or 2 in the examples we
out out
provide). During network training both J and W are modified. In addition to the trained
connectionsdescribedby J,wealsoincluderandomconnectionsdefinedbyanother N×N
matrix, Jf. The elements of this matrix are chosen randomly from a Gaussian distribution
with mean µ/N and variance g2/N and are not modified during training. The values of µ
f
and g are given below for the different examples we present. This random connectivity
f
produces chaotic spiking in the network (van Vreeswijk & Sompolinsky, 1998; Brunel,
2000),whichweuseasasourceofspikingirregularityandtrial-to-trialvariability.Weuse
theparameterg tocontrolthelevelofthisvariability.
f
When a neuron in the network fires an action potential, it contributes both fast and slow
synaptic currents to other network neurons. These currents are described by the two N-
dimensional vectors, s(t) and f(t). When neuron i in the network fires an action potential,
componentiofboth s(t)and f(t)isincrementedby1,otherwise
ds(t) df(t)
τ = −s(t) and τ = −f(t). (1)
s f
dt dt
3
Thetwotimeconstantsdeterminethedecaytimesoftheseslowandfastsynapticcurrents,
andwesetτ = 100msandτ = 2ms.Thesynapsesthataremodifiedbytrainingareallof
s f
the slow type, while the random synapses are fast. For example, the output of the network
istheproductof s(t)andtheoutputweightmatrixW,Ws(t)(Figure1a).
The membrane potentials of the model neurons, collected together into a N-component
vectorV,obeytheequation
dV (cid:16) (cid:17)
τ = V −V +g Js(t)+ Jff(t)+UF (t) +I, (2)
m rest in
dt
with τ = 20 ms. For a case with N inputs, U is an N × N matrix (we consider N = 1
m in in in
and 2) with elements drawn independently from a uniform distribution between -1 and 1.
I is a bias current set equal to 10 mV. It is increased between trials in the examples of
Figures3and4,representinga“holding”input.Eachneuronfiresanactionpotentialwhen
its membrane potential reaches a threshold V = −55 mV and is then reset to V =
th reset
V = −65 mV. Following an action potential, the membrane potential is held at the reset
rest
potential for a refractory period of 2 ms unless stated otherwise. The parameter g controls
thestrengthoftheinputstoeachneuron,andweprovideitsvalueforthedifferentexamples
below.
We can now specify the goal and associated challenges of network training. The goal is
to modify the entries of J andW so that the network performs the task specified by F (t)
in
and F (t),meaningthat
out
Ws(t) ≈ F (t) (3)
out
when the network responds to F (t) (with the approximation being as accurate as possi-
in
ble). Equation 3 stipulates that s(t) must provide a basis for the function F (t). If it does,
out
it is straightforward to compute the optimal W by minimizing the squared difference be-
tween the two sides of equation 3, averaged over time. This can be done either recursively
(Haykin,2002)orusingastandardbatchleast-squaresapproach.
Determiningtheoptimal J issignificantlymorechallengingbecauseoftherecurrentnature
ofthenetwork. J mustbechosensothattheinputtothenetworkneurons, Js(t),generates
a pattern of spiking that produces s(t). The circularity of this constraint is what makes
recurrentnetworklearningdifficult.Thedifferencebetweentheeasyproblemofcomputing
W and the difficult problem of computing J is that, in the case of W, we have the target
F (t) in equation 3 that specifies what signal W should produce. For J, it is not obvious
out
whattheinputitgeneratesshouldbe.
Suppose that we did have targets analogous to F (t) but for computing J (we call them
out
auxiliary target functions and denote them by the N-component vector F (t)). Then, J,
J
likeW, could be determined by a least-squares procedure, that is, by minimizing the time-
averagedsquareddifferencesbetweenthetwosidesof
Js(t) ≈ F (t). (4)
J
4
Therearestabilityissuesassociatedwiththisprocedure,thatwediscussbelow,howeverthe
main challenge in this approach is to determine the appropriate auxiliary target functions.
Our solution to this problem is to obtain them from a continuous-variable model. More
generally,ifwecantrainorotherwiseidentifyanothermodelthatimplementsasolutionto
atask,wecanusesignalsgeneratedfromthatmodeltotrainourspikingnetwork.
Usingcontinuousvariablemodelstodetermineauxiliarytargetfunctions
Equations 4 and 3, respectively, summarize two key features of the vector of functions
F (t): 1) They should correspond to the inputs of a recurrently connected dynamic sys-
J
tem, and 2) They should provide a basis for the network output F (t). To satisfy the
out
first of these requirements, we identify F (t) with the inputs of a recurrently connected
J
continuous-variable “rate” network. These networks have been studied intensely (Som-
polinsky, Crisanti & Sommers, 1988; Rajan, Abbott & Sompolinsky, 2010) and have been
trained to perform a variety of tasks (Jaeger & Haas, 2004; Sussillo & Abbott, 2009; Laje
& Buonomano, 2013; Sussillo, 2014). To satisfy the second condition, we use the desired
spikingnetworkoutput, F (t),asaninput totheratenetwork.Thisallowsustoobtainthe
out
auxiliarytargetfunctionswithouthavingtotrainthecontinuousvariablenetwork.
The continuous-variable model we use is a randomly connected network of N˜ firing-rate
units(throughoutweusetildestodenotequantitiesassociatedwiththecontinuous-variable
network). Like the spiking networks, these units receive the input F (t) and, as mentioned
in
above, they also receive F (t) as an input. The continuous-variable model is described by
out
an N˜-componentvector x(t)thatsatisfiestheequation
dx(t)
τ = −x(t)+g˜J˜H(x(t))+u˜F (t)+U˜F (t), (5)
x out in
dt
where τ = 10 ms, H is a nonlinear function (we use H(·) = tanh(·)), and J˜, u˜, and U˜ are
x
matrices of dimension N˜ × N˜, N˜ × N and N˜ × N , respectively. The elements of these
out in
matricesarechosenindependentlyfromaGaussiandistributionofzeromeanandvariance
1/N for J,andauniformdistributionbetween-1and1foru˜ andU˜ unlessstatedotherwise.
Wesetg˜ = 1.2exceptwhereidentifiedotherwise.
To be sure that signals from this driven network are appropriate for training the spiking
model, the continuous-variable network, driven by the target output, should be capable of
producing a good approximation of F (t). To check this, we can test whether an N × N˜
out out
matrixcanbefound(byleastsquares)thatsatisfiesW˜ H(x(t)) ≈ F (t)toasufficientdegree
out
of accuracy. Provided J˜ and u˜ are appropriately scaled, this can be done for a wide range
oftasks(Sussillo,2014).
The auxiliary target functions F (t) that we seek are generated from the inputs to the neu-
J
rons in the continuous-variable network. There is often, however, a mismatch between the
dimensionsof F (t),whichis N,andoftheinputstothecontinuous-variablemodel,which
J
isN˜.Todealwiththis,weintroduceanN×N˜ matrixu,withelementsdrawnindependently
5
(cid:112)
fromauniformdistributionovertherange± 3/N˜,andwrite
(cid:16) (cid:17)
F (t) = u g˜J˜H(x(t))+u˜F (t) . (6)
J out
We leave out the input term proportion to F (t) in this expression because the spiking
in
network receives the input F (t) directly. This set of target functions satisfies both of the
in
requirements listed at the beginning of this section and, as we show in the following ex-
amples,theyallowfunctionalspikingnetworkstobeconstructedbyfindingconnections J
that satisfy equation 4. We do this initially by a recursive least squares algorithm (Haykin,
2002),butlaterwediscusssolvingthisproblembybatchleastsquaresinstead.
Examplesoftrainednetworks
Theproceduredescribedabovecanbeusedtoconstructnetworksthatperformavarietyof
tasks. We present three examples that range from tasks inspired by problems of relevance
toneurosciencetomodelingexperimentaldata.
Our first example is an autonomous oscillation task that requires the network to generate a
self-sustained,temporallycomplexoutput(Figure2).F (t)forthistaskisaperiodicfunc-
out
tion created by summing sine functions with frequencies of 1, 2, 3, and 5 Hz. We require
that the network generates this output autonomously, therefore F (t) = 0 for this exam-
in
ple.Complex,autonomousoscillatorydynamicsareafeatureofneuralcircuitsinvolvedin
repetitivemotoractssuchlocomotion(Marder,2000).
Initially J = 0, so the activity of the network is determined by the random synaptic input
prodvide by Jf, and the neurons exhibit irregular spiking (Figure 2a-c). In this initial con-
figuration, the average Fano factor, computed using 100 ms bins, is 0.5, and the average
firingrateacrossthenetworkis5Hz.Followingthetrainingprocedure,thelearnedpostsy-
naptic currents Js(t) closely match their respective auxiliary target functions (Figure 2d),
and the network output similarly matches the target F (t) (Figure 2e). Residual chaotic
out
spiking due to Jf (Figure 2c) and the fact that we are approximating a continuous function
byasumofdiscontinuousfunctionscauseunavoidabledeviations.Nevertheless,anetwork
of3,000LIFneuronsfiringatanaveragerateof6.5HzwithanaverageFanofactorof0.25
performs this task with normalized post-training error of 5% (this error is the variance of
thedifferencebetweenWs(t)and F (t)dividedbythevarianceof F (t)).
out out
Because the output for this first task can be produced by a linear dynamical system, previ-
ous methods could also have been used to construct a functioning spiking network (Elia-
smith, 2005; Boerlin, Machens & Dene`ve, 2013). However, this is no longer true for the
following examples. In addition, it is worth noting that the network we have constructed
generates its output as an isolated periodic attractor of a nonlinear dynamical system. The
otherprocedures,inparticularthatofBoerlin,Machens&Dene`ve(2013),createnetworks
that reproduce the linear dynamics that generates F (t). This results in a system that can
out
produce not only Ws(t) ≈ F (t), but also Ws(t) ≈ αF (t) over a continuous range of α
out out
values.ThisoftenresultsinaslowdriftintheamplitudeofWs(t)overtime.Thepointhere
6
pre-training post-training
a
V
m
0
1
b
ns
o
ur
e
n
0
0
2
c
V
m
0
4
d
V
m
0
1
e
1
300 ms 300 ms
Figure2:Networkactivitybeforeandaftertrainingforanautonomousoscillationtask.(a)Membranevolt-
age and spiking activity of two example neurons. (b) Raster plot of 200 neurons. (c) Random recurrent
input Jff(t) for two example neurons. (d) Auxiliary target function F (t) (black) and learned recurrent in-
J
put Js(t) (red) for two example neurons. (e) Target output F (t) (black) and the generated output Ws(t)
out
(red)overoneperiod.Parameterforthisexample:N =3000,g=7mV,µ=−57,g =17,andN˜ =1000.
f
is that our procedure solves a different problem than previous procedures, despite the fact
that it generates the same output. The previous procedures were designed to duplicate the
lineardynamicsthatproduce F (t),whereasourprocedureduplicates F (t)uniquely.
out out
ThesecondtaskwepresentisatemporalXORtaskthatrequiresthenetworktocategorize
theinputitreceivesonagiventrialandreportthisdecisionthroughitsoutput.Eachtrialfor
thistaskconsistsofasequenceoftwopulsesappearingasthenetworkinput F (t)(Figure
in
3).Eachpulsehasanamplitudeof0.3,anditsdurationcanbeeithershort(100ms)orlong
(300 ms). The two pulses are separated by 300 ms, and after an additional 300 ms delay
the network must respond with either a positive or a negative pulse (with a shape given by
1/2 cycle of a 1 Hz sinewave). The rule for choosing a positive or negative output is an
exclusive OR function of the input sequence (short-short → −, short-long → +, long-short
→ +, long-long → −). The time between trials and the input sequence on each trial are
7
{100 ms, 300 ms} {300 ms, 300 ms}
{100 ms, 100 ms} {300 ms, 100 ms}
0 time (ms) 2100 0 time (ms) 2100
Figure 3: Temporal XOR task. The input F (t) (black) consists of two pulses that are either short or long
in
induration.TheoutputWs(t)(red)shouldreportanXORfunctionofthecombinationofpulses.Membrane
potentials of 10 example neurons (blue) are shown for the 4 different task conditions. Parameters for this
example:N =3000,g=10mV,µ=−40,g =12andN˜ =1000.
f
chosenrandomly.
A network of 3,000 LIF neurons with an average firing rate of 7 Hz can perform this task
correctlyon95%oftrials.Asintheautonomousoscillationtask,individualneuronspiking
activityvariesfromtrial-to-trialduetotheeffectof Jf.TheFanofactorcomputedacrossall
neurons,allanalysistimes,andalltaskconditionsis0.26.Thistaskrequiresintegrationof
eachinputpulse,storageofamemoryofthefirstpulseatleastuntilthetimeofthesecond
pulse, memory of the decision during the delay period before the output is produced, and
classificationaccordingtotheXORrule.
GeneratingEMGactivityduringreaching
Wenowturntoanexamplebasedondatafromanexperimentalstudy,withthespikingnet-
work generating outputs that match electromyograms (EMGs) recorded in 2 arm muscles
8
of a non-human primate performing a reaching task (Churchland, Lara & Cunningham,
2016).Inthistask,atrialbeginswiththeappearanceofatargetcueatoneofeightpossible
reach directions (task conditions). After a short delay, during which the arm position must
be held fixed, a “go” cue appears, instructing a reach to the target. The time between trials
andthesequenceofreachdirectionsarevariedrandomly.
To convey information about target location to the model, we use two network inputs de-
noted by a two-component vector F (t) and with amplitudes 2cos(θ) and 2sin(θ) where
in
the angle θ specifies the reach direction (Figure 4a left). The input is applied for 500 ms
and,whenitterminates,thenetworkisrequiredtogeneratetwooutputs(thus F (t)isalso
out
two-dimensional) that match trial-averaged and smoothed EMG recordings from the ante-
rior and posterior deltoid muscles during reaches to the specified target (Figure 4a right &
e).
A network of 5,000 neurons with an average firing rate of 6 Hz performs this task with a
normalized post-training error of 7% (Figure 4e), consistent with another modeling study
that used a firing-rate network (Sussillo et. al., 2015). The activity of the trained network
exhibits several features consistent with recordings from neurons in motor cortex. Individ-
ual neurons show a large amount of spiking irregularity that is variable across trials and
conditions (Figure 4b). The Fano factor computed across all neurons and all task condi-
tions drops during task execution (Figure 4d), consistent with observations across a wide
range of cortical areas (Churchland et. al., 2010). This network shows that EMGs can be
generatedbyanetworkwitharealisticdegreeofspikingvariability.
Another interesting feature of the network activity is the variety of different neural re-
sponses. Individual neurons display tuning during the input period, the output period, or
acrossmultipleepochswithdifferenttunings(Figure4f).Asinrecordingsfrommotorcor-
tex, consistently tuned neurons represent a minority of the network; most neurons change
theirtuningduringthetask(Figure4g).
To examine the population dynamics of this model, we performed PCA on filtered (4 ms
Gaussianfilter)single-trialspiketrainsfromaT×NCdatamatrix,whereT isthenumberof
timessampledduringatrial(2900),N isthenumberofneurons(5000),andCisthenumber
of reach conditions (8). We computed the eigenvectors of the T × T covariance matrix
obtained from these “data”, generating temporal PCs (Figure 5a & b). Each temporal PC
representsatemporalfunctionthatisstronglyrepresentedacrossthepopulationandacross
all reach conditions, albeit in different ways across the population on a single trial for
each condition. Two important features emerge from this analysis. First, more prominent
single-trialPCs(Figure5a)arerelativelyconsistentfromtrial-to-trial,whilelessprominent
PCs(Figure5b)varymore.Second,theprominentPCsfluctuateonaslowertimescalethan
the less prominent PCs. These two features indicate that the more prominent PCs form the
basis of the network output on a single trial, while the less prominent PCs reflect network
noise. This can be verified by reconstructing the network output using increasing numbers
of temporal PCs and calculating the fraction of the output variance captured (Figure 5c
9
a
e
b
f
c
z
H
0
3
d g1.0
1.0
C
F A
F
0.5 0.0
200 ms move onset 200 ms move onset
Figure 4: Producing EMG during reaching. (a) Task design. A two-dimensional input (left) is applied to
thenetworkfor500mstospecifyareachdirectionafterwhichthenetworkmustproduceoutputmatching
thecorrespondingEMGactivitypatternsrecordedfromtwoarmmuscles(right).Eachcolorrepresentsthe
activityforaspecificdirection.(b)Rasterplotshowingtheactivityofasingleneuronacrossalltrials(each
row)forallconditions(differentcolors).TheFanofactorforthisneuronis1.2.(c)Firingrateoftheneuron
shownin(b).Eachcolorrepresentsthetrial-averagedfiringrateforasinglecondition.(d)TheFanofactoras
afunctionoftimecomputedacrossallneuronsandconditions.(e)Ws(t)forbothoutputsandallconditions
(differentcolors)onasingletrial.(f)Firingratesforfournetworkneurons.Someneuronsaretunedduring
input or output periods exclusively (bottom two plots), while most are tuned during both periods (top two
plots). (g) Trial-averaged firing rate autocorrelation, averaged across all neurons and all conditions. By the
time of movement, the autocorrelation function is near zero, indicating that the tuning between input and
movement periods is, on average, uncorrelated. The time bar in all panels represents 200 ms, and the dot
denotesmovementonset.Parametersforthisexample:N =5000,µ=−112,g =33,g=7.2mV,N˜ =1000,
f
g˜ =1.4.Also,inthisexample,theelementsofu˜ werechosenrandomlyanduniformlyovertherange±0.76,
and the range for U˜ was ±0.25. The refractory period was set to zero in this example. EMGs were filtered
witha4msGaussianwindow.
10
