Predicting non-linear dynamics by stable local learning in
a recurrent spiking neural network
Aditya Gilra1,*

arXiv:1702.06463v2 [q-bio.NC] 26 Apr 2017

1

Wulfram Gerstner1

School of Computer and Communication Sciences, and Brain-Mind Institute, School of Life
Sciences, École Polytechnique Fédérale de Lausanne (EPFL), Lausanne 1015, Switzerland.
∗
Correspondence: aditya.gilra@epfl.ch

Abstract
Brains need to predict how the body reacts to motor commands. It is an open question
how networks of spiking neurons can learn to reproduce the non-linear body dynamics
caused by motor commands, using local, online and stable learning rules. Here, we present
a supervised learning scheme for the feedforward and recurrent connections in a network of
heterogeneous spiking neurons. The error in the output is fed back through fixed random
connections with a negative gain, causing the network to follow the desired dynamics, while
an online and local rule changes the weights. The rule for Feedback-based Online Local
Learning Of Weights (FOLLOW) is local in the sense that weight changes depend on the
presynaptic activity and the error signal projected onto the postsynaptic neuron. We provide examples of learning linear, non-linear and chaotic dynamics, as well as the dynamics
of a two-link arm. Using the Lyapunov method, and under reasonable assumptions and
approximations, we show that FOLLOW learning is stable uniformly, with the error going
to zero asymptotically.

1

Introduction

Over the course of life, we learn many motor tasks such as holding a pen, chopping vegetables,
riding a bike or playing tennis. To control and plan such movements, the brain must implicitly or
explicitly learn forward models [Conant and Ashby 1970] that predict how our body responds to
neural activity in brain areas known to be involved in motor control (Figure 1A). More precisely,
the brain must acquire a representation of the dynamical system formed by our muscles, our
body, and the outside world in a format that can be used to plan movements and initiate
corrective actions if the desired motor output is not achieved [Pouget and Snyder 2000; Wolpert
and Ghahramani 2000; Lalazar and Vaadia 2008]. Visual and / or proprioceptive feedback from
spontaneous movements during pre-natal [Khazipov et al. 2004] and post-natal development
[Petersson et al. 2003] or from voluntary movements during adulthood [Wong et al. 2012; Hilber
and Caston 2001] are important to learn how the body moves in response to neural motor
commands [Lalazar and Vaadia 2008; Wong et al. 2012; Sarlegna and Sainburg 2009; Dadarlat
et al. 2015], and how the world reacts to these movements [Davidson and Wolpert 2005; Zago
et al. 2005; 2009; Friston 2008]. We wondered whether a non-linear dynamical system, such as a
forward predictive model of a simplified arm, can be learned and represented in a heterogeneous
network of spiking neurons by adjusting the weights of recurrent connections.
Supervised learning of recurrent weights to predict or generate non-linear dynamics, given
command input, is known to be difficult in networks of rate units, and even more so in networks
1

of spiking neurons [Abbott et al. 2016]. Ideally, in order to be biologically plausible, a learning
rule must be online i.e. constantly incorporating new data, as opposed to batch learning where
weights are adjusted only after many examples have been seen; and local i.e. the quantities
that modify the weight of a synapse must be available locally at the synapse as opposed to
backpropagation through time (BPTT) [Rumelhart et al. 1986] or real-time recurrent learning
(RTRL) [Williams and Zipser 1989] which are non-local in time or in space, respectively [Pearlmutter 1995; Jaeger 2005]. Even though Long-Short-Term-Memory (LSTM) units [Hochreiter
and Schmidhuber 1997] avoid the vanishing gradient problem [Bengio et al. 1994; Hochreiter
et al. 2001] in recurrent networks, the corresponding learning rules are difficult to interpret
biologically.
Our approach toward learning of recurrent spiking networks is situated at the crossroads of
reservoir computing [Jaeger 2001; Maass et al. 2002; Legenstein et al. 2003; Maass and Markram
2004; Jaeger and Haas 2004; Joshi and Maass 2005; Legenstein and Maass 2007], FORCE learning
[Sussillo and Abbott 2009; 2012; DePasquale et al. 2016; Thalmeier et al. 2016; Nicola and
Clopath 2016], and adaptive control theory [Morse 1980; Narendra et al. 1980; Slotine and
Coetsee 1986; Slotine and Weiping Li 1987; Narendra and Annaswamy 1989; Sastry and Bodson
1989; Ioannou and Sun 2012]. In contrast to the original reservoir scheme [Jaeger 2001; Maass
et al. 2002] or neural network implementations of control theory [Sanner and Slotine 1992;
DeWolf et al. 2016] where learning was restricted to the read-out or feedforward connections, we
focus on a learning rule for the recurrent connections; and in contrast to FORCE learning where
recurrent synaptic weights have to change rapidly during the initial phase of learning [Sussillo
and Abbott 2009; 2012], we aim for a learning rule that works in the biologically more plausible
setting of slow synaptic changes. While previous work has shown that linear dynamical systems
can be represented and learned with local online rules in recurrent spiking networks [MacNeil and
Eliasmith 2011; Bourdoukan and Denève 2015], for non-linear dynamical systems the recurrent
weights in spiking networks have typically been computed offline [Eliasmith 2005].
Here, we propose a scheme for how a recurrently connected network of heterogeneous deterministic spiking neurons may learn to mimic a low-dimensional non-linear dynamical system,
with a local and online learning rule. The proposed learning rule is supervised, and requires
access to the error in observable outputs. The output errors are fed back with random, but fixed
feedback weights. Given a set of fixed error-feedback weights, the learning rule is synaptically
local and combines presynaptic activity with the local postsynaptic error variable.

2

Results

A forward predictive model (Fig. 1A) takes, at each time step, a motor command ~u(t)
as input and predicts the next observable state ~x̂(t + ∆t) of the system. In the numerical
implementation, we consider ∆t = 1ms, but for the sake of notational simplicity we drop the ∆t
in the following. The predicted system state ~x̂ (e.g., the position and velocity of the hand) is
assumed to be low-dimensional with dimensionality Nd (4-dimensional for a two-link arm). The
motor command ~u(t) is used to generate target movements such as “lift your arm to a location”,
with its dimensionality Nc typically smaller than the dimensionality Nd of the system state.
In our neural network model, the low dimensional motor command drives the spiking activity
of a command representation layer of 3,000 to 5,000 integrate-and-fire neurons (Fig. 1B) via
connections with fixed random weights. These neurons project, via plastic feedforward connections, to a recurrent network of also 3000 to 5000 integrate-and-fire neurons. We assume

2

A

B

C

D

Figure 1: Schematic for learning a forward model A. During learning, random motor commands (motor
babbling) cause movements of the arm, and are also sent to the forward predictive model, which must learn to
predict the positions and velocities (state variables) of the arm. The deviation of the predicted state from the
reference state, obtained by visual and proprioceptive feedback, is used to learn the forward predictive model
with architecture shown in B. B. Command input ~u projected onto neurons with random weights effkα . The spike
ff
trains of these command representation neurons Slff are sent via plastic feedforward weights wil
into the neurons
of the recurrent network having plastic weights wij (plastic weights in green). Readout weights dαi decode the
filtered spiking activity of the recurrent network as the predicted state x̂α (t). The deviations of the predicted
state from the reference state is fed back into the recurrent network (red and blue for different values of index
α) with encoding weights keiα . C. A cartoon depiction of feedforward, recurrent and error currents entering a
neuron i in the recurrent network. The error current enters the apical dendrite and triggers intra-cellular signals
available at the synapses for weight updates, isolated from the somatic current (here shown via basal dendrites).
D. A few spike trains of neurons of the recurrent network from the non-linear oscillator example are plotted.

3

that the predicted state x̂ is linearly decoded from the activity of the recurrent network. If we
denote the
P spikeR ttrain of neuron i by Si (t), then component α of the predicted system state is
x̂α (t) = i dαi −∞ Si (s)κ(t − s)ds where dαi are the readout weights. The integral represents
a convolution with a low-pass filter κ(t) with a time constant of 20 ms, and will be denoted by
(S ∗ κ)(t) in the following.
The predicted output is compared with the observable system state x(t), e.g., the position
and velocity of the hand deduced from visual and proprioceptive input. A nonlinear dynamical
control system will evolve generally as d~x/dt = ~h(~x(t), ~u(t)), but we first consider a simpler
dynamical system given by a set of coupled differential equations
dxα (t)
= fα (~x(t)) + gα (~u(t)),
dt

(1)

where ~x with components xα (where α = 1, . . . , Nd ) is the vector of observable state variables,
~u(t) ∈ RNc is the motor command input, and f~ and ~g are vectors whose components are arbitrary
non-linear functions fα and gα respectively.
Parameters of the spiking model neurons vary between one neuron and the next, both in
the command representation layer and the recurrent network, which yields different frequencycurrent curves for different neurons (Supplementary Fig. S1). Since arbitrary low-dimensional
functions can be approximated by linear decoding from a basis of nonlinear functions, such as
neuronal tuning curves [Funahashi 1989; Girosi and Poggio 1990; Hornik et al. 1989; Sanner and
Slotine 1992; Eliasmith and Anderson 2004], we may expect that a suitable synaptic plasticity
rule can tune the feedforward weights onto, and the lateral weights within, the recurrent network
so as to approximate the role of the functions ~g and f~ in equation (1), respectively, while the
read-out weights are kept fixed.
To enable weight tuning, we make four assumptions. First, we assume that, during the learning
phase, a random time-dependent motor command input ~u(t) is given to both the muscle-body
reference system described by equation 1 and to the spiking network. The random input produces
trajectories in the observable state variable corresponding to motor babbling [Petersson et al.
2003]. Second, we assume that each component x̂α of the output predicted by the spiking network
is compared to the actual observable output produced by the reference system of equation 1
and their difference α is calculated, similar to supervised learning schemes such as perceptron
learning [Rosenblatt 1961]. Third, we assume that this difference α = xα − x̂α is projected back
to neurons in the recurrent network through fixed random feedback weights with a large
P gain.
More precisely, our third assumption is that neuron i receives a total error input Ii = k α eiα α
with feedback weights keiα , where k is fixed at a large constant value. Fourth, we assume that
the readout weights dαi have been pre-learned by standard learning schemes [Voegtlin 2006;
Burbank 2015], so as to form an auto-encoder loop of gain k with the fixed random feedback
weights keiα , i.e. an arbitrary value α sent via the error feedback weights to the recurrent
network and read out, from its N neurons, via the decoding weights gives back (approximately)
kα . From adaptive control theory [Narendra and Annaswamy 1989; Ioannou and Sun 2012]
we may expect that the negative feedback arising from assumptions three and four, drives the
neurons to generate a coarse activity pattern that leads to a close-to-ideal observable output
(small errors), even at the very start of learning. Note that a vanishing error (α = 0 for all
components), after a sufficiently long learning time, indicates that the neuronal network has
autonomously generated the desired output so that feedback is no longer required.

4

While error feedback is on, the change in synaptic weights wilff and wij on the feedforward and
recurrent connections, respectively, is:
ẇilff = η (Ii ∗ κ )(Slff ∗ κ)(t),
ẇij = η (Ii ∗ κ )(Sj ∗ κ)(t),

(2)

where η is the learning rate, and κ is an exponentially decaying filter kernel with a time constant
of 200 ms. For a postsynaptic neuron i, the error term Ii ∗ κ is the same for all its synapses,
while the presynaptic contribution is synapse-specific. We call the learning scheme ‘Feedbackbased Online Local Learning Of Weights’ (FOLLOW), since the predicted state ~x̂ follows the
true state ~x from the start of learning. Under precise mathematical conditions, the FOLLOW
scheme converges to a stable solution (Methods subsections 4.4 and 4.5).
We emphasize that the learning rule of equation (2) uses an error α ≡ xα −x̂α in the observable
state rather than an error in the derivative that would appear if descending down the gradient
of a loss function (see Supplementary subsection 7.2) [Eliasmith 2005; MacNeil and Eliasmith
2011]. Furthermore, it is a local learning rule since all quantities needed on the right-hand-side
of equation 2 could be available at the location of the synapse in the postsynaptic neuron. For
a potential local implementation scheme, let us imagine that the postsynaptic error current
Ii arrives in the apical dendrite where it stimulates messenger molecules that quickly diffuse
or are actively transported into the soma and basal dendrites where synapses from feedfoward
and feedback input could be located, as depicted in Figure 1C. Consistent with the picture
of a messenger molecule, we low-pass filter the error current with an exponential filter κ of
time constant 200 ms, much longer than the synaptic time constant of 20 ms of the filter κ.
Simultaneously, filtered information about presynaptic spike arrival Sj ∗ κ is available at each
synapse, possibly in the form of glutamate bound to the postsynaptic receptor or by calcium
triggered signaling chains localized in the postsynaptic spines. Thus the combination of effects
caused by presynaptic spike arrival and error information available in the postsynaptic cell drives
weight changes, in loose analogy to standard Hebbian learning. For a critical evaluation of the
notion of ‘local rule’, we refer to the Discussion.

Spiking networks learn target dynamics via FOLLOW learning
In order to check whether the FOLLOW scheme would enable the network to learn various
dynamical systems, we studied three systems describing a nonlinear oscillator, low-dimensional
chaos and simulated arm movements (additional examples in Supplementary Figures S5, S6
and Methods). In all simulations, we started with vanishingly small feedforward and recurrent
weights (tabula rasa), but assumed pre-learned readout weights matched to the error feedback
weights. For each of the three dynamical systems, we had a learning phase and a testing
phase. During each phase, we provided time-varying input to both the network (Fig. 1B) and
the reference system. During the learning phase, rapidly changing control signals mimicked
spontaneous movements (motor babbling) while synaptic weights were updated according to the
FOLLOW learning rule (2).
During learning, the mean squared error, where the mean was taken over the number of
dynamical dimensions Nd and over a duration of a few seconds, decreased (Fig. 2D). We stopped
the learning phase i.e. weight updating, when the mean squared error approximately plateaued
as a function of learning time (Fig. 2D). At the end of the learning phase, we switched the error
feedback off (‘open loop’) and provided different test inputs to both the reference system and the
recurrent spiking network. A successful forward predictive model should be able to predict the
state variables in the open-loop model over a finite time horizon (corresponding to the planning
5

horizon of a short action sequence) and in the closed-loop mode (with error feedback) without
time limit.
Non-linear oscillator
Our FOLLOW learning scheme enabled a network with 3000 neurons in the recurrent network
and 3000 neurons in the motor command representation layer to approximate the non-linear 2dimensional van der Pol oscillator (Fig. 2). We used a superposition of random steps as input,
with amplitudes drawn uniformly from an interval, changing on two time scales 50 ms and 4 s
(see Methods).
During the 4 seconds before learning started, we blocked error feedback. Because of our
initialization with zero feedforward and recurrent weights, the output x̂ decoded from the network
of spiking neurons remained constant at zero while the reference system performed the desired
oscillations. Once the error feedback with large gain was turned on, the feedback forced the
network to roughly follow the reference. Thus, with feedback, the error dropped to a very low
value, immediately after the start of learning (Fig. 2B,C). During learning, the error dropped
even further over time (Fig. 2D). After having stopped learning at 50000 s, we found the weight
distribution to be uni-modal with a few very large weights (Fig. 2G). In the open-loop testing
phase without error feedback, a sharp square pulse as initial input on different 4 s long pedestal
values caused the network to track the reference as shown in Figure 2Aii-Cii panels. For some
values of the constant pedestal input, the phase of the output of the recurrent network differed
from that of the reference (Fig. 2Bii), but the shape of the non-linear oscillation was well
predicted as indicated by the similarity of the trajectories in state space (Fig. 2Cii).
The spiking pattern of neurons of the recurrent network changed as a function of time, with
inter-spike intervals of individual neurons correlated with the output, and varying over time (Fig.
2H,I). The distributions of firing rates averaged over a 250 ms period with fairly constant output,
and over a 16 s period with time-varying output, were long-tailed, with the mean across neurons
maintained at approximately 37 Hz (Fig. 2E,F). The latter distribution had lesser number of
very low- and very high-firing neurons compared to the former, consistent with the expectation
that the identity of low-firing and high-firing neurons changed over time for time-varying output
(Fig. 2E,F). We repeated this example experiment (‘van der Pol oscillator’) with a network of
equal size but with neurons that had firing rates increased by a factor of 2, so that some neurons
could reach a maximal rate of 400 Hz (Supplementary Fig. S1). We found that the learning
time was reduced by around 5 times (Supplementary Fig. S2) but with qualitatively the same
behaviour. Hence, for all further simulations, we set neuronal parameters to enable firing rates
up to 400 Hz (Supplementary Fig. S1B).
Chaotic Lorenz system
Our FOLLOW scheme also enabled a network with 5000 neurons each in the command representation layer and recurrent network, to learn the 3-dimensional non-linear chaotic Lorenz
system (Fig. 3). We considered a paradigm where the command input remained zero so that
the network had to learn the autonomous dynamics characterized in chaos theory as a ’strange
attractor’ [Lorenz 1963]. During the testing phase without error feedback minor differences led
to different trajectories of the network and the reference which show up as large fluctuations of
3 (t) (Fig. 3A-C). Such a behavior is to be expected for a chaotic system where small changes in
initial condition can lead to large changes in the trajectory. Importantly, however, the activity
of the spiking network exhibits qualitively the same underlying strange attractor dynamics, as
seen from the butterfly shape [Lorenz 1963] of the attractor in configuration space, and the tent
6

|

u2

Ai

0.20
0.00
0.15
feedback off | feedback on
5
0
5

| noise

end

Testing
Aii pulse on pedestal

feedback on |feedback off

feedback off

Bii

x2 , x̂ 2

0

time (s)

8

0

time (s)

50000

10 -2

49980

E 1500

time (s)

49988

rate (Hz)

200

time (s)

1.25

0
3

0

time (s)

8
¦

4

x1 , x̂ 1

4

0

rate (Hz)

200

1
0.50

time (s)

1.25

F1500

­

err 2

®

count

Nd , t

D

5
0
4

Cii 3

¦

count

Ci

error ²2

x2 ,x̂ 2

Bi

start

Learning

H

0

0

50

0.3 weight (arb·10−3 )

0.3

0

I

50
neuron #

45
density

G

0
neuron #

10 -6

1
0.50

Figure 2: Learning non-linear dynamics via FOLLOW: the van der Pol oscillator. A-C. Control input,
output, and error, before and at the start of learning; in the last 4 s of learning; and during testing without
error feedback (demarcated by the vertical red lines). Weight updating and error current feedback were both
turned on after the vertical red line on the left at the start of learning, and turned off after the vertical red
line in the middle at the end of learning. A. Second component of the input u2 . B. Second component of the
learned dynamical variable x̂2 (red) decoded from the network, and the reference x2 (blue). After the feedback
was turned on, the output tracked the reference. The output continued to track the reference, even after the end
of the learning phase, when feedback and learning were turned off. The output tracked the reference even with a
very different input (Bii). C. Second component of the error 2 = x2 − x̂2 between the reference and the output.
Cii Trajectory (x1 (t), x2 (t)) is the phase plane for reference (red,magenta) and prediction (blue,cyan) during two
different intervals as indicated by ? and  in Bii. D. Mean squared error per dimension averaged over 4 s blocks,
on a log scale, during learning with feedback on. E. Histogram of firing rates of neurons in the recurrent network
averaged over 0.25 s (interval marked in green in H) when output was fairly constant (mean across neurons was
36.95 Hz). F. As in E, but averaged over 16 s (mean across neurons was 37.36 Hz). G. Histogram of weights
after learning. A few strong weights |wij | > 0.3 are out of bounds and not shown here. H. Spike trains of 50
randomly-chosen neurons in the recurrent network (alternating colors for guidance of eye only). I. Spike trains
of H, reverse-sorted by first spike time after 0.5 s, with output component x̂2 overlaid for timing comparison.

7

u3

Ai

| startLearning end | zero

0.06

Testing
Aii

zero

0.00
0.06
feedback on |feedback off

feedback off

Bii

0
25

Cii 20

30
0
30

0

40

time (s)

14900

14940

time (s)

x1

x2

x̂3
x̂ 2

x̂1

0

time (s)

40

20

x1 , x̂ 1

25

max n (x3 )

20

0
15

E 20

x3

D

feedback off | feedback on

x2 , x̂ 2

error ²3

Ci

20

max n + 1 (x3 )

x3 ,x̂ 3

Bi

2

2

Figure 3: Learning chaotic dynamics via FOLLOW: the Lorenz system.
Layout and legend of panels A-C are analogous to Figure 2A-C. D. The trajectories of the reference (left panel)
and the learned network (right panel) are shown in state space for 40 s with zero input during the testing phase,
forming the well-known Lorenz attractor. E. Tent map, i.e. local maximum of the third component of the
reference signal (blue) / network output (red) is plotted versus the previous local maximum, for 800 s of testing
with zero input. The reference is plotted with filtering in panels A-C, but unfiltered for the strange attractor
(panel D left) and the tent map (panel E blue).

8

map [Lorenz 1963] of successive maxima versus the previous maxima (Fig. 3D,E). The tent
map generated from our network dynamics (Fig. 3E) has lower values for the larger maxima
compared to the reference tent map. However, very large outliers like those seen in a network
trained by FORCE [Thalmeier et al. 2016] are absent. Since we expected that the observed
differences are due to the filtering of the reference by an exponentially-decaying filter, we repeated learning without filtering the Lorenz reference signal (Supplementary Fig. S3), and found
that the mismatch is reduced, but a doubling appeared in the tent map which had been almost
imperceptible with filtering (Fig. 3E).

FOLLOW enables learning a two-link planar arm model under gravity
To turn to a task closer to real life, we next wondered if a spiking network can also learn the
dynamics of a two-link arm via the FOLLOW scheme. We used a two-link arm model adapted
from [Li 2006] as our reference. The two links in the model correspond to the upper and fore
arm, with the elbow joint in between and the shoulder joint at the top. The arm moved in
the vertical plane under gravity, while torques were applied directly at the two joints, so as to
coarsely mimic the action of muscles. To avoid full rotations, the two joints were constrained to
vary in the range from −90◦ to +90◦ where the resting state is at 0◦ (see Methods).
The dynamical system representing the arm is four-dimensional with the state variables being
the two joint angles and two angular velocities. The network must integrate the torques to
obtain the angular velocities which in turn must be integrated for the angles. Learning these
dynamics is difficult due to these sequential integrations involving nonlinear functions of the
state variables and the input. Because of the coupling of angles, torques, and moments, the
dynamics of the two-link arm belong to a more general class of nonlinear differential equations,
rather than equation (1). In particular, the general equations cannot be approximated by the
combined feedforward and recurrent network architecture used until now, as the control input
is not simply added to, but needs to be combined with the recurrent term. Thus, we use a
modified network architecture with only the recurrent network, as described in Supplementary
Section 7.3, and schematized in Supplementary Figure S4.
Similar to the previous examples, random input torque with amplitudes of short and long
pulses changing each 50 ms and 1 s, respectively, was provided to each joint during the learning
phase. The input was linearly interpolated between consecutive values drawn every 50 ms. In the
closed loop scenario with error feedback, the trajectory converges rapidly to the target trajectory
(Fig. 4). The general FOLLOW scheme learns to reproduce the arm dynamics even without
error feedback for a few seconds during the test phase (Fig. 4 and Supplementary Videos 1 and
2), which corresponds to the time horizon needed for the planning of short arm movements.
To assess the generalization capacity of the network, we fixed the parameters postlearning,
and tested the network in the open-loop setting on a reaching task and an acrobot-inspired
swinging task [Sutton 1996]. In the reaching task, torque was provided to both joints to enable
the arm-tip to reach beyond a specific (x, y) position from rest. The arm dynamics of the
reference model and the network are illustrated in Figure 4D and animated in Supplementary
Video 1. We also tested the learned network model of the 2-link arm on an acrobot-like task
i.e. a gymnast swinging on a high-bar [Sutton 1996], with the shoulder joint analogous to the
hands on the bar, and the elbow joint to the hips. The gymnast can only apply small torques at
the hip and none at the hands, and must reach beyond a specified (x, y) position, by swinging.
Thus, during the test, we provided input only at the elbow joint, with a time course that could
make the reference reach beyond a specific (x, y) position from rest by swinging. The control
9

u2

Ai

0.15
0.00
0.15
feedback off |feedback on
1.0
0.0
0.8
0.2
0.0

Testing
Aii acrobot-like task

feedback on |feedback off

feedback off

Bii

θ2 , θ̂ 2
ω2 , ω̂ 2

Bi

| start Learning end | noise

0.6

0

D

2
time (s)
reaching task

Cii 0.2 0
θ1 , θ̂ 1

error ²2

Ci

19995

time (s)

E

time (s)

0.0
0.2

19997
0.8
θ2 , θ̂ 2
acrobot-like task

2

1.0

0.0 s

0.3 s

0.6 s

gravity

reference
0.0 s

0.7 s

1.1 s

1.6 s

2.3 s

network

Figure 4: Learning arm dynamics via FOLLOW. Layout and legend of panels A-C are analogous to Figure
2A-C except that: in panel A, the control input (torque) on the elbow joint is plotted; in panel B, reference
and decoded angle θ2 , θ̂2 (solid) and angular velocity ω2 , ω̂2 (dotted) are plotted, for the elbow joint; in panel C,
the error θ2 − θ̂2 in the elbow angle is plotted. Aii-Cii. The control input was chosen to perform a swinging
acrobot-like task by applying small torque only on the elbow joint. Cii. The shoulder angle θ1 (t) is plotted
versus the elbow angle θ2 (t) for the reference (blue) and the network (red) for the full duration in Aii-Bii. The
green arrow shows the starting direction. D. Reaching task. Snapshots of the configuration of the arm, reference
in blue (top panels) and network in red (bottom panels) subject to torques in the directions shown by the circular
arrows. After 0.6 s, the tip of the forearm reaches the cyan target. Gravity acts downwards in the direction of
the arrow. E. Acrobot-inspired swinging task (visualization of panels of Aii-Cii). Analogous to D, except that
the torque is applied only at the elbow. To reach the target, the arm swings forward, back, and forward again.

10

A

C

10

test after 0 s of learning

10 -9

1

­

err 2

®

Nd , t

10 -2

1

1

10

0
R2

0
10

100000

R 2 = 0.925

0
1

time (s)

neuron number

learned weights (arb×10−4 )
recurrent
feedforward

B

0

10

1

10
1
10

= 0.990

0
10
−4
realizable network weights (arb×10 )

test after 10,000 s of learning

1
21.1

test after 100,000 s of learning

time (s)

21.2

Figure 5: Convergence of error, weights and spike times for a realizable reference network.
A. We ran our FOLLOW scheme on a network for learning one of two different implementations of the reference
van der Pol oscillator: (1) differential equations, versus (2) a network realized using FOLLOW learning for
10,000s. We plot the evolution of the mean squared error, mean over number of dimensions Nd and over 4 s
time blocks, from the start to 100,000 s of learning. With the weights starting from zero, mean squared error for
the differential equations reference (1) is shown in black, while that for the realizable network reference (2) is in
red. B. The feedforward weights (top panel) and the recurrent weights (bottom panel) at the end of 100,000s
of learning, are plotted versus the corresponding weights of the realizable target network. The coefficient of
determination i.e the R2 value of the fit to the identity line (y = x) is also displayed for each panel. A value of
R2 = 1 denotes perfect equality of weights to those of the realizable network. Some weights fall outside the plot
limits. C. After 0 s, 10,000 s, and 100,000 s of the learning protocol against the realizable network as reference,
we show spike trains of a few neurons in the recurrent network (red) and the reference network (blue) in the top,
middle and bottom panels respectively, from test simulations while providing the same control input and keeping
error feedback on.

input and the dynamics (Figure 4A-C right panels, Figure 4E and Supplementary Video 2) show
that the network can perform the task in open-loop condition suggesting that it has learned the
inertial properties of the arm model, necessary for this simplified acrobot task.

The FOLLOW scheme implicitly learns spike timings
In Methods subsections 4.4 and 4.5, we show that the FOLLOW learning scheme is Lyapunov
stable and that the error tends to zero under certain reasonable assumptions and approximations.
Two important assumptions of the proof are that the weights remain bounded and that the
desired dynamics are realizable by the network architecture, i.e. there exist feedforward and
recurrent weights that enable the network to mimic the reference dynamics perfectly. However,
in practice the realizability is limited by at least two constraints. First, even in networks of N
rate neurons with non-linear tuning curves, the nonlinear function f of the reference system in
equation (1) can in general only be approximated with a finite error [Funahashi 1989; Girosi
and Poggio 1990; Hornik et al. 1989; Sanner and Slotine 1992; Eliasmith and Anderson 2004]
which can be interpreted as a form of frozen noise, i.e. even with the best possible setting of the
weights, the network predicts, for most values of the state variables, a next state which is slightly
different than the one generated by the reference differential equation. Second, since we work
with spiking neurons, we expect on top of this frozen noise the effect of shot noise caused by
quasi-random spiking. Both noise sources may potentially cause drift of the weights [Narendra
and Annaswamy 1989; Ioannou and Sun 2012] which in turn can make the weights leave the
11

bounded regime. Ameliorative techniques from adaptive control are discussed in Supplementary
subsection 7.4. In our simulations, we did not find any effect of drift of weights on the error
during a learning time up to 100,000 s (Fig. 5A), 10 times longer than that required for learning
this example (Supplementary Fig. S2).
To highlight the difference between a realizable reference system and nonlinear differential
equations as a reference system, we used a spiking network with fixed weights as the reference.
For both the spiking reference network and the to-be-trained learning network we used the same
architecture, the same number of neurons, and the same neuronal parameters as in Supplementary Figure S2 for the learning of the van der Pol oscillator. However, instead of using the
differential equations of the van der Pol oscillator as a reference, we now used as a reference the
spiking network that was the final result after 10,000 s of FOLLOW learning in Supplementary
Figure S2, i.e., a spiking approximation of the van der Pol oscillator. The read-out and feedback
weights of the learning network had the same parameters as those of the spiking reference network, but the feedforward and recurrent weights of the learning network were initialized to zero
and updated, during the learning phase, with the FOLLOW rule. We ran FOLLOW learning
against the reference network for 100,000 s (Supplementary Fig. S2).
We emphasize that, analogous to the earlier simulations, the feedback error α was lowdimensional and calculated from the decoded outputs — as opposed to supervised learning
schemes for spiking neurons which are usually based on a direct comparison of spike times
[Gütig and Sompolinsky 2006; Pfister et al. 2006; Florian 2012; Mohemmed et al. 2012; Gütig
2014; Memmesheimer et al. 2014; Gardner and Grüning 2016]. We found that, with the realizable network as a reference, learning was more rapid than with the original van der Pol oscillator
as a reference. Interestingly, spike trains of the learning networks became similar to the spike
trains in the reference network, even though matching of spike times was not used as an optimization criterion (Fig. 5C). In particular, a few neurons fired only two or three spikes at very
precise moments in time. For example, after learning, the spikes of neuron i = 9 in the learning
network were tightly aligned with the spike times of the neuron with the same index i in the
spiking reference network. Similarly, neuron i = 8 that was inactive at the beginning of learning
was found to be active, and aligned with the spikes of the reference network, after 100,000 s of
learning.
Moreover, network weights became very similar, though not completely identical, to the
weights of the realizable reference network (Figure 5B), which suggests that the theorem for
convergence of parameters from adaptive control [Ioannou and Sun 2012; Narendra and Annaswamy 1989] should carry over to our learning scheme also.
Our results with the spiking reference network suggest that the error is reduced to a value
close to zero for a realizable or closely-approximated system (Methods subsection 4.5).

Learning is robust to sparse connectivity and noisy decoding weights
So far, our spiking networks had all-to-all connectivity. We next tested whether sparse connectivity [Markram et al. 2015; Brown and Hestrin 2009] of the feedforward and recurrent connections was sufficient for learning low-dimensional dynamics. We ran the van der Pol oscillator
learning protocol with the connectivity varying from 0.1 (10 percent connectivity) to 1 (full connectivity). Connections that were absent after the sparse initialization could not appear during
learning, while the existing sparse connections were allowed to evolve according to FOLLOW
learning. As shown in Figure 6A, we found that learning was slower with sparser connectivity;
12

B

10 -5

10 4 s

10 4 s, χ = 2

­

err 2

C

1 · 10 4 s
2 · 10 4 s

®

Nd , t

A 10-4

10 -6
0.0

0.5
connectivity

1.0

0

10
18
decoder noise param. χ

0.5

0.0
0.5
decoder noise param. ξ

Figure 6: Learning error versus sparse connectivity and noisy decoding weights.
We ran the van der Pol oscillator learning protocol for 10,000 s for different parameter values and measured the
mean squared error, over the last 400 s before the end of learning, mean over number of dimensions Nd and time.
A. We evolved only a fraction of the feedforward and recurrent connections, randomly chosen as per a specific
connectivity, according to FOLLOW learning, while keeping the rest zero. The round dots show the mean squared
error for different connectivity after a 10,000 s learning protocol (default connectivity = 1 is starred); while the
square dots show the same after a 20,000 s protocol. B. We multiplied the original decoding weights (that form
an auto-encoder with the error encoders) by a random factor (1+uniform(−χ, χ)) drawn for each weight. The
mean squared error at the end of a 10,000s learning protocol for increasing values of χ is plotted (default χ = 0 is
starred). C. We multiplied the original decoding weights by a random factor (1+uniform(−χ + ξ, χ + ξ)), fixing
χ = 2, drawn for each weight. The mean squared error at the end of a 10,000 s learning protocol, for a few values
of ξ on either side of zero, is plotted.

but with twice the learning time, a sparse network with about 25% connectivity reached similar
performance as the fully connected network with standard learning time.
The read-out weights have been pre-learned until now, so that, in the absence of recurrent
connections, error feedback weights and decoding weights formed an auto-encoder. We sought
to relax this requirement. Simulations showed that with completely random read-out weights,
the system did not learn to reproduce the target dynamical system. However, if the read-out
weights had some overlap with the auto-encoder, learning was still possible (Figure 6C). If for a
feedback error ~, the error encoding followed by output decoding yields k(1 + ξ)~ + n(~), where n
is an arbitrary function, and ξ is sufficiently greater than −1 so that the effective gain k(1 + ξ)
remains large enough, then the term linear in error can still drive the output close to the desired
one (see Methods).
To check this intuition in simulations, we incorporated multiplicative noise on the decoders
by multiplying each decoding weight of the auto-encoder by one plus γ, where for each weight γ
was drawn independently from a uniform distribution between −χ + ξ and χ + ξ. We found that
the system was still able to learn the van der Pol oscillator up to χ ∼ 5 and ξ = 0, or χ = 2 and
ξ variable (Figure 6B,C). Negative values of ξ result in a lower overlap with the auto-encoder
leading to the asymmetry seen in Figure 6C.
In conclusion, the FOLLOW learning scheme is robust to multiplicative noise on the decoding
weights. Alternatively, decoder noise can also be studied as frozen noise in approximating the
reference, causing a drift of the learned weights, possibly controlled by weight decay (Supplementary subsection 7.4).

3

Discussion

The FOLLOW learning scheme enables a spiking neural network to function as a forward
predictive model that mimics a non-linear dynamical system activated by one or several timevarying inputs. The learning rule is supervised, local, and comes with a proof of stability.
13

It is supervised because the FOLLOW learning scheme uses error feedback where the error
is defined as the difference between predicted output and the actual observed output. Error
feedback forces the output of the system to mimic the reference, an effect that is widely used in
adaptive control theory [Narendra and Annaswamy 1989; Ioannou and Sun 2012].
The learning rule is local in the sense that it combines information about presynaptic spike
arrival with an abstract quantity that we imagine to be available in the postsynaptic neuron.
In contrast to standard Hebbian learning, the variable representing this postsynaptic quantity
is not the postsynaptic firing rate, spike time, or postsynaptic membrane potential, but the
error current projected by feedback connections onto the postsynaptic neuron, similar in spirit
to modern biological implementation of approximated BackPropagation [Roelfsema and van
Ooyen 2005], [Lillicrap et al. 2016] or local versions of FORCE [Sussillo and Abbott 2009] learning
rules. We emphasize that the postsynaptic quantity is different from the postsynaptic membrane
potential or the total postsynaptic current which would also include input from feedforward and
recurrent connections. The separation of the error current from the currents at feedforward
and recurrent synapses could be spatial (such as suggested in Fig. 1C) or chemical if the error
current projects onto synapses that trigger a signaling cascade that is different from that at
other synapses. Importantly, whether it is a spatial or chemical separation, the signals triggered
by the error currents need to be available throughout the postsynaptic neuron.
The learning rule is provenly stable with errors converging asymptotically to zero under a few
assumptions (Methods subsection 4.4). The first assumption is that error encoding feedback
weights and output decoding read-out weights form an auto-encoder. This requirement can be
met if both sets of weights are learned at an early developmental stage, e.g., using mirrored
STDP [Burbank 2015]. The second assumption is that the reference dynamics f (~x) is realizable.
This requirement can be approximately met by having a recurrent network with a large number
N of neurons with different parameters [Eliasmith and Anderson 2004]. The third assumption
is that the state variables ~x(t) are observable. While currently we calculate the feedback error
directly from the state variables as a difference between reference and predicted state, we could
soften this condition and calculate the difference in a higher-dimensional space with variables ~y (t)
as long as ~y = K(~x) is an invertible function of ~x(t) (Supplementary section 7.3). The fourth
assumption is that the system dynamics be slower than synaptic dynamics. Indeed, typical
reaching movements extend over hundreds of millisconds or a few seconds whereas neuronal
spike transmission delays and synaptic time constants can be as short as a few milliseconds.
In our simulations, neuronal and synaptic time constants were set to 20 ms, yet the network
dynamics evolved on the time scale of hundreds of milliseconds or a few seconds, even in the
open-loop condition when error feedback was switched off (Figs. 2 and 4). The fifth assumption
is that weights stay bounded. Indeed, in biology, synaptic weights should not grow indefinitly.
Algorithmically, a weight decay term in the learning rule can suppress the growth of large weights
(see also Supplementary subsection 7.4), though we did not need to implement a weight decay
term in our simulations.
Simulations with the FOLLOW learning scheme have demonstrated that strongly nonlinear
dynamics can be learned in a recurrent spiking neural network using a local online learning rule
that does not require rapid weight changes. Previous work has mainly focused on a limited
subset of these aspects. For example, Eliasmith and colleagues used a learning rule derived from
stochastic gradient descent, in a network structure comprising heterogeneous spiking neurons
with error feedback [MacNeil and Eliasmith 2011], but did not demonstrate learning non-linear
dynamics (Supplementary subsection 7.2). The mathematical framework of Denève and colleagues with error feedback in a homogeneous spiking network with fast and slow connections
14

yields a rule similar to ours, but has been applied to linear dynamics only [Bourdoukan and
Denève 2015]. It will be interesting to see if their approach to learning the autoencoder and
maintaining balance could be used in our heterogeneous network.
Reservoir computing models exploit recurrent networks of non-linear units in an activity
regime close to chaos where temporal dynamics is rich [Jaeger 2001; Maass et al. 2002; Legenstein et al. 2003; Maass and Markram 2004; Jaeger and Haas 2004; Joshi and Maass 2005;
Legenstein and Maass 2007]. While typical applications of reservoir computing are concerned
with tasks involving a small set of desired output trajectories (such as switches or oscillators),
our FOLLOW learning enables a recurrent network with a single set of parameters to mimic a
dynamical system over a broad range of time-dependent inputs with a large family of different
trajectories in the output.
Whereas initial versions of reservoir computing focused on learning the readout weights, applications of FORCE learning to recurrent networks of rate units made it possible to also learn the
recurrent weights [Sussillo and Abbott 2009; 2012]. However, in the case of a multi-dimensional
target, multi-dimensional errors were typically fed to distinct parts of the network, as opposed to
the distributed encoding used in our network. Moreover, the time scale of plasticity in FORCE
learning is faster than the time scale of the dynamical system which is unlikely to be consistent
with biology. Modern applications of FORCE learning to spiking networks [DePasquale et al.
2016; Thalmeier et al. 2016; Nicola and Clopath 2016] inherit these issues.
Adaptive control of non-linear systems using continuous rate neurons [Sanner and Slotine
1992; Slotine and Weiping Li 1987; Slotine and Coetsee 1986] and or spiking neurons [DeWolf
et al. 2016] has primarily focused on learning the weights on a feedforward processing path, but
not on the weights in a recurrent network (note that adaptive control systems incorporate an
error feedback loop but in the current context, we distinguish between the error feedback loop
and the recurrent network). Recurrent networks of rate units have occasionally been used for
control [Zerkaoui et al. 2009], but trained either via real-time recurrent learning or the extended
Kalman filter which are non-local in space, or via backpropagation through time which is offline
[Pearlmutter 1995]. Optimal control methods [Hennequin et al. 2014] or stochastic gradient
descent [Song et al. 2016] have also been applied in recurrent networks of neurons, but with
limited biological plausibility of the published learning rules. As an alternative to supervised
schemes, biologically plausible forms of reward-modulated Hebbian rules on the output weights
of a reservoir have been used to learn periodic pattern generation and abstract computations
[Hoerzer et al. 2014; Legenstein et al. 2010], but how such modulated Hebbian rules could be
used in predicting nonlinear dynamics for time-dependent control input remains open.
We found that the FOLLOW learning scheme does not require full connectivity but also works
with biologically more plausible sparse connectivity. Furthermore, it is robust to multiplicative
noise in the output decoders, analogous to recent results on approximate error backpropagation
in artificial neural networks [Lillicrap et al. 2016]. Since the low-dimensional output and all
neural currents are spatially averaged over a large number of synaptically-filtered spike trains,
neurons in the FOLLOW network do not necessarily need to fire at rates higher than the inverse
of the synaptic time scale.
Our present implementation of the FOLLOW learning scheme in spiking neurons violates
Dale’s law because synapses originating from the same presynaptic neuron can have positive
or negative weights, but in a different context extensions incorporating Dale’s law have been
suggested [Parisien et al. 2008]. Neurons in cortical networks are also seen to maintain a balance
15

of excitatory and inhibitory incoming currents [Denève and Machens 2016]. It would be interesting to investigate a more biologically plausible extension of FOLLOW learning that maintains
Dale’s law; works in the regime of excitatory-inhibitory balance, possibly using inhibitory plasticity [Vogels et al. 2011]; pre-learns the autoencoder, potentially via mirrored STDP [Burbank
2015]; and possibly implements spatial separation between different compartments [Urbanczik
and Senn 2014]. Further directions worth pursuing include learning multiple different dynamical
transforms within one recurrent network, without interference; hierarchical learning with stacked
recurrent layers; and learning the inverse model of motor control so as to generate the control
input given a desired state trajectory.

4
4.1

Methods
Simulation software

All simulation scripts were written in python (https://www.python.org/) for the Nengo 2
simulator [Stewart et al. 2009] (http://www.nengo.ca/) with minor custom modifications to
support sparse weights. We ran the model using the Nengo GPU backend (https://github.
com/nengo/nengo_ocl) for speed. The script for plotting the figures was written in python
using the matplotlib module (http://matplotlib.org/). These simulation and plotting scripts
will be made available online at https://github.com/adityagilra, once this article is peerreviewed and published.

4.2

Network architecture and parameters

For all numerical simulations, we used deterministic leaky integrate and fire (LIF) neurons
[Eliasmith and Anderson 2004]. The voltage Vl of each LIF neuron indexed by l, was a low-pass
filter of its current Jl :
dVl
= −Vl + Jl ,
τm
dt
with τm the membrane time constant, set at 20 ms. The neuron fired when the voltage Vl
crossed a threshold θ = 1 from below, after which the voltage was reset to zero for a refractory
period τr of 2 ms. If the voltage went below zero, it was clipped to zero. The output spike train
of neuron l is denoted Slff (t) in the command representation layer and Sl (t) in the recurrent
network. Mathematically, a spike train is a sequence of events, modeled as a sum of Dirac
delta-functions. To get the synaptic input current, spike trains were filtered with an exponential
kernel κ(t) ≡ exp(−t/τs )/τs with a time constant of τs = 20 ms.
The low-dimensional motor command input ~u, with Nc components uα was projected onto a
layer of command-representing neurons. The current into a neuron with index l was
X
Jlff = νlff
effkα uα + bffl ,
(3)
α

where effkα were fixed random weights, while bffl and νlff were neuron-specific constants for bias
and gain respectively (see below). We use Greek letters for the indices of low-dimensional
variables and Latin letters for neuronal indices, with summations going over the full range of
the indices. The number of neurons N in the command-representation layer was much larger
than the dimensionality of the input i.e N  Nc .

16

A recurrent network of another N neurons received feedforward input from the commandrepresentation neurons, input from within the recurrent network, and error feedback. The input
current to a neuron with index i in the recurrent network was
!
X
X
X
ff
ff
Ji =
wil (Sl ∗ κ)(t) +
wij (Sj ∗ κ)(t) +
keiα (α ∗ κ)(t) νi + bi ,
(4)
α

j

l

where wilff and wij were the feedforward and recurrent weights, respectively, which were both
subject to our synaptic learning rule, whereas keiα were fixed error feedback weights (see below).
The spike trains travelling along the feedforward path Slff and those within the recurrent network
Sj were both low-pass filtered (convolution denoted by ∗) at the synapses with an exponential
filter κ with a time constant of 20 ms. The parameters bi and νi were neuron specific constants
for the bias and gain, respectively. The constant k > 0 was the feedback gain, and α ≡ xα − x̂α
was the output error. Here too, the number of neurons N in the recurrent network was much
larger than the dimensionality Nd of the represented variable x̂ i.e N  Nd .
Fixed, random encoding weights effkα and error feedback weights (without factor k) eiα were
uniformly chosen on an Nd -dimensional hypersphere (Nengo 2 default) of radius 1/R1 and 1/R2
respectively. The norm of the variables, say yα represented in the command representation
layer
P
(subscript 1) or the recurrent network (subscript 2) could range from (0, R1,2 ) for | α eiα yα | < 1,
where R1 and R2 (Table 1) will be called the representation radii of the two ensembles of N
neurons each. Typically R1 was smaller than (around 0.2 of) R2 as the input from the command
representation layer was integrated in the recurrent network.
P
The intercept of the transfer function ai ( α eiα yα ), i.e. firing rate vs. input curve, for each
neuron (Fig. S1), was chosen uniformly from
P (−1, 1); and its maximum firing rate for an input
within the representation radius, i.e. for α eiα yα = 1, was chosen uniformly from (100, 200)
Hz for the simulations in Figure 2 and from (200, 400) Hz for all other simulations. From this
desiredP
intercept and maximum rate for the transfer function of each LIF neuron, with input
zi = γi α eiα yα + bi , its bias bi and gain γi were calculated from the static rate equation for the
LIF neuron
ai = 1/(τr + τm ln(1 − 1/zi )).
(5)
Ideally, these encoders, biases and gains would be learnt during development from the input
statistics, but here we set them randomly.
Decoders for the output ~x̂, i.e. the linear readout weights dβi from the recurrently connected
network, were computed algorithmically to form an auto-encoder with respect to error-feedback
weights eiα (factor k on these weights gives us gain k). To do this, we randomly selected N error
(p)
vectors, that we used as training points for optimization, β where 1 ≤ p ≤ N is the label of
the training sample and β is the index of the vector component. Since the observable system
is Nd dimensional, we chose the training points randomly from an Nd -dimensional hypersphere
of radius R1 . We applied each of the error vectors statically as input for the error feedback
(p)
connections and calculated the activity ai of neuron i for error vector p using the static equation
(5). The decoders dβi acting on these activities should yield back the encoded points thus forming
an auto-encoder. A squared-error loss function
L=

P X X
N
X
p

β

i

17

!2
(p)

(p)

dβi ai − β

(6)

with L2 regularization was used for this linear regression (default in Nengo 2) [Eliasmith and
Anderson 2004]. Biologically plausible learning rules exist for auto-encoders [Burbank 2015;
Voegtlin 2006], but we simply calculated and set the decoding weights as if they had already
been learned.
Classical three-layer (input-hidden-output-layer) auto-encoders come in two different flavours,
viz. compressive or expansive, which have the dimensionality of the hidden layer smaller or
larger respectively, than that of the input and output layers. Instead of a three-layer feedfoward
network, our auto-encoder forms a loop from the neurons in the recurrent network via readout
weights to the output and from there via error-encoding weights to the input. Since the autoencoder is in the loop, we expect that it works both as a compressive one (defined in equation
(13): from neurons in the recurrent network over the output back to the neurons) and as
an expansive one (from the output through the neurons in the recurrent network back to the
PN
(p)
output). Rather than constraining the low-dimensional input β and output
to be
i dβi ai
P
equal (expansive auto-encoder), we can enforce the high dimensional input β eiβ β and output
PN
(p)
to be equal (compressive auto-encoder) for computing the decoders of the autoi,β ejβ dβi ai
encoder. Thus the squared-error loss becomes:
!2
P X X
N
X
X
(p)
(p)
L0 =
ejβ (
dβi ai − β )
p

=

j

P X
X
p

=

i

X

N
X

j

+

i

j

e2jβ

i

X

j

β

!

i

X

(p)

(p)

dβi ai − β )

ejγ (

N
X

γ

l

!
(p)

dγl al − (p)
γ )

!2
(p)

(p)

dβi ai − β

ejβ ejγ (

N
X

(p)
dβi ai

i

β,γ,β6=γ

P X X
N
X
p

N
X

β

P X
X
p

ejβ (

β

P XX
X
p

≈

β

−

N
X

(p)
β )(

l

!
(p)
dγl al

− (p)
γ )

!2
(p)

(p)

dβi ai − β

,

P P
where in the approximation, we consider that the term involving j β,γ,β6=γ ejβ ejγ tends to
P 2
zero as ejβ and ejγ are independent; and further j ejβ = 1. Thus, this loss function is approximately the same as the squared-error loss function in equation (6) used for the expansive
auto-encoder, showing that for an auto-encoder embedded in a loop, the expansive and compressive descriptions are equivalent. We employed a large number of random low-dimensional
inputs when constraining the expansive auto-encoder.
The FOLLOW learning rule i.e. equation (2) was applied on the feedforward and recurrent
ff
weights, namely wji
and wji . The error for our learning rule was the error β = xβ − x̂β in
the observable output ~x, not the error in the desired function f~(~x). The observable reference
state ~x was obtained by integrating the differential equations of the dynamical system. The
synaptic time constant τs was 20 ms in all synapses, including that for calculating the error and
for feeding the error back to the neurons. The error used for the weight update was filtered by
a 200 ms decaying exponential.
The command input vector ~u(t) to the network was Nc -dimensional (Nc = Nd for all our
simulations) and time-varying. During the learning phase, input changed on two different time
18

scales. The fast value of each command component was switched every 50 ms to a level u0α
chosen uniformly between (−ζ1 , ζ1 ) and this number was added to a more slowly changing input
variable ūα (called ’pedestal’ in the main part of the paper) which changed with a period Tperiod
as indicated in each of the figures. Here ūα is the component of a vector of length ζ2 with a
randomly chosen direction. The value of component α of the command is then uα = ūα + u0α .
Parameter values for the network and input for each dynamical system are provided in Table 1.
Further details are noted in subsection 4.3.
During the testing phase without error feedback, the network reproduced the reference trajectory of the dynamical system for a few seconds, in response to the same kind of input as during
learning. We also tested the network on a different input not used during learning as shown in
Figures 2 and 4.

4.3

Equations and parameters for the example dynamical systems

The equations and input modifications for each dynamical system are detailed below. Time
derivatives are in units of s−1 .
The equations for a linear decaying oscillators system (Supplementary Fig. S5) were
ẋ1 = u1 /0.02 + (−0.2x1 − x2 )/0.05
ẋ2 = u2 /0.02 + (x1 − 0.2x2 )/0.05.
For this linear dynamical system, we tested the learned network on a ramp for 2 s followed
by a step to a constant non-zero value. A ramp can be viewed as a preparatory input before
initiating an oscillatory movement, in a similar spirit to that observed in (pre-)motor cortex
[Churchland et al. 2012]. For such input too, the network tracked the reference for a few seconds
(Fig. S5A-C).
The equations for the van der Pol oscillator system were
ẋ1 = u1 /0.02 + x2 /0.125

ẋ2 = u2 /0.02 + 2(1 − x21 )x2 − x1 /0.125.
Each component of the pedestal input ūα was scaled differently for the van der Pol oscillator as
reported in Table 1.
The equations for the chaotic Lorenz system were
ẋ1 = u1 /0.02 + 10(x2 − x1 )
ẋ2 = u2 /0.02 − x1 x3 − x2
ẋ3 = u3 /0.02 + x1 x2 − 8(x3 + 28)/3.
In our equations above, x3 of the Lorenz equations was represented by an output variable
x˜3 = x3 − 28 so as to have observable variables that vary around zero. This does not change the
system dynamics, just its representation in the network. For the Lorenz system, only a pulse
at the start for 250 ms, chosen from a random direction of norm ζ1 , was provided to set off the
system, after which the system followed autonomous dynamics.

19

Number of neurons/layer
Tperiod (s)
Representation radius R1
Representation radius R2
Learning pulse ζ1
Learning pedestal ζ2

Linear
2000
2
0.2
1
R1 /6
R2 /16

van der Pol
3000
4
0.2
5
R1 /6, R1 /2
R1 /6, R1 /2

Lorenz
5000
20
6
30
R2 /10
0

Arm
Non-linear feedforward
5000
2000
2
2
0.2
1
1
R2 /0.3 R1 /0.6
R2 /0.3 R2 /1.6

Table 1: Network parameters for example systems

Our FOLLOW scheme also learned non-linear feedforward transforms, as demonstrated in
Supplementary Figures S6 and S7. For the non-linear feedforward case, we used the linear system
as the reference, but with the input transformed nonlinearly by gα (~u) = 10((uα /0.1)3 − uα /0.4).
Thus, the equations of the reference were:
ẋ1 = 10((u1 /0.1)3 − u1 /0.4) + (−0.2x1 − x2 )/0.05
ẋ2 = 10((u2 /0.1)3 − u2 /0.4) + (x1 − 0.2x2 )/0.05.
The input to the network remained ~u. Thus, the feedforward weights had to learn the non-linear
transform ~g (~u) while the recurrent weights learned the linear system.
In the example of learning arm dynamics, we used a two-link model for an arm moving in
the vertical plane with damping under gravity (see for example http://www.gribblelab.org/
compneuro/5_Computational_Motor_Control_Dynamics.html and https://github.com/studywolf/
control/tree/master/studywolf_control/arms/two_link), with parameters from [Li 2006].
The differential equations for the four state variables, namely the shoulder and elbow angles
θ~ = (θ1 , θ2 )T and the angular velocities ω
~ = (ω1 , ω2 )T , given input torques ~τ = (τ1 , τ2 )T were:

with

˙
θ~ = ω
~

(7)



−1
˙ω
~
~
~
~ = M (θ)
~τ − C(θ, ω
~ ) − B~ω − gD(θ)

(8)



2
2
2
d
+
2d
cos
θ
+
m
s
+
m
s
d
+
d
cos
θ
+
m
s
1
2
2
1
2
3
2
2
2
1
2
2
~ =
M (θ)
d3 + d2 cos θ2 + m2 s22
d3 + m2 s22




−θ̇2 (2θ̇1 + θ̇2 )
b11 b12
~
d2 sin θ2 , B =
C(θ, ω
~) =
,
b21 b22
θ̇12


(m1 s1 + m2 l1 ) sin θ1 + m2 s2 sin(θ1 + θ2 )
~
D(θ) =
,
m2 s2 sin(θ1 + θ2 )
d1 = I1 + I2 + m2 l12 , d2 = m2 l1 s2 , d3 = I2 ,

where mi is the mass, li the length, si the distance from the joint center to the center of the
mass, and Ii the moment of inertia, of link i; M is the moment of inertia matrix; C contains
centripetal and Coriolis terms; B is for joint damping; and D contains the gravitational terms.
Here, the state variable vector ~x = [θ1 , θ2 , ω1 , ω2 ], but the effective torque τ was obtained from
the input torque ~u as below.
To avoid any link from rotating full 360 degrees, we provided an effective torque τα to the
arm, by subtracting a term proportional to the input torque uα , if the angle crossed ±90 degrees

20

and uα was in the same direction:


uα > 0
uα σ̃(θα )
τ α = uα − 0
uα = 0 ,


uα σ̃(−θα ) uα < 0
where σ̃(θ) increases linearly from 0 to 1 as θ goes from π/2 to 3π/4 as below:


θ ≤ π/2
0
σ̃(θ) = (θ − π/2)/(π/4) 3π/4 > θ > π/2


1
θ ≥ 3π/4
The parameter values were as per Model 1 of the human arm in section 3.1.1 of the PhD thesis
of Li [Li 2006] from the Todorov lab; namely m1 = 1.4 kg, m2 = 1.1 kg, l1 = 0.3 m, l2 = 0.33 m,
s1 = 0.11 m, s2 = 0.16 m, I1 = 0.025 kg m2 , I2 = 0.045 kg m2 , and b11 = b22 = 0.05,
b12 = b21 = 0.025. Acceleration due to gravity was set at g = 9.81 m/s2 . For the arm, we did
not filter the reference variables for calculating the error.
The input torque ~u(t) for learning the two-link arm was generated, not by switching the
pulse and pedestal values sharply, every 50 ms and Tperiod as for the others, but by linearly
interpolating in-between to avoid oscillations from sharp transitions, due to the feedback loop
in the input in the general network (Supplementary Fig. S4).
The input torque ~u and the variables ω
~ , θ~ obtained on integrating the arm model above
were scaled by 0.02, 0.05 and 1./2.5 respectively, and then used as the reference for the spiking
network. Effectively, we scaled the input torques to cover one-fifth of the representation radius,
the angular velocities one-half, and the angles full, as each successive variable was the integral
of the previous one.

4.4

Derivation and proof of stability of the FOLLOW learning scheme

We derive the FOLLOW learning rules, while simultaneously proving the stability of the
scheme. We assume that: (1) the feedback {keiα } and readout weights {dαj } form an autoencoder with gain k; (2) given the gains and biases of the spiking LIF neurons, there exist
feedforward and recurrent weights that make the network follow the reference dynamics perfectly
(in practice, the dynamics is only approximately realizable by our network, see Supplementary
subsection 7.4 for a discussion); (3) the state ~x of the dynamical system is observable; (4) the
intrinsic time scales of the reference dynamics are much larger than the synaptic time scale, the
time scale of the error feedback loop, and the time scale of learning; (5) the feedforward and
recurrent weights remain bounded; and (6) the input ~u and reference output ~x remain bounded.
The proof proceeds in three major steps: using the auto-encoder assumption to write the
evolution equation of the low-dimensional output state variable in terms of the recurrent and
feedforward weights; showing that output follows the reference the error feedback loop; and
obtaining the evolution equation for the error and using it in the time-derivative of a Lyapunov
function V , to show that V̇ ≤ 0 for uniform stability, similar to proofs in adaptive control theory
[Narendra and Annaswamy 1989; Ioannou and Sun 2012].

21

Role of network weights for low-dimensional output
The filtered low-dimensional output of the recurrent network is given by
X
x̂α =
dαj (Sj ∗ κ)(t)

(9)

j

where dαj are the readout weights. Since κ is an exponential filter with time constant τs , equation
(9) can also be written as
X
τs x̂˙ α (t) = −x̂α (t) +
dαj Sj (t),
(10)
j

We convolve this equation with kernel κ, multiply by the error feedback weights, and sum over
the output components α
X
X
X
X
τs
eiα (x̂˙ α ∗ κ)(t) = −
eiα (x̂α ∗ κ)(t) +
eiα
dαj (Sj ∗ κ)(t).
(11)
α

α

α

j

We would like to write this evolution equation in terms of the recurrent and feedforward weight
in the networks.
To do this, we exploit assumptions (1) and (4). Having shown the equivalence of the compressive and expansive descriptions of our auto-encoder in the error-feedback loop (Methods
subsection 4.2), we formulate our nonlinear auto-encoder as compressive, which starts with a
high-dimensional set of inputs Ij ≡ (Jj − bj )/νj (where Jj is the current into neuron j having
gain νj and bias bj , cf. equations (3) and (4)); transforms these nonlinearly into filtered spike
trains Sj [Ij ] ∗ κ; decodes
these filtered spike trains into a low-dimensional representation ~z with
P
components zα = j dαj (Sj [Ij ] ∗ κ); and blows-up the dimensionality back to the original one,
via weights keiα , to get inputs:
XX
X
keiα zα = k
eiα dαj (Sj [Ij ] ∗ κ).
(12)
Ii0 =
α

α

j

The spike train Sj is a functional of the input Ij . Sj depends on neuron index j as neurons
have different gains νj and biases bj . Using assumption (1) we expect that the final inputs Ii0
are approximately k times the initial inputs Ii :
XX
k
eiα dαj (Sj [Ij ] ∗ κ) ≈ kIi .
(13)
α

j

Our assumption (4) says that the state variables of the reference dynamics change slowly
compared to neuronal dynamics. Due to the spatial averaging (sum over j in equation (13))
over a large number of neurons, individual neurons do not necessarily have to fire at a rate higher
than the inverse of the synaptic time scale, while we can still assume that the total round trip
input Ii0 on the left hand side of equation (13) is varying only on the slow time scale. Therefore,
we used firing rate equations to compute mean outputs given static input when pre-calculating
the readout weights (Methods subsection 4.2).
Inserting the approximate equation (13) in equation (11) we find
X
X
τs
eiα (x̂˙ α ∗ κ)(t) ≈ −
eiα (x̂α ∗ κ)(t) + Ii (t).
α

α

22

(14)

We replace Ii ≡ (Ji − bi )/νi , using the current Ji from equation (4) for neuron i of the recurrent
network, to obtain
X
X
X
τs
eiα (x̂˙ α ∗ κ)(t) ≈ −
eiα (x̂α ∗ κ)(t) +
wij (Sj ∗ κ)(t)
α

α

+

j

X
l

wilff (Slff ∗ κ)(t) +

X
α

keiα (α ∗ κ)(t).

(15)

Thus, the change of the low-dimensional output x̂α ∗ κ depends on the network weights, which
need to be learned. This finishes the first step of the proof.
Error-feedback loop ensures that output follows reference
Because of assumption (2), we may assume that there exists a recurrent network of spiking
neurons that represents the desired dynamics of equation (1) without any error feedback. This
second network serves as a target during learning and has variables and parameters indicated
with an asterisk. In particular, the second network has feedforward weights wilff* and recurrent
∗
weights wij
. We write an equation similar to equation (14) for the output x∗α of the target
network:
X
X
X
∗
τs
eiα (ẋ∗α ∗ κ)(t) = −
eiα (x∗α ∗ κ)(t) +
wij
(Sj∗ ∗ κ)(t)
α

α

j

+

X
l

wilff* (Slff* ∗ κ)(t),

(16)

where (Slff* ∗κ)(t) and (Sj∗ ∗κ)(t) are defined as the filtered spike trains of neurons in the realizable
target network. We emphasize that this target network does not need error feedback because
its output is, by definition, always correct. In fact, the readout from the spike trains Sj∗ gives
the target output which we denote by ~x∗ . The weights of the target network are constant and
their actual values are unimportant. They are mere mathematical devices to demonstrate stable
learning of the first network which has adaptable weights. For the first network, we choose the
same number of neurons and the same neuronal parameters as for the second network; moreover,
the weights from the command input to the representation layer and the readout weights from
the recurrent network to the output are identical for both networks. Thus, the only difference is
that the feedforward and recurrent weights of the target network are realized, while for the first
network they need to be learnt.
In view of potential generalization, we note that any nonlinear dynamical system is approximately realizable due to the expansion in a high-dimensional non-linear basis that is effectively
performed by the recurrent network (see Supplementary subsection 7.1). Approximative weights
(close to the ideal ones) could in principle also be calculated algorithmically as in Supplementary subsection 7.2. In the following we exploit assumption (2) and assume that the dynamics
is actually (and not only approximately) realized by the target network.
Our assumption (3) states that the output is observable. Therefore the error component α can
be computed directly via a comparison of the true output ~x of the reference with the output ~x̂ of
the network: α = xα − x̂α . (In view of potential generalizations, we remark that the observable
output need not be the state variables themselves, but could be a higher-dimensional nonlinear
function of the state variables, as shown for the general scheme in Supplementary section 7.3.)

23

As the second step of the proof, we now show that the error feedback loop enables the first
network
the target network under assumptions (4)-(6). More precisely, we want to show
P to follow P
that α eiα x̂α ≈ α eiα x∗α for each neuron index i. To do so, we use assumption (4) and exploit
that (i) learning is slow compared to the network dynamics so the weights of the first network
can be considered momentarily constant; (ii) the reference dynamics is slower than the synaptic
and feedback loop time scales, so the reference output xα can be assumed momentarily constant.
Thus, we have a separation of time scales in equation (15): for a given input (transmitted via the
feedforward weights) and a given target value x∗α , the network dynamics settles on the fast time
scale τs to a momentary fixed point x̂† which we find by setting the derivative on the left-hand
side of equation (15) to zero:
X
X
X
X
0=−
eiα (x̂†α ∗ κ)(t) +
wij (Sj ∗ κ)(t) +
wilff (Slff ∗ κ)(t) +
keiα ((x∗α − x̂†α ) ∗ κ)(t).
α

j

α

l

We rewrite this equation in the form
X
α

eiα (x̂†α ∗ κ)(t) =

k
k+1

X
α

eiα (x∗α ∗ κ)(t) +

1
k+1

!
X
j

wij (Sj ∗ κ)(t) +

X
l

wilff (Slff ∗ κ)(t) .

(17)
We choose the feedback gain for the error much larger than 1 (k  1), such that k/(k + 1) ≈ 1.
Using our assumption (5) that the feedforward and recurrent weights are bounded, and since the
filtered spike trains remain bounded due to refractory period, the term in parentheses multiplying
1/(k+1) remains bounded, by say B1 . Choosing k  B1 , the second term can be made negligible.
Thus, to obtain x̂†α ≈ x∗α , we set k  1 and k  B1 .
To show that the fixed point is stable at the fast synaptic time scale, we calculate the Jacobian
J =P
[Jil ], for the dynamical system given by equation (15). We introduce auxiliary variables
yi ≡ α eiα x̂α to rewrite equation (15) with the new variables in the form ẏi = Fi (~y ); and then
we take derivative
of its right hand side to obtain the elements of the Jacobian matrix at the
P
fixed point α eiα x̂†α :
P
Z t
∂ j wij (Sj ∗ κ)(t)
∂Fi (~y )
= −(k + 1)δil
κ(τ )dτ +
,
Jil ≡
P
∂yl
∂yl
−∞
yi = α eiα x̂†α
P
where δil is the Kronecker delta function. We note that j wij (Sj ∗ κ) is a spatially and temporally averaged measure of the population activity in the network with appropriate weighting
factors wij . We assume that the population activity varies smoothly with input, which is equivalent to requiring that on the time scale τs , the network fires asynchronously, i.e. there are no
precisely timed population spikes. Then we can take the second term to be bounded, by say B2 .
If we make k  B2 , then the Jacobian matrix J has large negative values on the diagonal and
negligible values off-diagonal. Effectively, the Jacobian has negative eigenvalues, rendering the
momentary fixed point asymptotically stable.
Thus, we have shown that if the initial state of the first network is close to the initial state of the
target network, e.g. both start from rest, then on the slow time scale of the
of
P system dynamics
P
the reference ~x∗ , the first network follows the target network at all times, α eiα x̂α ≈ α eiα x∗α .
With these constraints on each neuron, and since the readout weights, error-encoding weights
and neuronal parameters are the same for the first and second network, the actual filtered spike
trains of the recurrent neurons in the first network will be approximately the same as those of
the target network, so that (Si ∗ κ)(t) can be used instead of (Si∗ ∗ κ)(t) in (16). Moreover,
the filtered spike trains (Slff ∗ κ)(t) of the command representation layer in the first network are
always the same as those in the target network, since they are driven by the same command
input ~u and the command encoding weights are, by construction, the same for both networks.
24

Stability of learning via Lyapunov’s method
We now turn to the third step of the proof and consider the temporal evolution of the error
α = xα − x̂α . We exploit that the network dynamics is realized by the target network and insert
equations (15) and (16) so as to find
X
X
−τs
eiα (˙α ∗ κ)(t) = τs
eiα ((x̂˙ α − ẋα ) ∗ κ)(t)
α

α

≈ τs
≈

X

X
j

α

eiα ((x̂˙ α − ẋ∗α ) ∗ κ)(t)

X


∗
wij − wij
(Sj ∗ κ)(t) +
wilff − wilff* (Slff ∗ κ)(t)

+ (k + 1)

X
α

≡

X
j

(18)

l

eiα (α ∗ κ)(t)

ψij (Sj ∗ κ)(t) +

X
l

φil (Slff ∗ κ)(t) + (k + 1)

X
α

eiα (α ∗ κ)(t),

In the second line, we have replaced the reference output by the target network output; and in
the third line we have used equations (15) and (16), and replaced the filtered spike trains of the
target network by those of the first network, exploiting the insights from the previous paragraph.
∗
In the last line, we have introduced abbreviations ψij ≡ wij − wij
and φil ≡ wilff − wilff* .
In order to show that the absolute value of the error decreases over time with an appropriate
learning rule, we consider the candidate Lyapunov function:
V (˜, ψ, φ) =

11 X
1X 2 1 1 X
˜i +
(ψij )2 +
(φil )2 ,
2 i
2 η̃1 i,j
2 η̃2 i,l

(19)

P
where ˜i ≡ τs α eiα (α ∗ κ) and η̃1 , η̃2 > 0 are positive constants. The Lyapunov function is
positive semi-definite V (˜, ψ, φ) ≥ 0, with the equality to zero only at (˜, ψ, φ) = (0, 0, 0). It has
continuous first-order partial derivatives.
Furthermore, V is radially unbounded since
V (˜, ψ, φ) > |(˜, ψ, φ)|2 /(4 max(1, η̃1 , η̃2 )),
and decrescent since
V (˜, ψ, φ) < |(˜, ψ, φ)|2 / min(1, η̃1 , η̃2 ),
P
P
P
where |(˜, ψ, φ)|2 ≡ i (˜i )2 + i,j (ψij )2 + i,k (φil )2 and min / max take the minimum / maximum
of their respective arguments.
Using Lyapunov’s direct method, we need to prove, apart from the above conditions, the
property V̇ ≤ 0 for uniform global stability (which implies that bounded orbits remain bounded,
so the error remains bounded); or the stronger property V̇ < 0 for asymptotic global stability
(see for example [Narendra and Annaswamy 1989; Ioannou and Sun 2012]). Taking the time

25

derivative of V , and replacing ˜˙i i.e. τs
V̇ =

X
i

P

α eiα (˙α

∗ κ) from (18), we have:

1 X
1 X
ψij ψ̇ij +
˜i ˜˙i +
φil φ̇il
η̃1 i,j
η̃2 i,l
!

≈−

X
i

˜i

X
j

ψij (Sj ∗ κ)(t) +

X
l

φil (Slff ∗ κ)(t) + (k + 1)

1 X
1 X
ψij ψ̇ij +
φil φ̇il
+
η̃1 i,j
η̃2 i,l


X
1
=
ψij −˜i (Sj ∗ κ)(t) + ψ̇ij
η̃1
i,j


X
X
1
ff
+
φil −˜i (Sl ∗ κ)(t) + φ̇il − (k + 1)
˜2i /τs .
η̃2
i
i,k

X
α

eiα (α ∗ κ)(t)
(20)

If we pick as a learning rule
ψ̇ij = η̃1 ˜i (Sj ∗ κ)(t)

φ̇il = η̃2 ˜i (Slff ∗ κ)(t),

(21)

then
V̇ = −(k + 1)

X
i

˜2i /τs ≤ 0

choosing k > −1, which is subsumed under k  1 for the error feedback. The condition (21)
with η1 ≡ η̃1 τs and η2 ≡ η̃2 τs , and κ replaced by a longer filtering kernel κ , is the learning rule
used in the main text, equation (2).
Thus, in the (˜, ψ, φ)-system given by equations (18) and (21), we have proven the global
uniform stability of the fixed point (˜, ψ, φ) = (0, 0, 0), which is effectively (, ψ, φ) = (0, 0, 0),
choosing η1 , η2 > 0 and k  max(1, B1 , B2 ), under assumptions (1)-(6).
This ends our proof. So far, we have shown that the system is Lyapunov stable i.e. bounded
orbits remain bounded, and not asymptotically stable. Indeed, with bounded firing rates and
fixed readout weights, the output will remain bounded, as will the error (for a bounded reference).
However, here, we also derived the FOLLOW learning rule, and armed with the inequality for the
time derivative of the Lyapunov function in terms of the error, we further show in the following
Methods subsection 4.5 that the error ~ goes to zero asymptotically, so that even without error
feedback, ~x̂ reproduces the dynamics of ~x after learning.
A major caveat of this proof is that under assumption (2) the dynamics are realizable by our
network. In a real application this might not be the case. Approximation errors arising from
a mismatch between the best possible network and the actual target dynamics are currently
ignored. The adaptive control literatue has shown that errors in approximating the reference
dynamics appear as frozen noise and can cause runaway drift of the parameters [Narendra and
Annaswamy 1989; Ioannou and Sun 2012]. In our simulations with a large number of neurons,
the approximations of a non-realizable reference dynamics (e.g., the Van der Pol oscillator) were
sufficiently good, and thus the expected drift was possibly slow, and did not cause the error
to rise during typical timescales of learning. A second caveat is our assumption (5). While
the input is under our control and can therefore be kept bounded, some additional bounding is
26

needed to stop weights from drifting. Various techniques to address such model-approximation
noise and bounding weights have been studied in the robust adaptive control literature (e.g.
[Ioannou and Tsakalis 1986; Slotine and Coetsee 1986; Narendra and Annaswamy 1989; Ioannou
and Fidan 2006; Ioannou and Sun 2012]). We discuss this issue and briefly mention some of
these ameliorative techniques in Supplementary sub-section 7.4.
To summarize, the FOLLOW learning rule (21) on the feedforward or recurrent weights has
two terms: (i) a filtered presynaptic firing trace (Slff ∗P
κ)(t) or (Sj ∗ κ)(t) that is available locally
at each synapse; and (ii) a projected filtered error α eiα (α ∗ κ)(t) used for all synapses in
neuron i that is available as a current in the postsynaptic neuron i due to error feedback, see
equation (4). Thus the learning rule can be classified as local. Moreover, it uses an error in the
observable ~x, not in its time-derivative. While we have focused on spiking networks, the learning
scheme can be easily used for non-linear rate units by replacing the filtered spikes (Si ∗ κ)(t)
by the output of the rate units r(t). The rate units must also low-pass filter their input with
kernel κ(t). Our proof is valid for arbitrary dynamical transforms f~(~x) +~g (~u) as long as they are
realizable in a network. The proof does not involve gradient descent but shows uniform global
stability using Lyapunov’s method.

4.5

Proof of error tending to zero asymptotically

P
2
In section 4.4 we showed uniform
P global stability using V̇ = −(k + 1) i (˜i ) ≤ 0, with
k  max(1, B1 , B2 ) and ˜i ≡ τs α ejα (α ∗ κ). This only means that bounded errors remain
bounded. Here, we show more importantly that the error tends to zero asymptotically with
time. We adapt the proof in section 4.2 of [Ioannou and Sun 2012], to our spiking network.
Here, we want to invoke a special case of Barbălat’s lemma: if f, f˙ ∈ L∞ and f ∈ Lp for some
p ∈ [1, ∞), then f (t) → 0 as t → ∞. Recall the definitions: function f ∈ Lp when ||x||p ≡
1/p
R∞
|f (τ )|p dτ
exists (is finite); and similarly function f ∈ L∞ when ||x||∞ ≡ supt≥0 |f (τ )|
0
exists (is finite).
Since V is positive semi-definite (V ≥ 0) and is a non-increasing function of time (V̇ ≤ 0), its
limt→∞ V = V∞ exists and is finite. Using this, the following limit exists and is finite:
Z ∞
XZ ∞
1
−1
2
V̇ (τ )dτ =
(V (0) − V∞ ).
(˜i (τ )) dτ =
k
+
1
k
+
1
0
0
i
Since each term in the above sum
is finite ∀i, and thus ˜i ∈ L2 ∀i.

P

i

is positive semi-definite,

R∞
0

(˜i (τ ))2 dτ also exists and

To show that ˜i , ˜˙i ∈ L∞ ∀i, consider equation (18). First, using the assumptions that the
input ~u(t) is bounded and the reference dynamics is stable, we have that the reference output
~x(t) is bounded. Since network output ~x̂ is also bounded due to saturation of firing rates (as
are the filtered spike trains), the error (each component) is bounded i.e. ˜i ∈ L∞ ∀i. If we
also bound the weights from diverging during learning, then ψij , φil ∈ L∞ ∀i, j, k. With these
reasonable assumptions, all terms on the right hand side of the equation (18) for ˜˙i are bounded,
hence ˜˙i ∈ L∞ ∀i.

27

Since ˜i ∈ L2 ∀i and ˜i , ˜˙i ∈ L∞ ∀i, invoking Barbălat’s lemma as above, we have ˜i → 0 ∀i as
t → ∞. We have shown that the error tends to zero asymptotically under assumptions (1)-(6).
In practice, the error shows fluctuations on a short time scale while the mean error over a longer
time scale reduces and then plateaus, possibly due to approximate realizability, imperfections in
the error-feedback, and spiking shot noise (cf. Fig. 5).
We do not further require the convergence of parameters to ideal ones for our purpose, since
the error tending to zero, i.e. network output matching reference, is functionally sufficient for the
forward predictive model. In the adaptive control literature [Ioannou and Sun 2012; Narendra
and Annaswamy 1989], the parameters are shown to converge to ideal ones if input excitation is
“persistent”, loosely that it excites all modes of the system. It should be possible to adapt the
proof to our spiking network, as suggested by simulations (cf. Fig. 5), but is not pursued here.

5

Acknowledgements

We thank Johanni Brea, Samuel Muscinelli and Laureline Logiaco for helpful discussions and
comments on a previous version of the manuscript. We thank Chris Stock, Tilo Schwalger, Olivia
Gozel, and Dane Corneil for comments on the manuscript. Financial support was provided by
the European Research Council (Multirules, grant agreement no. 268 689) and by the Swiss
National Science Foundation (grant agreement no. CRSII 147636).

6

References

Abbott, L. F., DePasquale, B., and Memmesheimer, R.-M. (2016). Building functional networks
of spiking model neurons. Nature Neuroscience, 19(3):350–355.
Bengio, Y., Simard, P., and Frasconi, P. (1994). Learning long-term dependencies with gradient
descent is difficult. IEEE Transactions on Neural Networks, 5(2):157–166.
Bourdoukan, R. and Denève, S. (2015). Enforcing balance allows local supervised learning in
spiking recurrent networks. In Cortes, C., Lawrence, N. D., Lee, D. D., Sugiyama, M., Garnett,
R., and Garnett, R., editors, Advances in Neural Information Processing Systems 28, pages
982–990. Curran Associates, Inc.
Brown, S. P. and Hestrin, S. (2009). Intracortical circuits of pyramidal neurons reflect their
long-range axonal targets. Nature, 457(7233):1133–1136.
Burbank, K. S. (2015). Mirrored STDP Implements Autoencoder Learning in a Network of
Spiking Neurons. PLoS Computational Biology, 11(12).
Burnod, Y., Grandguillaume, P., Otto, I., Ferraina, S., Johnson, P. B., and Caminiti, R. (1992).
Visuomotor transformations underlying arm movements toward visual targets: A neural network model of cerebral cortical operations. Journal of Neuroscience, 12(4):1435–1453.
Churchland, M. M., Cunningham, J. P., Kaufman, M. T., Foster, J. D., Nuyujukian, P., Ryu,
S. I., and Shenoy, K. V. (2012). Neural population dynamics during reaching. Nature,
487(7405):51–56.
Conant, R. C. and Ashby, W. R. (1970). Every good regulator of a system must be a model of
that system. Intl. J. Systems Science, pages 89–97.

28

Dadarlat, M. C., O’Doherty, J. E., and Sabes, P. N. (2015). A learning-based approach to
artificial sensory feedback leads to optimal integration. Nature Neuroscience, 18(1):138–144.
Davidson, P. R. and Wolpert, D. M. (2005). Widespread access to predictive models in the
motor system: A short review. Journal of Neural Engineering, 2(3):S313.
Denève, S. and Machens, C. K. (2016). Efficient codes and balanced networks. Nature Neuroscience, 19(3):375–382.
DePasquale, B., Churchland, M. M., and Abbott, L. F. (2016). Using Firing-Rate Dynamics to
Train Recurrent Networks of Spiking Model Neurons. arXiv:1601.07620 [q-bio].
DeWolf, T., Stewart, T. C., Slotine, J.-J., and Eliasmith, C. (2016). A spiking neural model of
adaptive arm control. Proc. R. Soc. B, 283(1843):20162134.
Eliasmith, C. (2005). A Unified Approach to Building and Controlling Spiking Attractor Networks. Neural Computation, 17(6):1276–1314.
Eliasmith, C. and Anderson, C. H. (2004). Neural Engineering: Computation, Representation,
and Dynamics in Neurobiological Systems. MIT Press.
Florian, R. V. (2012). The Chronotron: A Neuron That Learns to Fire Temporally Precise Spike
Patterns. PLOS ONE, 7(8):e40233.
Friston, K. (2008).
4(11):e1000211.

Hierarchical Models in the Brain.

PLOS Computational Biology,

Funahashi, K.-I. (1989). On the approximate realization of continuous mappings by neural
networks. Neural Networks, 2(3):183–192.
Gardner, B. and Grüning, A. (2016). Supervised Learning in Spiking Neural Networks for Precise
Temporal Encoding. PLOS ONE, 11(8):e0161335.
Gerstner, W., Kistler, W. M., Naud, R., and Paninski, L. (2014). Neuronal Dynamics: From
Single Neurons to Networks and Models of Cognition. Cambridge University Press, Cambridge,
United Kingdom, 1st edition.
Girosi, F. and Poggio, T. (1990). Networks and the best approximation property. Biological
Cybernetics, 63(3):169–176.
Gütig, R. (2014). To spike, or when to spike? Current opinion in neurobiology, 25C:134–139.
Gütig, R. and Sompolinsky, H. (2006). The tempotron: A neuron that learns spike timing–based
decisions. Nature Neuroscience, 9(3):420–428.
Hennequin, G., Vogels, T. P., and Gerstner, W. (2014). Optimal Control of Transient Dynamics
in Balanced Networks Supports Generation of Complex Movements. Neuron, 82(6):1394–1406.
Hilber, P. and Caston, J. (2001). Motor skills and motor learning in Lurcher mutant mice during
aging. Neuroscience, 102(3):615–623.
Hochreiter, S., Bengio, Y., Frasconi, P., and Schmidhuber, J. (2001). Gradient Flow in Recurrent
Nets: The Difficulty of Learning Long-Term Dependencies.
Hochreiter, S. and Schmidhuber, J. (1997). Long Short-Term Memory. Neural Computation,
9(8):1735–1780.
29

Hoerzer, G. M., Legenstein, R., and Maass, W. (2014). Emergence of Complex Computational
Structures From Chaotic Neural Networks Through Reward-Modulated Hebbian Learning.
Cerebral Cortex, 24(3):677–690.
Hornik, K., Stinchcombe, M., and White, H. (1989). Multilayer feedforward networks are universal approximators. Neural Networks, 2(5):359–366.
Ioannou, P. and Fidan, B. (2006). Adaptive Control Tutorial. SIAM, Society for Industrial and
Applied Mathematics, Philadelphia, PA.
Ioannou, P. and Sun, J. (2012). Robust Adaptive Control. Dover Publications, Mineola, New
York, first edition.
Ioannou, P. and Tsakalis, K. (1986). A robust direct adaptive controller. IEEE Transactions on
Automatic Control, 31(11):1033–1043.
Jaeger, H. (2001). The ”echo state” approach to analysing and training recurrent neural networks. Technical report.
Jaeger, H. (2005). A Tutorial on Training Recurrent Neural Networks, Covering BPPT, RTRL,
EKF and the ”Echo State Network” Approach.
Jaeger, H. and Haas, H. (2004). Harnessing Nonlinearity: Predicting Chaotic Systems and
Saving Energy in Wireless Communication. Science, 304(5667):78–80.
Joshi, P. and Maass, W. (2005). Movement Generation with Circuits of Spiking Neurons. Neural
Computation, 17(8):1715–1738.
Khazipov, R., Sirota, A., Leinekugel, X., Holmes, G. L., Ben-Ari, Y., and Buzsáki, G. (2004).
Early motor activity drives spindle bursts in the developing somatosensory cortex. Nature,
432(7018):758–761.
Lalazar, H. and Vaadia, E. (2008). Neural basis of sensorimotor learning: Modifying internal
models. Current Opinion in Neurobiology, 18(6):573–581.
Legenstein, R., Chase, S. M., Schwartz, A. B., and Maass, W. (2010). A Reward-Modulated
Hebbian Learning Rule Can Explain Experimentally Observed Network Reorganization in a
Brain Control Task. Journal of Neuroscience, 30(25):8400–8410.
Legenstein, R. and Maass, W. (2007). Edge of chaos and prediction of computational performance for neural circuit models. Neural Networks, 20(3):323–334.
Legenstein, R., Markram, H., and Maass, W. (2003). Input prediction and autonomous movement analysis in recurrent circuits of spiking neurons. Reviews in the Neurosciences, 14(12):5–19.
Li, W. (2006). Optimal Control for Biological Movement Systems. PhD thesis, University of
California, San Diego.
Lillicrap, T. P., Cownden, D., Tweed, D. B., and Akerman, C. J. (2016). Random synaptic
feedback weights support error backpropagation for deep learning. Nature Communications,
7:13276.
Lorenz, E. N. (1963). Deterministic Nonperiodic Flow. Journal of the Atmospheric Sciences,
20(2):130–141.

30

Maass, W. and Markram, H. (2004). On the computational power of circuits of spiking neurons.
Journal of Computer and System Sciences, 69(4):593–616.
Maass, W., Natschläger, T., and Markram, H. (2002). Real-Time Computing Without Stable
States: A New Framework for Neural Computation Based on Perturbations. Neural Computation, 14(11):2531–2560.
MacNeil, D. and Eliasmith, C. (2011). Fine-Tuning and the Stability of Recurrent Neural
Networks. PLoS ONE, 6(9):e22885.
Markram, H., Muller, E., Ramaswamy, S., Reimann, M. W., Abdellah, M., Sanchez, C. A.,
Ailamaki, A., Alonso-Nanclares, L., Antille, N., Arsever, S., Kahou, G. A. A., Berger, T. K.,
Bilgili, A., Buncic, N., Chalimourda, A., Chindemi, G., Courcol, J.-D., Delalondre, F., Delattre, V., Druckmann, S., Dumusc, R., Dynes, J., Eilemann, S., Gal, E., Gevaert, M. E.,
Ghobril, J.-P., Gidon, A., Graham, J. W., Gupta, A., Haenel, V., Hay, E., Heinis, T., Hernando, J. B., Hines, M., Kanari, L., Keller, D., Kenyon, J., Khazen, G., Kim, Y., King, J. G.,
Kisvarday, Z., Kumbhar, P., Lasserre, S., Le Bé, J.-V., Magalhães, B. R. C., Merchán-Pérez,
A., Meystre, J., Morrice, B. R., Muller, J., Muñoz-Céspedes, A., Muralidhar, S., Muthurasa,
K., Nachbaur, D., Newton, T. H., Nolte, M., Ovcharenko, A., Palacios, J., Pastor, L., Perin,
R., Ranjan, R., Riachi, I., Rodrı́guez, J.-R., Riquelme, J. L., Rössert, C., Sfyrakis, K., Shi, Y.,
Shillcock, J. C., Silberberg, G., Silva, R., Tauheed, F., Telefont, M., Toledo-Rodriguez, M.,
Tränkler, T., Van Geit, W., Dı́az, J. V., Walker, R., Wang, Y., Zaninetta, S. M., DeFelipe, J.,
Hill, S. L., Segev, I., and Schürmann, F. (2015). Reconstruction and Simulation of Neocortical
Microcircuitry. Cell, 163(2):456–492.
Memmesheimer, R.-M., Rubin, R., Ölveczky, B. P., and Sompolinsky, H. (2014). Learning
Precisely Timed Spikes. Neuron.
Mohemmed, A., Schliebs, S., Matsuda, S., and Kasabov, N. (2012). Span: Spike pattern association neuron for learning spatio-temporal spike patterns. International Journal of Neural
Systems, 22(4):1250012.
Morse, A. (1980). Global stability of parameter-adaptive control systems. IEEE Transactions
on Automatic Control, 25(3):433–439.
Narendra, K., Lin, Y.-H., and Valavani, L. (1980). Stable adaptive controller design, part II:
Proof of stability. IEEE Transactions on Automatic Control, 25(3):440–448.
Narendra, K. S. and Annaswamy, A. M. (1989). Stable Adaptive Systems. Prentice-Hall, Inc.
Nicola, W. and Clopath, C. (2016). Supervised Learning in Spiking Neural Networks with
FORCE Training. arXiv:1609.02545 [q-bio].
Parisien, C., Anderson, C. H., and Eliasmith, C. (2008). Solving the Problem of Negative
Synaptic Weights in Cortical Models. Neural Computation, 20(6):1473–1494.
Pearlmutter, B. A. (1995). Gradient calculations for dynamic recurrent neural networks: A
survey. IEEE Transactions on Neural Networks, 6(5):1212–1228.
Petersson, P., Waldenström, A., Fåhraeus, C., and Schouenborg, J. (2003). Spontaneous muscle
twitches during sleep guide spinal self-organization. Nature, 424(6944):72–75.
Pfister, J.-P., Toyoizumi, T., Barber, D., and Gerstner, W. (2006). Optimal spike-timingdependent plasticity for precise action potential firing in supervised learning. Neural Computation, 18(6):1318–1348.
31

Poggio, T. (1990). A theory of how the brain might work. Cold Spring Harbor Symposia on
Quantitative Biology, 55:899–910.
Pouget, A. and Sejnowski, T. J. (1997). Spatial Transformations in the Parietal Cortex Using
Basis Functions. Journal of Cognitive Neuroscience, 9(2):222–237.
Pouget, A. and Snyder, L. H. (2000). Computational approaches to sensorimotor transformations. Nature Neuroscience, 3:1192–1198.
Roelfsema, P. R. and van Ooyen, A. (2005). Attention-gated reinforcement learning of internal
representations for classification. Neural Computation, 17(10):2176–2214.
Rosenblatt, F. (1961). Principles of neurodynamics. perceptrons and the theory of brain mechanisms. Technical report.
Rumelhart, D., Hinton, G., and Williams, R. (1986). Learning Internal Representations by
Error Propagation. In Parallel Distributed Processing: Explorations in the Microstructure of
Cognition, Vol 1, volume 1. MIT Press Cambridge, MA, USA.
Sanner, R. M. and Slotine, J. J. E. (1992). Gaussian networks for direct adaptive control. IEEE
Transactions on Neural Networks, 3(6):837–863.
Sarlegna, F. R. and Sainburg, R. L. (2009). The Roles of Vision and Proprioception in the
Planning of Reaching Movements. Advances in experimental medicine and biology, 629:317–
335.
Sastry, S. and Bodson, M. (1989). Adaptive Control: Stability, Convergence, and Robustness.
Prentice-Hall, Inc., Upper Saddle River, NJ, USA.
Seung, H. S., Lee, D. D., Reis, B. Y., and Tank, D. W. (2000). Stability of the Memory of Eye
Position in a Recurrent Network of Conductance-Based Model Neurons. Neuron, 26(1):259–
271.
Slotine, J.-J. E. and Coetsee, J. A. (1986). Adaptive sliding controller synthesis for non-linear
systems. International Journal of Control, 43(6):1631–1651.
Slotine, J.-J. E. and Weiping Li (1987). On the Adaptive Control of Robot Manipulators. The
International Journal of Robotics Research, 6(3):49–59.
Song, H. F., Yang, G. R., and Wang, X.-J. (2016). Training Excitatory-Inhibitory Recurrent
Neural Networks for Cognitive Tasks: A Simple and Flexible Framework. PLOS Comput Biol,
12(2):e1004792.
Stewart, T. C., Tripp, B., and Eliasmith, C. (2009). Python scripting in the Nengo simulator.
Frontiers in Neuroinformatics, 3.
Sussillo, D. and Abbott, L. (2012). Transferring Learning from External to Internal Weights in
Echo-State Networks with Sparse Connectivity. PLoS ONE, 7(5):e37372.
Sussillo, D. and Abbott, L. F. (2009). Generating Coherent Patterns of Activity from Chaotic
Neural Networks. Neuron, 63(4):544–557.
Sutton, R. S. (1996). Generalization in reinforcement learning: Successful examples using sparse
coarse coding. Advances in Neural Information Processing Systems 8, pages 138–1044.
Thalmeier, D., Uhlmann, M., Kappen, H. J., and Memmesheimer, R.-M. (2016). Learning
Universal Computations with Spikes. PLOS Comput Biol, 12(6):e1004895.
32

Urbanczik, R. and Senn, W. (2014). Learning by the Dendritic Prediction of Somatic Spiking.
Neuron, 81(3):521–528.
Voegtlin, T. (2006). Temporal Coding using the Response Properties of Spiking Neurons. pages
1457–1464.
Vogels, T. P., Sprekeler, H., Zenke, F., Clopath, C., and Gerstner, W. (2011). Inhibitory
Plasticity Balances Excitation and Inhibition in Sensory Pathways and Memory Networks.
Science, 334(6062):1569 –1573.
Williams, R. J. and Zipser, D. (1989). A Learning Algorithm for Continually Running Fully
Recurrent Neural Networks. Neural Computation, 1(2):270–280.
Wolpert, D. M. and Ghahramani, Z. (2000). Computational principles of movement neuroscience.
Nature Neuroscience, 3:1212–1217.
Wong, J. D., Kistemaker, D. A., Chin, A., and Gribble, P. L. (2012). Can proprioceptive training
improve motor learning? Journal of Neurophysiology, 108(12):3313–3321.
Zago, M., Bosco, G., Maffei, V., Iosa, M., Ivanenko, Y. P., and Lacquaniti, F. (2005). Fast
Adaptation of the Internal Model of Gravity for Manual Interceptions: Evidence for EventDependent Learning. Journal of Neurophysiology, 93(2):1055–1068.
Zago, M., McIntyre, J., Senot, P., and Lacquaniti, F. (2009). Visuo-motor coordination and
internal models for object interception. Experimental Brain Research, 192(4):571–604.
Zerkaoui, S., Druaux, F., Leclercq, E., and Lefebvre, D. (2009). Stable adaptive control with
recurrent neural networks for square MIMO non-linear systems. Engineering Applications of
Artificial Intelligence, 22(4–5):702–717.

33

7

Supplementary Material

7.1

Decoding

Consider only the command representation layer without the subsequent recurrent network.
Assume, following [Eliasmith and Anderson 2004], we wish to decode an arbitrary output ~v (~u)
corresponding to the ~u encoded in the command representation layer, from the spike trains Slff (t)
of the neurons, by synaptically filtering and linearly weighting the trains with decoding weights
(~v )
dαl :
X (~v)
v̂α (~u) =
dαl (Slff ∗ κ)(t),
(22)
l

R∞
Rt
where ∗ denotes convolution (Slff ∗ κ)(t) ≡ −∞ Slff (t0 )κ(t − t0 )dt0 = 0 Slff (t − t0 )κ(t0 )dt0 , and
κ(t) ≡ exp(−t/τs )/τs is a normalized filtering kernel.
(~v )

We can obtain the decoders dαi by minimizing the loss function
*
!2 +
X
X (~v)
L=
vα (~u) −
dαl hSlff ∗ κit
α

l

(23)

~
u

with respect to the decoders. The average h·i~u over ~u guarantees that the same constant decoders
are used over the whole range of constant inputs ~u. The time average h·it denotes an analytic rate
computed for each constant input for a LIF neuron. Linear regression with a finite set of constant
inputs ~u was used to obtain the decoders (see Methods). With these decoders, if the input ~u
P (~v)
varies slowly compared to the synaptic time constant τs , we have v̂α = l dαl (Slff ∗κ)(t) ≈ vα (~u).
Any function of the input ~v (~u) can be approximated with appropriate linear decoding weights
from the high-dimensional basis of non-linear tuning curves of heterogeneous neurons with
different biases, encoding weights and gains, schematized in Figure 1B. With a large enough
number of such neurons, the function is expected to be approximated to arbitrary accuracy.
While this has not been proven rigorously for spiking neurons, this has theoretical underpinnings
from theorems on universal function approximation using non-linear basis functions [Funahashi
1989; Hornik et al. 1989; Girosi and Poggio 1990] sucessful usage in spiking neural network
models by various groups [Seung et al. 2000; Eliasmith and Anderson 2004; Eliasmith 2005;
Pouget and Sejnowski 1997], and biological plausibility [Poggio 1990; Burnod et al. 1992; Pouget
and Sejnowski 1997].
(~v )
dαl

Here, the neurons that are active at any given time operate in the mean driven regime i.e. the
instantaneous firing rate increases with the input current [Gerstner et al. 2014]. The dynamics
is dominated by synaptic filtering, and the membrane time constant does not play a significant
role [Eliasmith and Anderson 2004; Eliasmith 2005; Seung et al. 2000; Abbott et al. 2016]. Thus,
the decoding weights derived from equation (23) with stationary input are good approximations
even in the time-dependent case, as long as the input varies on a time scale slower than the
synaptic time constant.

7.2

Online learning based on a loss function and its shortcomings

The dynamical system given by equation (1) is to be mimicked by our spiking network implementing a different dynamical system with an extra error feedback term as in equation (15).
This can be interpreted as:
τs x̂˙ α = −x̂α + f˘α (~x̂) + ğα (~u) + kα .
34

(24)

Comparing with the reference equation (1), after learning we want that f˘α (~x̂) and ğα (~u) should
approximate f˜α (~x̂) ≡ τs fα (~x̂) + x̂α and g̃α (~u) ≡ τs gα (~u) respectively (upto a constant term
distributed between the two functions). In our simulations, we will usually start with zero
feedforward and recurrent weights, so that initially f˘(~x̂) = 0 = ğα (~u).]
Assuming that the time scales of dynamics are slower than synaptic time scale τs , we can
approximate the requisite feedforward and recurrent weights, by minimizing the following loss
functions respectively, with respect to the weights [Eliasmith and Anderson 2004]:
*
!2 +
X X
X
effkα g̃α (~u) −
wjl hSlff ∗ κit
,
(25)
Lff =
α

j

Lrec =

*
X X
j

α

l

x

!2 +
ejα f˜α (~x) −

X
i

wji hSi ∗ κit

.

(26)

x

Using these loss functions, we can precalculate the weights required for any dynamical system
numerically, similarly to the calculation of decoders in Supplementary subsection 7.1.
We now derive rules for learning the weights online based on stochastic gradient descent of
these loss functions, similar to [MacNeil and Eliasmith 2011], and point out some shortcomings.
The learning rule for the recurrent weights by gradient descent on the loss function given by
equation (26) is
1 ∂Lrec
dwji
=− η
dt
2 ∂wji
*
!
+
X
X
≈η
ejβ f˜β (~x) −
wji (Si ∗ κ)(t) (Si ∗ κ)(t)
i

β

≡η

D

(f˜)
j (Si

∗ κ)(t)

E

(27)

x

.
x

In the second line, the effect of the weight change on the filtered spike trains is assumed small and
neglected, using a small learning rate η. With requisite dynamics
slower than P
synaptic τs , and
P
with large enough number of neurons, we have approximated i wji hSi ∗κit (t) ≈ i wji (Si ∗κ)(t).
The third line defines an error in the projected f~˜(~x), which is the supervisory signal.
If we assume that the learning rate is slow, and the input samples the range of x uniformly,
then we can remove the averaging over x, similar to stochastic gradient descent.

(f˜)

where j

dwji
(f˜)
≈ −ηj (Si ∗ κ)(t).
(28)
dt
P

P
˜
≡
x) − i wji (Si ∗ κ)(t) . This learning rule is the product of a projected
β ejβ fβ (~
(f˜)

and the filtered presynaptic spike train (Si ∗ κ)(t). However, this
projected error in the unobservable f~˜ is not available to the postsynaptic neuron, making the
learning rule non-local. A similar issue arises in the feedforward case.
multi-dimensional error j

35

In mimicking a dynamical system, we want only the observable output of the dynamical system
i.e. ~x to be used in a supervisory signal, not a term involving the unknown f (~x) appearing in
the derivative ~x˙ . Even if this derivative is computed from the observable ~x, it will be noisy.
Furthermore, this derivative cannot be obtained by differentiating the observable versus time,
if the observable is not directly the output, but an unknown non-linear function of it, which
however our FOLLOW learning can handle (see Supplementary subsection 7.3). Thus, an online
rule using just the observable error can learn only an integrator for which f (x) ∼ x [MacNeil
and Eliasmith 2011].
Indeed learning both the feedforward and recurrent weights simultaneously using gradient descent on these loss functions, requires two different and unavailable error currents to be projected
into the postsynaptic neuron to make the rule local.

7.3

General dynamical system

General dynamical systems of the form
d~x(t) ~
= h(~x(t), ~u(t)),
dt
~ x(t))
~y (t) = K(~
can be learned using a different network configuration than the Figure 1A configuration used
for systems of the form (1). Here, the state variable is ~x, but the observable which serves as the
reference to the network is ~y . We handle this general system in two steps, first by absorbing the
transformation equation of the observable (second equation) into the first equation, and then
using an augmented variable and a different network structure.
Consider the transformation equation for the observable. The dimensionality of the relevant
variables: (1) the state variables (say joint angles and velocities) ~x; (2) the observables represented in the brain (say sensory representations of the joint angles and velocities) ~y ; and (3)
the control input (motor command) ~u, can be different from each other, but must be small
compared to the number of neurons. Furthermore, we require the observable ~y to not lose infor~ must be invertible, so ~y will have at least the same dimension as
mation compared to ~x, i.e. K
~x.
The time evolution of the observable is
ẏβ =

X ∂Kβ (~x)
α

∂xα

ẋα =

X ∂Kβ (~x)
∂xα

α

hα (~x, ~u) ≡ pβ (~y , ~u).

~ is invertible, so that ~x = K
~ −1 (~y ). So we essentially need
The last step follows since function K
to learn ẏβ = pβ (~y , ~u).
Having solved the observable transformation issue, we use ~x now for our observable, consistent
with the main text. The dynamical system to be learned is now ẋβ = hβ (~x, ~u). Consider an
augmented vector of the state and input variables x̃γ ≡ [xα , uβ ], where the index γ runs over
the values for indices α and β serially. Its time derivative is x̃˙ γ = [ẋα , u̇β ] = [hα (~x̃), u̇β ]. The α
components are of the same form as (1), but the β components involve u̇β which is not specified
as a function hβ (~x̃). However, we do not need u̇β as only the state variables need to be predicted,
not the input.
36

Our network for the general dynamical system, as shown in Supplementary Figure S4, decodes
the augmented ~x̃ˆ i.e. ~x̂ and ~û both as output. The error in the augmented ~x̃ is fed back to the
network. Since the hα (~x̃) is a mixed function of the state and input variables, ~u is no longer fed
into the network via a command representation layer, rather it enters only via the error in the
augmented variable. Once learning is complete, the error feedback in ~x̂ can be stopped, but the
error in ~û must still be fed back, so that ~u functions as a motor command.

7.4

Approximation error causes drift in weights

A frozen noise term ξ(~x(t)) due to the approximate decoding from non-linear tuning curves
of neurons, by the feedforward weights, recurrent weights and output decoders, will appear
additionally in equation (18). If this frozen noise has a non-zero mean over time as ~x(t) varies,
leading to a non-zero mean error, then it causes a drift in the weights due to the error-based
learning rules in equations (2), and possibly a consequent increase in error. Note that the
stability and error tending to zero proofs assume that this frozen noise is negligible.
Multiple strategies with contrasting pros and cons have been proposed to counteract this
parameter drift in the robust adaptive control literature [Ioannou and Sun 2012; Narendra and
Annaswamy 1989; Ioannou and Fidan 2006]. These include a weight leakage / regularizer term
switched slowly on, when a weight crosses a threshold [Ioannou and Tsakalis 1986; Narendra
and Annaswamy 1989], or a dead zone strategy with no updation of weights once the error is
lower than a set value [Slotine and Coetsee 1986; Ioannou and Sun 2012]. In our simulations,
the error continued to drop even over longer than typical learning time scales (Figure 5), and
so, we did not implement these strategies.
In practice, the learning can be stopped once error is low enough, while the error feedback
can be continued, so that the learned system does not deviate too much from the observed one.

37

7.5

Supplementary Figures

Firing rate (Hz)

A 200

0

B 400

1

Input u

0

1

1

Input u

1

Figure S1: Tuning curves of heterogeneous neurons The neurons in the command representation layer and
the recurrent network were heterogeneous with different tuning curves
P (random gains/encoders/biases). A,B.
Illustrative tuning curves of a few neurons, to projected input ui = α eiα uα into neuron i, are plotted, with
maximal firing rate of A. 200 Hz as for the network of Figure 2; and B. 400 Hz as for all other simulations (cf.
Supplementary Figure S2).

38

|

u2

Ai

0.20
0.00
0.15
feedback off | feedback on
5
0
4

| noise

end

Testing
Aii pulse on pedestal

feedback on |feedback off

feedback off

Bii

x2 , x̂ 2

0

10 -2

time (s)

8

time (s)

10000

9980

E 1500

time (s)

9988

rate (Hz)

400

time (s)

1.25

0
3

0

time (s)

8
¦

4

x1 , x̂ 1

4

0

rate (Hz)

400

1
0.50

time (s)

1.25

F1500

­

err 2

®

count

Nd , t

D

5
0
4

Cii 3

¦

count

Ci

error ²2

x2 ,x̂ 2

Bi

start

Learning

0

0

0

50

0.6 weight (arb·10−3 )

0.6

0

I

50
neuron #

25

H

density

G

0

neuron #

10 -6

1
0.50

Figure S2: Learning van der Pol oscillator dynamics via FOLLOW with double firing rates. Layout
and legend of panels A-I are analogous to Figure 2A-I, except that the gains of neurons were doubled yielding
approximately double firing rates, but one-fifth the learning time.

39

u3

Ai

| startLearning end | zero

0.06

Testing
Aii

zero

0.00
0.06
feedback on |feedback off

feedback off

Bii

0
25
30

Cii 20

0

0

30

0

40

time (s)

14900

14940

time (s)

x̂3

x2

x1

x̂ 2

x̂1

20

E 20

x3

D

feedback off | feedback on

x2 , x̂ 2

error ²3

Ci

20

0

time (s)

40

20

x1 , x̂ 1

25

max n (x3 )

20

max n + 1 (x3 )

x3 ,x̂ 3

Bi

2

2

Figure S3: Learning the Lorenz system without filtering the reference variables.
Panels A-E are interpreted similar to Figure 3A-E, except that the reference signal was not filtered in computing
the error. Without filtering, the tent map for the network output (panel E red) shows a doubling, but the
mismatch for large maxima is reduced, compared to with filtering in Figure 3.

recurrent network

predicted state
& input
ˆ γ (t) = (~x̂, ~û)
x̃

Si

wij

(x̃)

dγi

keiγ
error ˜γ

−
+

reference state
& input
x̃γ (t) = (~x, ~u)

Figure S4: Network configuration for learning a general dynamical system
ˆγ = (~x̂, ~û). Compare this configuration with that in
The recurrent network decodes the augmented vector x̃
Figure 1A. The input is no longer provided directly to the network but only via the augmented error signal. The
augmented error signal is fed back into the recurrent layer. Only the recurrent weights need to be learned, as
there is no feedforward layer here. After learning is turned off, the error in input must continue to be fed-back
to the neurons, to serve as a motor command; but the error in the state variables need not be.

40

Ai

|

u1

0.10

start

Learning

| noise

end

Testing
ramp and step
Aii

0.00
0.06

x1 ,x̂ 1

Bi

feedback off | feedback on
0.6
0.0

feedback on |feedback off

feedback off

Bii

0.8

Ci

x2 , x̂ 2

0

time (s)

E

®

err 2
­

10 -9

9990

time (s)

9994

time (s)

10000

2
0

weight (arb·10−3 )

10 -3

4

0

time (s)

10000

7

0.0
0.3

F 12

time (s)

4

0.2

x1 , x̂ 1

0.5

1

weight (arb·10−3 )

1

density

error ²1

0.0
0.1

Nd , t

D

Cii0.4

0.5

0

0

0

Figure S5: Learning linear dynamics via FOLLOW: 2D decaying oscillator. Layout and legend of panels
A-D are analogous to Figure 2A-D, except an error component is shown in the right-most panel of C. E. A few
randomly selected weights are shown evolving during learning. F. Histogram of weights after learning. A few
strong weights |wij | > 1 are out of bounds and not shown here (cf. panel E).

41

Ai

|

u2

0.35

start

Learning

| noise

end

Testing
ramp and step
Aii

0.00
0.30
feedback off | feedback on
0.8

x2 ,x̂ 2

Bi

0.0
0.4

Cii0.10

0.15
0.00
0

time (s)

E

time (s)

9994

time (s)

10000

3

®

err 2

0

­

10 -9

9990

weight (arb·10−3 )

10 -4

4

0

time (s)

10000

1

0

time (s)

4

0.3

x1 , x̂ 1

0.3

0.00
0.35

F 70
density

0.30

Nd , t

D

feedback off

Bii

x2 , x̂ 2

error ²2

Ci

feedback on |feedback off

0

0

0.3 weight (arb·10−3 )

0.3

Figure S6: Learning nonlinear feedforward transformation with linear recurrent dynamics via FOLLOW. Panels A-F are interpreted similar to Figure S5A-F, except in panel A, along with the input u2 (t) (blue)
to layer 1, the required non-linear transform g2 (~u(t))/20 is also plotted (cyan); and in panels E and F, the evolution and histogram of weights are those of the feedforward weights, that perform the non-linear transform on
the input.

42

non-lin-sys wts (arb×10−4 )

A 50 feedforward B 50 recurrent
R 2 = -13.354

R 2 = 0.942

0
50

0
50

50
0
50
50−4
lin-sys wts (arb×10 )

0

50

Figure S7: Feedforward weights are uncorrelated, while recurrent ones are correlated, when learning
same recurrent dynamics but with different feedforward transforms.
The linear decaying oscillators system was learned for 10,000s with either a non-linear or a linear feedforward
ff
transform. A. The learned feedforward weights wij
were plotted for the system with the non-linear feedforward
transform, versus the corresponding feedforward weights for the system with the linear feedforward transform.
The feedforward weights for the two systems do not fit the identity line (coefficient of determination R2 is
negative; R2 is not the square of a number and can be negative) showing that the learned feedforward transform
is different in the two systems as expected. B. Same as A, but for the recurrent weights in the two systems. The
weights fit the identity line with an R2 close to 1 showing that the learned recurrent transform is similar in the
two systems as expected. Some weights fall outside the plot limits.

43

