Predicting nonlinear dynamics by stable local learning in a recurrent spiking neural network
Abstract
The brain needs to predict how the body reacts to motor commands, but how a network of spiking neurons can learn nonlinear body dynamics using local, online and stable learning rules is unclear. 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. The rule for Feedbackbased 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, nonlinear and chaotic dynamics, as well as the dynamics of a twolink arm. Under reasonable approximations, we show, using the Lyapunov method, that FOLLOW learning is uniformly stable, with the error going to zero asymptotically.
https://doi.org/10.7554/eLife.28295.001Introduction
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 Ross 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 prenatal (Khazipov et al., 2004) and postnatal 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 nonlinear 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 nonlinear dynamics, given command input, is known to be difficult in networks of rate units, and even more so in networks of spiking neurons (Abbott et al., 2016). Ideally, in order to be biologically plausible, a learning rule must be online that is constantly incorporating new data, as opposed to batch learning where weights are adjusted only after many examples have been seen; and local that is 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 realtime recurrent learning (RTRL) (Williams and Zipser, 1989) which are nonlocal in time or in space, respectively (Pearlmutter, 1995; Jaeger, 2005). Even though LongShortTermMemory (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), function and dynamics approximation (Funahashi, 1989; Hornik et al., 1989; Girosi and Poggio, 1990; Sanner and Slotine, 1992; Funahashi and Nakamura, 1993; Pouget and Sejnowski, 1997; Chow and XiaoDong Li, 2000; Seung et al., 2000; Eliasmith and Anderson, 2004; Eliasmith, 2005) and adaptive control theory (Morse, 1980; Narendra et al., 1980; Slotine and Coetsee, 1986; Weiping Li et al., 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) where learning was restricted to the readout connections, we focus on a learning rule for the recurrent connections. Whereas neural network implementations of control theory (Sanner and Slotine, 1992; DeWolf et al., 2016) modified adaptive feedback weights without a synaptically local interpretation, we modify the recurrent weights in a synaptically local manner. Compared 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 Denève, 2015), for nonlinear 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 lowdimensional nonlinear 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 errorfeedback weights, the learning rule is synaptically local and combines presynaptic activity with the local postsynaptic error variable.
Results
A forward predictive model (Figure 1A) takes, at each time step, a motor command $\overrightarrow{u}(t)$ as input and predicts the next observable state $\overrightarrow{\widehat{x}}(t+\mathrm{\Delta}t)$ of the system. In the numerical implementation, we consider $\mathrm{\Delta}t=1$ms, but for the sake of notational simplicity we drop the $\mathrm{\Delta}t$ in the following. The predicted system state $\overrightarrow{\widehat{x}}$ (e.g., the vector of joint angles and velocities of the arm) is assumed to be lowdimensional with dimensionality ${N}_{d}$ (4dimensional for a twolink arm). The motor command $\overrightarrow{u}(t)$ is used to generate target movements such as ‘lift your arm to a location’, with a dimensionality ${N}_{c}$ of the command typically smaller than the dimensionality ${N}_{d}$ of the system state.
The actual state of the reference system (e.g., actual joint angles and velocities of the arm) is described by a nonlinear dynamical system, which receives the control input $\overrightarrow{u}(t)\in {\mathbb{R}}^{{N}_{c}}$ and evolves according to a set of coupled differential equations
where $\overrightarrow{x}$ with components ${x}_{\alpha}$ (where $\alpha =1,\mathrm{\dots},{N}_{d}$) is the vector of observable state variables, and $\overrightarrow{h}$ is a vector whose components are arbitrary nonlinear functions ${h}_{\alpha}$. For example, the observable system state $\overrightarrow{x}(t)$ could be the joint angles and velocities of the arm deduced from visual and proprioceptive input (Figure 1A). We show that, with training, the forward predictive model learns to make the error
between the actual state $\overrightarrow{x}(t)$ and the predicted state $\overrightarrow{\widehat{x}}(t)$ negligible.
Network architecture for learning the forward predictive model
In our neural network model (Figure 1B), the motor command $\overrightarrow{u}(t)$ drives the spiking activity of a command representation layer of 3000 to 5000 leaky integrateandfire neurons via connections with fixed random weights. These neurons project, via plastic feedforward connections, to a recurrent network of also 3000 to 5000 integrateandfire neurons. We assume that the predicted state $\overrightarrow{\widehat{x}}$ is linearly decoded from the activity of the recurrent network. Denoting the spike train of neuron $i$ by ${S}_{i}(t)$, the component $\alpha $ of the predicted system state is
where ${d}_{\alpha i}$ are the readout weights. The integral represents a convolution with a lowpass filter
with a time constant ${\tau}_{s}=20$ ms, and is denoted by $({S}_{i}*\kappa )(t)$.
The current into a neuron with index $l$ ($l=1,\mathrm{\dots},N$), in the command representation layer comprising $N$ neurons, is
where $e}_{l\alpha}^{\mathrm{f}\mathrm{f}$ are fixed random weights, while $b}_{l}^{\mathrm{f}\mathrm{f}$ is a neuronspecific constant for bias (see Methods) (Eliasmith and Anderson, 2004). We use Greek letters for the indices of lowdimensional variables (such as command) 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 is much larger than the dimensionality of the input, that is $N\gg {N}_{c}$.
The input current to a neuron with index $i$ ($i=1,\mathrm{\dots},N$) in the recurrent network is
where $w}_{il}^{\mathrm{f}\mathrm{f}$ and ${w}_{ij}$ are the feedforward and recurrent weights, respectively, which are both subject to our synaptic learning rule, whereas $k{e}_{i\alpha}$ are fixed error feedback weights (see below). The spike trains travelling along the feedforward path $S}_{l}^{\mathrm{f}\mathrm{f}$ and those within the recurrent network ${S}_{j}$ are both lowpass filtered (convolution denoted by $*$) at the synapses with the exponential filter $\kappa $ defined above. The constant parameter ${b}_{i}$ is a neuron specific bias (see Methods). The constant $k>0$ is the gain for feeding back the output error. The number of neurons $N$ in the recurrent network is much larger than the dimensionality ${N}_{d}$ of the represented variable $\widehat{x}$, that is $N\gg {N}_{d}$.
For all numerical simulations, we used deterministic leaky integrate and fire (LIF) neurons. The voltage ${V}_{l}$ of each LIF neuron indexed by $l$, was a lowpass filter of its driving current ${J}_{l}$:
with a membrane time constant, of ${\tau}_{m}=20$ ms. The neuron fired when the voltage ${V}_{l}$ crossed a threshold $\theta =1$ from below, after which the voltage was reset to zero for a refractory period ${\tau}_{r}$ of 2 ms. If the voltage went below zero, it was clipped to zero. Mathematically, the spike trains ${S}_{l}^{\mathrm{f}\mathrm{f}}(t)$ in the command representation layer and ${S}_{l}(t)$ in the recurrent network, are a sequence of events, modelled as a sum of Dirac deltafunctions.
Biases and input weights of the spiking neurons vary between one neuron and the next, both in the command representation layer and the recurrent network, yielding different frequency versus input curves for different neurons (Figure 1—figure supplement 1). Since arbitrary lowdimensional functions can be approximated by linear decoding from a basis of nonlinear functions (Funahashi, 1989; Girosi and Poggio, 1990; Hornik et al., 1989), such as neuronal tuning curves (Sanner and Slotine, 1992; Seung et al., 2000; Eliasmith and Anderson, 2004), we may expect that suitable feedforward weights onto, and lateral weights within, the recurrent network can be found that approximate the role of the function $\overrightarrow{h}$ in Equation (1). In the next subsection, we propose an error feedback architecture along with a local and online synaptic plasticity rule that can train these feedforward and recurrent weights to approximate this role, while the readout weights are kept fixed, so that the network output mimics the dynamics in Equation (1).
Negative error feedback via autoencoder enables local learning
To enable weight tuning, we make four assumptions regarding the network architecture. The initial two assumptions are related to input and output. First, we assume that, during the learning phase, a random timedependent motor command input $\overrightarrow{u}(t)$ is given to both the musclebody reference system described by Equation (1) and to the spiking network. The random input generates irregular trajectories in the observable state variables, mimicking motor babbling (Meltzoff and Moore, 1997; Petersson et al., 2003). Second, we assume that each component ${\widehat{x}}_{\alpha}$ of the output predicted by the spiking network is compared to the actual observable output ${x}_{\alpha}$ produced by the reference system of Equation (1) and their difference (the output error ${\u03f5}_{\alpha}$; Equation (2)) is calculated, similar to supervised learning schemes such as perceptron learning (Rosenblatt, 1961).
The final two assumptions are related to the error feedback. Our third assumption is that the readout weights ${d}_{\alpha i}$ have been prelearned, possibly earlier in development, in the absence of feedforward and recurrent connections, so as to form an autoencoder of gain $k$ with the fixed random feedback weights $k{e}_{i\alpha}$. Specifically, an arbitrary value ${\u03f5}_{\alpha}$ 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{\u03f5}_{\alpha}$. Thus, we set the decoding weights so as to minimize the squared error between the decoded output and required output $k\overrightarrow{\u03f5}$ for a set of randomly chosen vectors $\overrightarrow{\u03f5}$ while setting feedforward and recurrent weights to zero (see Methods). We used an algorithmic learning scheme here, but we expect that these decoding weights can also be prelearned by biologically plausible learning schemes (D'Souza et al., 2010; Urbanczik and Senn, 2014; Burbank, 2015).
Fourth, we assumed that the error ${\u03f5}_{\alpha}={x}_{\alpha}{\widehat{x}}_{\alpha}$ is projected back to neurons in the recurrent network through the abovementioned fixed random feedback weights. From the third term in Equation (6) and Figure 1B–C, we define a total error input that neuron $i$ receives:
with feedback weights $k{e}_{i\alpha}$, where $k$ is fixed at a large constant positive value.
The combination of the autoencoder and the error feedback implies that the output stays close to the reference, as explained now. In open loop that is without connecting the output $\overrightarrow{\widehat{x}}$ and the reference $\overrightarrow{x}$ to the error node, an input $\overrightarrow{\u03f5}$ to the network generates an output $\overrightarrow{\widehat{x}}=k\overrightarrow{\u03f5}$ due to the autoencoder of gain $k$. In closed loop, that is with the output and reference connected to the error node (Figure 1B), the error input is $\overrightarrow{\u03f5}=\overrightarrow{x}\overrightarrow{\widehat{x}}$, and the network output $\overrightarrow{\widehat{x}}$ settles to:
that is approximately the reference $\overrightarrow{x}$ for large positive $k$. The fedback residual error $\overrightarrow{\u03f5}=\overrightarrow{x}/(k+1)$ drives the neural activities and thence the network output. Thus, feedback of the error causes the output ${\widehat{x}}_{\alpha}$ to approximately follow ${x}_{\alpha}$, for each component $\alpha $, as long as the error feedback time scale is fast compared to the reference dynamical system time scale, analogous to negative error feedback in adaptive control (Narendra and Annaswamy, 1989; Ioannou and Sun, 2012).
While error feedback is on, the synaptic weights $w}_{il}^{\mathrm{f}\mathrm{f}$ and ${w}_{ij}$ on the feedforward and recurrent connections, respectively, are updated as:
where $\eta $ is the learning rate (which is either fixed or changes on the slow time scale of minutes), and ${\kappa}^{\u03f5}$ is an exponentially decaying filter kernel with a time constant of 80 or 200 ms. For a postsynaptic neuron $i$, the error term ${I}_{i}^{\u03f5}*{\kappa}^{\u03f5}$ is the same for all its synapses, while the presynaptic contribution is synapsespecific.
We call the learning scheme ‘Feedbackbased Online Local Learning Of Weights’ (FOLLOW), since the predicted state $\overrightarrow{\widehat{x}}$ follows the true state $\overrightarrow{x}$ from the start of learning. Under precise mathematical conditions, we show in the Methods that the FOLLOW scheme converges to a stable solution, while simultaneously deriving the learning rule.
Because of the error feedback, with constant $k\gg 1$, the output is close to the reference from the start of learning. However, initially the error is not exactly zero, and this nonzero error drives the weight updates via Equation (10). After a sufficiently long learning time, a vanishing error (${\u03f5}_{\alpha}=0$ for all components) indicates that the neuronal network now autonomously generates the desired output, so that feedback is no longer required. In the Methods section, we show that not just the lowdimensional output $\overrightarrow{\widehat{x}}$, but also the spike trains ${S}_{i}(t)$, for $i=1,\mathrm{\dots},N$, are entrained by the error feedback to be close to the ideal ones required to generate $\overrightarrow{x}$.
During learning, the error feedback via the autoencoder in a loop serves two roles: (i) to make the error current available in each neuron, projected correctly, for a local synaptic plasticity rule, and (ii) to drive the spike trains to the target ones for producing the reference output. In other learning schemes for recurrent neural networks, where neural activities are not constrained by error feedback, it is difficult to assign credit or blame for the momentarily observed error, because neural activities from the past affect the present output in a recurrent network. In the FOLLOW scheme, the spike trains are constrained to closely follow the ideal time course throughout learning, so that the present error can be attributed directly to the weights, enabling us to change the weights with a simple perceptronlike learning rule (Rosenblatt, 1961) as in Equation (10), bypassing the credit assignment problem. In the perceptron rule, the weight change $\mathrm{\Delta}w\sim (\text{pre})\cdot \delta $ is proportional to the presynaptic input $(\text{pre})$ and the error $\delta $. In the FOLLOW learning rule of Equation (10), we can identify $({S}_{i}*\kappa )$ with $(\text{pre})$ and $({I}_{i}^{\u03f5}*{\kappa}^{\u03f5})$ with $\delta $. In Methods, we derive the learning rule of Equation (10) in a principled way from a stability criterion.
FORCE learning (Sussillo and Abbott, 2009, 2012; DePasquale et al., 2016; Thalmeier et al., 2016; Nicola and Clopath, 2016) also clamps the output and neural activities to be close to ideal during learning, by using weight changes that are faster than the time scale of the dynamics. In our FOLLOW scheme, clamping is achieved via negative error feedback using the autoencoder, which allows weight changes to be slow and makes the error current available locally in the postsynaptic neuron. Other methods used feedback based on adaptive control for learning in recurrent networks of spiking neurons, but were limited to linear systems (MacNeil and Eliasmith, 2011; Bourdoukan and Denève, 2015), whereas the FOLLOW scheme was derived for nonlinear systems (see Methods). Our learning rule of Equation (10) uses an error ${\u03f5}_{\alpha}\equiv {x}_{\alpha}{\widehat{x}}_{\alpha}$ in the observable state, rather than an error involving the derivative $d{x}_{\alpha}/dt$ in Equation (1), as in other schemes (see Appendix 1) (Eliasmith, 2005; MacNeil and Eliasmith, 2011). The reader is referred to Discussion for detailed further comparisons. The FOLLOW learning rule is local since all quantities needed on the righthandside of Equation (10) could be available at the location of the synapse in the postsynaptic neuron. For a potential implementation and prediction for errorbased synaptic plasticity, and 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 (Figure 2), lowdimensional chaos (Figure 3) and simulated arm movements (Figure 4) (additional examples in Figure 2—figure supplement 2, Figure 2—figure supplement 4 and Methods). In all simulations, we started with vanishingly small feedforward and recurrent weights (tabula rasa), but assumed prelearned 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 timevarying input to both the network (Figure 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 Equation (10).
During learning, the mean squared error, where the mean was taken over the number of dynamical dimensions ${N}_{d}$ and over a duration of a few seconds, decreased (Figure 2D). We stopped the learning phase that is weight updating, when the mean squared error approximately plateaued as a function of learning time (Figure 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 openloop model over a finite time horizon (corresponding to the planning horizon of a short action sequence) and in the closedloop mode (with error feedback) without time limit.
Nonlinear 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 nonlinear 2dimensional van der Pol oscillator (Figure 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 four seconds before learning started, we blocked error feedback. Because of zero error feedback and our initialization with zero feedforward and recurrent weights, the output $\widehat{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 ($k=10$) 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 (Figure 2B,C). During learning, the error dropped even further over time (Figure 2D). After having stopped learning at 5000 s ($\sim $2 hr), we found the weight distribution to be unimodal with a few very large weights (Figure 2G). In the openloop 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 (Figure 2Bii), but the shape of the nonlinear oscillation was well predicted as indicated by the similarity of the trajectories in state space (Figure 2Cii).
The spiking pattern of neurons of the recurrent network changed as a function of time, with interspike intervals of individual neurons correlated with the output, and varying over time (Figure 2H,I). The distributions of firing rates averaged over a 0.25 s period with fairly constant output, and over a 16 s period with timevarying output, were longtailed, with the mean across neurons maintained at approximately 12–13 Hz (Figure 2E,F). The distribution averaged over 16 s had a smaller number of neurons firing at very low and very high rates compared to the distribution over 0.25 s, consistent with the expectation that the identity of lowrate and highrate neurons changed over time for timevarying output (Figure 2E,F). We repeated this example experiment (‘van der Pol oscillator’) with a network of equal size but with neurons that had higher firing rates, so that some neurons could reach a maximal rate of 400 Hz (Figure 1—figure supplement 1). The reference was approximated better and learning time was shorter with higher rates (Figure 2—figure supplement 1 – 10,000 s with constant learning rate) compared to the low rates here (Figure 2 – 5,000 s with 20 times the learning rate after 1,000 s). Hence, for all further simulations, we set neuronal parameters to enable peak firing rates up to 400 Hz (Figure 1—figure supplement 1B).
We also asked whether merely the distribution of the learned weights in the recurrent layer was sufficient to perform the task, or whether the specific learned weight matrix was required. This question was inspired from 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), where the recurrent weights are random, and only the readout weights are learned. To answer this question, we implemented a perceptron learning rule on the readout weights initialized at zero, with the learned network’s output as the target, after setting the feedforward and/or recurrent weights to either the learned weights as is or after shuffling them. The readout weights could be approximately learned only for the network having the learned weights and not the shuffled ones (Figure 2—figure supplement 3), supporting the view that the network does not behave like a reservoir (Methods).
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 3dimensional nonlinear chaotic Lorenz system (Figure 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 ${\u03f5}_{3}(t)$ (Figure 3A–C). Such a behaviour 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 qualitatively the same underlying strange attractor dynamics, as seen from the butterfly shape (Lorenz, 1963) of the attractor in configuration space, and the tent map (Lorenz, 1963) of successive maxima versus the previous maxima (Figure 3D,E). The tent map generated from our network dynamics (Figure 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 exponentiallydecaying filter, we repeated learning without filtering the Lorenz reference signal (Figure 3—figure supplement 1), and found that the mismatch for large maxima reduced, but a doubling appeared in the tent map (Figure 3—figure supplement 1E) which had been almost imperceptible with filtering (cf. Figure 3E).
FOLLOW enables learning a twolink 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 twolink arm via the FOLLOW scheme. We used a twolink 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}^{\circ}$ to $+{90}^{\circ}$ where the resting state is at ${0}^{\circ}$ (see Methods).
The dynamical system representing the arm is fourdimensional 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. Still, our feedforward and recurrent network architecture (Figure 1B) with 5000 neurons in each layer was able to approximate these dynamics.
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 converged rapidly to the target trajectory (Figure 4). We found that the FOLLOW scheme learned to reproduce the arm dynamics even without error feedback for a few seconds during the test phase (Figure 4 and Video 1 and Video 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 post learning, and tested the network in the openloop setting on a reaching task and an acrobotinspired swinging task (Sutton, 1996). In the reaching task, torque was provided to both joints to enable the armtip 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 Video 1. We also tested the learned network model of the 2link arm on an acrobotlike task that is a gymnast swinging on a highbar (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 the target $(x,y)$ position from rest by swinging. The control input and the dynamics (Figure 4A–C right panels, Figure 4E and Video 2) show that the network can perform the task in openloop condition suggesting that it has learned the inertial properties of the arm model, necessary for this simplified acrobot task.
Feedback in the FOLLOW scheme entrains spike timings
In Methods, 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, that is 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 nonlinear tuning curves, the nonlinear function $\overrightarrow{h}$ 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, that is 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 pseudorandom 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 grow beyond any reasonable bound. Ameliorative techniques from adaptive control are discussed in Appendix 1. 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 (Figure 5A), 10 times longer than that required for learning this example (Figure 2—figure supplement 1).
To highlight the difference between a realizable reference system and nonlinear differential equations as a reference system, we used, in an additional simulation experiment, a spiking network with fixed weights as the reference. More precisely, instead of using directly the differential equations of the van der Pol oscillator as a reference, we now used as a reference a spiking approximation of the van der Pol oscillator, that is the spiking network that was the final result after 10,000 s ($\sim $3 hr) of FOLLOW learning in Figure 2—figure supplement 1. For both the spiking reference network and the tobetrained learning network we used the same architecture, the same number of neurons, and the same neuronal parameters as in Figure 2—figure supplement 1 for the learning of the van der Pol oscillator. The readout and feedback weights of the learning network also 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 ($\sim $28 hr) (Figure 5). With the realizable network as reference, learning was more rapid than with the original van der Pol oscillator as reference (Figure 5A).
We emphasize that, analogous to the earlier simulations, the feedback error ${\u03f5}_{\alpha}$ was lowdimensional and calculated from the decoded outputs. Nevertheless, the lowdimensional error feedback was able to entrain the network spike times to the reference spike times (Figure 5C). In particular, a few neurons learned to fire 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 ($\sim $28 hr) of learning. The spike trains were entrained by the lowdimensional feedback. With the feedback off, even the lowdimensional output, and hence the spike trains, diverged from the reference. It will be interesting to explore if this entrainment by lowdimensional feedback via an autoencoder loop can be useful in supervised spike train learning (Gütig and Sompolinsky, 2006; Pfister et al., 2006; Florian, 2012; Mohemmed et al., 2012; Gütig, 2014; Memmesheimer et al., 2014; Gardner and Grüning, 2016).
Our results with the spiking reference network suggest that the error is reduced to a value close to zero for a realizable or closelyapproximated system (Figure 5A) as shown in Methods, analogous to proofs in adaptive control (Ioannou and Sun, 2012; Narendra and Annaswamy, 1989). 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 should carry over to our learning scheme.
Learning is robust to sparse connectivity, noisy error or reference, and noisy decoding weights, but not to delays
So far, our spiking networks had alltoall 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 lowdimensional 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; 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.
We added Gaussian white noise to each component of the error, which is equivalent to adding it to each component of the reference, and ran the van der Pol oscillator learning protocol for 10,000 s for different standard deviations of the noise (Figure 6B). The learning was robust to noise with standard deviation up to around $0.001$, which must be compared with the error amplitude of the order of $0.1$ at the start of learning, and orders of magnitude lower later.
The readout weights have been prelearned until now, so that, in the absence of recurrent connections, error feedback weights and decoding weights formed an autoencoder. We sought to relax this requirement. Simulations showed that with completely random readout weights, the system did not learn to reproduce the target dynamical system. However, if the readout weights had some overlap with the autoencoder, learning was still possible (Figure 6C). If for a feedback error $\overrightarrow{\u03f5}$, the error encoding followed by output decoding yields $k(1+\xi )\overrightarrow{\u03f5}+\overrightarrow{n}(\overrightarrow{\u03f5})$, where $\overrightarrow{n}$ is a vector of arbitrary functions not having linear terms and small in magnitude compared to the first term, and $\xi $ is sufficiently greater than $1$ so that the effective gain $k(1+\xi )$ remains large enough, then the term that is 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 autoencoder by one plus $\gamma $, where for each weight $\gamma $ was drawn independently from a uniform distribution between $\chi +\xi $ and $\chi +\xi $. We found that the system was still able to learn the van der Pol oscillator up to $\chi \sim 5$ and $\xi =0$, or $\chi =2$ and $\xi $ variable (Figure 6B,C). Negative values of $\xi $ result in a lower overlap with the autoencoder leading to the asymmetry seen in Figure 6C. Thus, the FOLLOW learning scheme is robust to multiplicative noise on the decoding weights. Alternative approaches for other noise models are discussed in Appendix 1.
We also asked if the network could handle sensory feedback delays in the reference signal. Due to the strong limit cycle attractor of the van der Pol oscillator, the effect of delay is less transparent than for the linear decaying oscillator (Figure 2—figure supplement 2), so we decided to focus on the latter. For the linear decaying oscillator, we found that learning degraded rapidly with a few milliseconds of delay in the reference, that is if $\overrightarrow{x}(t\mathrm{\Delta})$ was provided as reference instead of $\overrightarrow{x}(t)$ (Figure 6E–F). We compensated for the sensory feedback delay by delaying the motor command input by identical $\mathrm{\Delta}$ (Figure 6G), which is equivalent to timetranslating the complete learning protocol, to which the learning is invariant, and thus the network would learn for arbitrary delay (Figure 6H). In the Discussion, we suggest how a forward model learned with a compensatory delay (Figure 6G) could be used in control mode to compensate for sensory feedback delays.
Discussion
The FOLLOW learning scheme enables a spiking neural network to function as a forward predictive model that mimics a nonlinear dynamical system activated by one or several timevarying inputs. The learning rule is supervised, local, and comes with a proof of stability.
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 implementations 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.
A possible implementation in a spatially extended neuron would be to imagine that the postsynaptic error current ${I}_{i}^{\u03f5}$ 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 7A. Consistent with the picture of a messenger molecule, we lowpass filtered the error current with an exponential filter ${\kappa}^{\u03f5}$ of time constant 80 ms or 200 ms, much longer than the synaptic time constant of 20 ms of the filter $\kappa $. Simultaneously, filtered information about presynaptic spike arrival ${S}_{j}*\kappa $ is available at each synapse, possibly in the form of glutamate bound to the postsynaptic receptor or by calcium triggered signalling 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.
The separation of the error current from the currents at feedforward or recurrent synapses could be spatial (such as suggested in Figure 7A) or chemical if the error current projects onto synapses that trigger a signalling 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. This leads us to a prediction regarding synaptic plasticity that, say in cortical pyramidal neurons, the plasticity of synapses that are driven by presynaptic input in the basal dendrites, should be modulated by currents injected in the apical dendrite or on stimulation of feedback connections.
The learning scheme is provenly stable with errors converging asymptotically to zero under a few reasonable assumptions (Methods). The first assumption is that error encoding feedback weights and output decoding readout weights form an autoencoder. This requirement can be met if, at an early developmental stage, either both sets of weights are learned using say mirrored STDP (Burbank, 2015), or the output readout weights are learned, starting with random encoding weights, via a biological perceptronlike learning rule (D'Souza et al., 2010; Urbanczik and Senn, 2014). A prelearned autoencoder in a highgain negative feedback loop is in fact a specific prediction of our learning scheme, to be tested in systemslevel experiments. The second assumption is that the reference dynamics $f(\overrightarrow{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 $\overrightarrow{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 higherdimensional space with variables $\overrightarrow{y}(t)$ as long as $\overrightarrow{y}=K(\overrightarrow{x})$ is an invertible function of $\overrightarrow{x}(t)$ (Appendix 1). The fourth assumption is that the system dynamics be slower than synaptic dynamics. Indeed, typical reaching movements extend over hundreds of milliseconds 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 openloop condition when error feedback was switched off (Figures 2 and 4). The fifth assumption is that weights stay bounded. Indeed, in biology, synaptic weights should not grow indefinitely. Algorithmically, a weight decay term in the learning rule can suppress the growth of large weights (see also Appendix 1), though we did not need to implement a weight decay term in our simulations.
One of the postulated uses of the forward predictive model is to compensate for delay in the sensory feedback during motor control (Wolpert and Miall, 1996; Wolpert et al., 1995) using the Smith predictor configuration (Smith, 1957). We speculate that the switch from the closedloop learning of forward model with feedback gain $k\gg 1$ to openloop motor prediction $k=0$ could also be used to switch delay lines: the system can have either a delay before the forward model as required for learning (Figure 7B), or after the forward model as required for the Smith predictor (Figure 7C). We envisage that FOLLOW learning of the forward model occurs in closed loop mode ($k\gg 1$) with a delay in the motor command path, as outlined earlier in Figure 6G and now embedded in the Smith predictor architecture in Figure 7B. After learning, the network is switched to motor control mode, with the forward predictive model in open loop ($k=0$), implementing the Smith predictor (Figure 7C). In this motor control mode, the motor command is fed with zero delay to the forward model. This enables to rapidly feed the estimated state back to the motor controller so as to take corrective actions, even before sensory feedback arrives. In parallel, available sensory feedback is compared with a copy of the forward model that has passed through a compensatory delay after the forward model (Figure 7C).
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 local 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 nonlinear dynamics (Appendix 1). Denève and colleagues used error feedback in a homogeneous spiking network with a rule similar to ours, for linear dynamics only (Bourdoukan and Denève, 2015), and while this article was in review, also for nonlinear dynamics (Alemi et al., 2017), but their network requires instantaneous lateral interactions and in the latter case, also nonlinear dendrites.
Reservoir computing models exploit recurrent networks of nonlinear 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 timedependent 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 multidimensional target, multidimensional 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 synaptic 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 nonlinear systems using continuous rate neurons (Sanner and Slotine, 1992; Weiping Li et al., 1987; Slotine and Coetsee, 1986) or spiking neurons (DeWolf et al., 2016) has primarily focused on learning parameters in adaptive feedback paths, rather than learning weights in a recurrent network, using learning rules involving quantities that do not appear in the pre or postsynaptic neurons, making them difficult to interpret as local to synapses. Recurrent networks of rate units have occasionally been used for control (Zerkaoui et al., 2009), but trained either via realtime recurrent learning or the extended Kalman filter which are nonlocal in space, or via backpropagation through time which is offline (Pearlmutter, 1995). Recent studies have used neural network techniques to train inverse models by motor babbling, to describe behavioral data in humans (Berniker and Kording, 2015) and song birds (Hanuschkin et al., 2013), albeit with abstract networks. 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 rewardmodulated 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 given timedependent control input remains open.
Additional features of the FOLLOW learning scheme are that it does not require full connectivity but also works with biologically more plausible sparse connectivity; and 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 lowdimensional output and all neural currents are spatially averaged over a large number of synapticallyfiltered spike trains, neurons in the FOLLOW network do not necessarily need to fire at rates higher than the inverse of the synaptic time scale. In conclusion, we used a network of heterogeneous neurons as in the Neural Engineering Framework (Eliasmith and Anderson, 2004), employed a prelearned autoencoder to enable negative feedback of error as in adaptive control theory (Morse, 1980; Narendra et al., 1980; Slotine and Coetsee, 1986; Weiping Li et al., 1987; Narendra and Annaswamy, 1989; Sastry and Bodson, 1989; Ioannou and Sun, 2012), and derived and demonstrated a local and online learning rule for recurrent connections that learn to reproduce nonlinear dynamics.
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 of excitatory and inhibitory incoming currents (Denè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 excitatoryinhibitory balance, possibly using inhibitory plasticity (Vogels et al., 2011); prelearns the autoencoder, potentially via mirrored STDP (Burbank, 2015); and possibly implements spatial separation between different compartments (Urbanczik and Senn, 2014). It would also be interesting for future work to see whether our model of an arm trained on motor babbling with FOLLOW, can explain aspects of human behavior in reaching tasks involving force fields (Shadmehr and MussaIvaldi, 1994), uncertainty (Körding and Wolpert, 2004; Wei and Körding, 2010) or noise (Burge et al., 2008). 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.
Methods
Simulation software
Request a detailed protocolAll simulation scripts were written in python (https://www.python.org/) for the Nengo simulator (Stewart et al., 2009) (http://www.nengo.ca/, version 2.4.0) 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 are available online at https://github.com/adityagilra/FOLLOW (Gilra, 2017). A copy is archived at https://github.com/elifesciencespublications/FOLLOW.
Network parameters
Initialization of plastic weights
Request a detailed protocolThe feedforward weights $w}_{il}^{\mathrm{f}\mathrm{f}$ from the command representation layer to the recurrent network and the recurrent weights ${w}_{ij}$ inside the network were initialized to zero.
Update of plastic weights
Request a detailed protocolWith the error feedback loop closed, that is with reference output $\overrightarrow{x}$ and predicted output $\widehat{\overrightarrow{x}}$ connected to the error node, and feedback gain $k=10$, the FOLLOW learning rule, Equation (10), was applied on the feedforward and recurrent weights, $w}_{il}^{\mathrm{f}\mathrm{f}$ and ${w}_{ij}$. The error for our learning rule was the error ${\u03f5}_{\alpha}={x}_{\alpha}{\widehat{x}}_{\alpha}$ in the observable output $\overrightarrow{x}$, not the error in the desired function $\overrightarrow{h}(\overrightarrow{x},\overrightarrow{u})$ (cf. Eliasmith, 2005; MacNeil and Eliasmith, 2011, Appendix 1). The observable reference state $\overrightarrow{x}$ was obtained by integrating the differential equations of the dynamical system. The synaptic time constant ${\tau}_{s}$ was 20 ms in all synapses, including those for calculating the error and for feeding the error back to the neurons (decaying exponential $\kappa $ with time constant ${\tau}_{s}$ in Equation (6)). The error used for the weight update was filtered by a 200 ms decaying exponential (${\kappa}^{\u03f5}$ in Equation (10)).
Random setting of neuronal parameters and encoding weights
Request a detailed protocolWe used leaky integrateandfire neurons with a threshold $\theta =1$ and time constant ${\tau}_{m}=20$ ms. After each spike, the voltage was reset to zero, and the neuron entered an absolute refractory period of ${\tau}_{r}=2$ ms. When driven by a constant input current $J$, a leaky integrateandfire neuron with absolute refractoriness fires at a rate $a=g(J)$ where $g$ is the gain function with value $g(J)=0$ for $J\le 1$ and
Our network was inhomogeneous in the sense that different neurons had different parameters as described below. The basic idea is that the ensemble of $N$ neurons, with different parameters, forms a rich set of basis functions in the ${N}_{c}$ or ${N}_{d}$ dimensional space of inputs or outputs, respectively. This is similar to tiling the space with radial basis functions, except that here we replace the radial functions by the gain functions of the LIF neurons (Equation (11)) each having different parameters (Eliasmith and Anderson, 2004). These parameters were chosen randomly once at the beginning of a simulation and kept fixed during the simulation.
For the command representation layer, we write the current $J$ into neuron $l$, in the case of a constant input $\overrightarrow{u}$, as
where $\nu}_{l}^{\mathrm{f}\mathrm{f}$ and $b}_{l}^{\mathrm{f}\mathrm{f}$ are neuronspecific gains and biases, and $\stackrel{~}{e}}_{l\beta}^{\mathrm{f}\mathrm{f}$ are ‘normalized’ encoding weights (cf. Equation (5)).
These random gains, biases and ‘normalized’ encoding weights must be chosen so that the command representation layer adequately represents the command input $\overrightarrow{u}$, whose norm is bounded in the interval $[0,{R}_{1}]$ (Table 1). First, we choose the ‘normalized’ encoding weight vectors on a hypersphere of radius $1/{R}_{1}$, so that the scalar product between the command vector and the vector of ‘normalized’ encoding weights, $\sum _{\beta}{\stackrel{~}{e}}_{l\beta}^{\mathrm{f}\mathrm{f}}{u}_{\beta}$, lies in the normalized range $[1,1]$. Second, the distribution of the gains sets the distribution of the firing rates in a target range. Third, we see from Equation (11) that the neuron starts to fire at the rheobase threshold $J=1$. The biases $b}_{l}^{\mathrm{f}\mathrm{f}$ randomly shift this rheobase threshold over an interval (see Figure 1—figure supplement 1). For the distributions used to set the fixed random gains and biases, see Table 1.
Analogously, for the recurrent network, we write the current into neuron $i$, for a constant ‘pseudoinput’ vector $\overrightarrow{\stackrel{~}{x}}$ being represented in the network, as
where ${\nu}_{i}$, ${b}_{i}$ are neuronspecific gains and biases, and ${\stackrel{~}{e}}_{i\alpha}$ are ‘normalized’ encoding weights. We call $\overrightarrow{\stackrel{~}{x}}$ a ‘pseudoinput’ for two reasons. First, the error encoding weights $k{e}_{i\alpha}$ are used to feed the error ${\u03f5}_{\alpha}=({x}_{\alpha}{\widehat{x}}_{\alpha})$ back to neuron $i$ in the network (cf. Equation (6)). However, ${\u03f5}_{\alpha}={x}_{\alpha}/(k+1)$, due to the feedback loop according to Equation (9). Thus, the ‘pseudoinput’ ${\stackrel{~}{x}}_{\alpha}=k{x}_{\alpha}/(k+1)$ has a similar range as $\overrightarrow{x}$, whose norm lies in the interval $[0,{R}_{2}]$ (see Table 1). Second, the neuron also gets feedforward and recurrent input. However, the feedforward and recurrent inputs get automatically adjusted during learning (starting from zero), so their absolute values do not matter for the initialization of parameters that we discuss here. Thus, we choose the ‘normalized’ encoding weight vectors on a hypersphere of radius $1/{R}_{2}$. For the distributions used to set the fixed random gains and biases, see Table 1.
Setting output decoding weights to form an autoencoder with respect to error encoding weights
Request a detailed protocolThe linear readout weights ${d}_{\alpha i}$ from the recurrently connected network were precomputed algorithmically so as to form an autoencoder with the error encoding weights ${e}_{i\alpha}$ (for $k=1$), while setting the feedforward and recurrent weights to zero (${w}_{l\beta}^{\mathrm{f}\mathrm{f}}=0$ and ${w}_{ij}=0$). To do this, we randomly selected $P$ error vectors ${\overrightarrow{\u03f5}}^{(p)}$, that we used as training samples for optimization, with sample index $p=1,\mathrm{\dots},P$, and having vector components ${\u03f5}_{\alpha}^{(p)}$, $\alpha =1,\mathrm{\dots},{N}_{d}$. Since the observable system is ${N}_{d}$dimensional, we chose the training samples randomly from within an ${N}_{d}$dimensional hypersphere of radius ${R}_{2}$. We applied each of the error vectors statically as input for the error feedback connections and calculated the activity
of neuron $i$ for error vector ${\overrightarrow{\u03f5}}^{(p)}$ using the static rate Equation (11). The decoders ${d}_{\alpha i}$ acting on these activities should yield back the encoded points thus forming an autoencoder. A squarederror loss function $\mathcal{L}$, with L2 regularization of the decoders,
setting $\lambda =P{\left(0.1\mathrm{max}\left(\left\{{a}_{i}^{(p)}\right\}\right)\right)}^{2}$ with number of samples $P=N$, was used for this linear regression (default in Nengo v2.4.0) (Eliasmith and Anderson, 2004; Stewart et al., 2009). Biologically plausible learning rules exist for autoencoders, either by training both encoding and decoding weights (Burbank, 2015), or by training decoding weights given random encoding weights (D'Souza et al., 2010; Urbanczik and Senn, 2014), but we simply calculated and set the decoding weights as if they had already been learned.
Compressive and expansive autoencoder
Request a detailed protocolClassical threelayer (inputhiddenoutputlayer) autoencoders 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 threelayer feedfoward network, our autoencoder forms a loop from the neurons in the recurrent network via readout weights to the output and from there via errorencoding weights to the input. Since the autoencoder is in the loop, we expect that it works both as a compressive one (from neurons in the recurrent network over the lowdimensional output back to the neurons) and as an expansive one (from the output through the neurons in the recurrent network back to the output).
Rather than constraining, as in Equation (15), the lowdimensional input ${\u03f5}_{\alpha}$ and roundtrip output ${\sum}_{i}{d}_{\alpha i}{a}_{i}\left(\overrightarrow{\u03f5}\right)$ to be equal for each component $\alpha $ (expansive autoencoder), we can alternatively enforce the high dimensional input ${I}_{j}$ (projection into neuron $j$ of lowdimensional input $\overrightarrow{\u03f5}$)
and roundtrip output ${I}_{j}^{\prime}\equiv {\sum}_{i,\alpha}{e}_{j\alpha}{d}_{\alpha i}{\stackrel{~}{g}}_{i}\left({I}_{i}\right)$, where ${\stackrel{~}{g}}_{i}\left({I}_{i}\right)\equiv {a}_{i}\left(\overrightarrow{\u03f5}\right)$, to be equal for each neuron $j$ in the recurrent network (compressive autoencoder) in order to optimize the decoding weights of the autoencoder. Thus, the squarederror loss for this compressive autoencoder becomes:
where in the approximation, we exploit that (i) the relative importance of the term involving ${\sum}_{p}^{P}{\sum}_{j}{\sum}_{\alpha ,\gamma ,\alpha \ne \gamma}{e}_{j\alpha}{e}_{j\gamma}$ tends to zero as $1/\sqrt{NP}$, since ${e}_{j\alpha}$ and ${e}_{j\gamma}$ are independent random variables; and (ii) ${\sum}_{j}{e}_{j\alpha}^{2}\approx c$ is independent of $\alpha $. Thus, the loss function of Equation (17) is approximately proportional to the squarederror loss function of Equation (15) (not considering the L2 regularization) used for the expansive autoencoder, showing that for an autoencoder embedded in a loop with fixed random encoding weights, the expansive and compressive descriptions are equivalent for those $N$dimensional inputs ${I}_{i}$ that lie in the ${N}_{d}$dimensional subspace spanned by $\{{e}_{i\alpha}\}$that is ${I}_{i}$ is of the form ${\sum}_{\alpha}{e}_{i\alpha}{\u03f5}_{\alpha}$ where ${\u03f5}_{\alpha}$ lies in a finite domain (hypersphere). We employed a large number $P=N$ of random low${N}_{d}$dimensional inputs when constraining the expansive autoencoder.
Command input
Request a detailed protocolThe command input vector $\overrightarrow{u}(t)$ to the network was ${N}_{c}$dimensional (${N}_{c}={N}_{d}$ for all systems except the arm) and timevarying. During the learning phase, input changed over two different time scales. The fast value of each command component was switched every 50 ms to a level ${u}_{\alpha}^{\prime}$ chosen uniformly between $({\zeta}_{1},{\zeta}_{1})$ and this number was added to a more slowly changing input variable ${\overline{u}}_{\alpha}$ (called ’pedestal’ in the main part of the paper) which changed with a period ${T}_{period}$. Here ${\overline{u}}_{\alpha}$ is the component of a vector of length ${\zeta}_{2}$ with a randomly chosen direction. The value of component $\alpha $ of the command is then ${u}_{\alpha}={\overline{u}}_{\alpha}+{u}_{\alpha}^{\prime}$. Parameter values for the network and input for each dynamical system are provided in Table 1. Further details are noted in the next subsection.
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.
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}$.
Linear system
Request a detailed protocolThe equations for a linear decaying oscillator system (Figure 2—figure supplement 2) are
For this linear dynamical system, we tested the learned network on a ramp of 2 s followed by a step to a constant nonzero 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 (Figure 2—figure supplement 2A–C).
van der Pol oscillator
Request a detailed protocolThe equations for the van der Pol oscillator system are
Each component of the pedestal input ${\overline{u}}_{\alpha}$ was scaled differently for the van der Pol oscillator as reported in Table 1.
Lorenz system
Request a detailed protocolThe equations for the chaotic Lorenz system (Lorenz, 1963) are
In our equations above, $Z$ of the original Lorenz equations (Lorenz, 1963) is represented by an output variable ${x}_{3}=Z28$ 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 ${\zeta}_{1}$, was provided to set off the system, after which the system followed autonomous dynamics.
Nonlinearly transformed input to linear system
Request a detailed protocolFor the above dynamical systems, the input adds linearly on the right hand sides of the differential equations. Our FOLLOW scheme also learned nonlinear feedforward inputs to a linear dynamical system, as demonstrated in Figure 2—figure supplement 4 and Figure 2—figure supplement 5. As the reference, we used the linear dynamical system above, but with its input transformed nonlinearly by ${g}_{\alpha}(\overrightarrow{u})=10({({u}_{\alpha}/0.1)}^{3}{u}_{\alpha}/0.4)$. Thus, the equations of the reference were:
The input to the network remained $\overrightarrow{u}$. Thus, effectively the feedforward weights had to learn the nonlinear transform $\overrightarrow{g}(\overrightarrow{u})$ while the recurrent weights learned the linear system.
Arm dynamics
Request a detailed protocolIn the example of learning arm dynamics, we used a twolink 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 $\overrightarrow{\theta}={({\theta}_{1},{\theta}_{2})}^{T}$ and the angular velocities $\overrightarrow{\omega}={({\omega}_{1},{\omega}_{2})}^{T}$, given input torques $\overrightarrow{\tau}={({\tau}_{1},{\tau}_{2})}^{T}$ are:
with
where ${m}_{i}$ is the mass, ${l}_{i}$ the length, ${s}_{i}$ the distance from the joint centre to the centre of the mass, and ${I}_{i}$ 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 $\overrightarrow{x}=[{\theta}_{1},{\theta}_{2},{\omega}_{1},{\omega}_{2}{]}^{T}$, but the effective torque $\overrightarrow{\tau}$ is obtained from the input torque $\overrightarrow{u}$ as follows.
To avoid any link from rotating full 360 degrees, we provide an effective torque ${\tau}_{\alpha}$ to the arm, by subtracting a term proportional to the input torque ${u}_{\alpha}$, if the angle crosses $\pm $90 degrees and ${u}_{\alpha}$ is in the same direction:
where $\stackrel{~}{\sigma}(\theta )$ increases linearly from 0 to 1 as $\theta $ goes from $\pi /2$ to $3\pi /4$:
The parameter values were taken from the human arm (Model 1) in section 3.1.1 of the PhD thesis of Li (Li, 2006) from the Todorov lab; namely $m}_{1}=1.4\phantom{\rule{thinmathspace}{0ex}}\mathrm{k}\mathrm{g$, $m}_{2}=1.1\phantom{\rule{thinmathspace}{0ex}}\mathrm{k}\mathrm{g$, $l}_{1}=0.3\phantom{\rule{thinmathspace}{0ex}}\mathrm{m$, $l}_{2}=0.33\phantom{\rule{thinmathspace}{0ex}}\mathrm{m$, $s}_{1}=0.11\phantom{\rule{thinmathspace}{0ex}}\mathrm{m$, $s}_{2}=0.16\phantom{\rule{thinmathspace}{0ex}}\mathrm{m$, $I}_{1}=0.025\phantom{\rule{thinmathspace}{0ex}}\mathrm{k}\mathrm{g}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{m}}^{2$, $I}_{2}=0.045\phantom{\rule{thinmathspace}{0ex}}\mathrm{k}\mathrm{g}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{m}}^{2$, and $b}_{11}={b}_{22}=0.05\phantom{\rule{thinmathspace}{0ex}}\mathrm{k}\mathrm{g}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{m}}^{2}/\mathrm{s$, $b}_{12}={b}_{21}=0.025\phantom{\rule{thinmathspace}{0ex}}\mathrm{k}\mathrm{g}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{m}}^{2}/\mathrm{s$. Acceleration due to gravity was set at $g=9.81\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}/{\mathrm{s}}^{2}$. For the arm, we did not filter the reference variables for calculating the error.
The input torque $\overrightarrow{u}(t)$ for learning the twolink arm was generated, not by switching the pulse and pedestal values sharply, every 50 ms and ${T}_{period}$ as for the others, but by linearly interpolating inbetween to avoid oscillations from sharp transitions.
The input torque $\overrightarrow{u}$ and the variables $\overrightarrow{\omega}$, $\overrightarrow{\theta}$ obtained on integrating the arm model above were scaled by $0.02$ (Nm)^{1}, $0.05$ (rad/s)^{1} and $1/2.5$ rad^{1} respectively, and then these dimensionless variables were used as the input and reference for the spiking network. Effectively, we scaled the input torques to cover onefifth of the representation radius ${R}_{2}$, the angular velocities onehalf, and the angles full, as each successive variable was the integral of the previous one.
Learning readout weights with recurrent weights fixed
Request a detailed protocolFor learning the readout weights after setting either the true or shuffled set of learned recurrent weights (Figure 2—figure supplement 3), we used a perceptron learning rule.
with learning rate ${\eta}_{r}=1e4$.
Derivation and proof of stability of the FOLLOW learning scheme
We derive the FOLLOW learning rule Equation (10), while simultaneously proving the stability of the scheme. We assume that: (1) the feedback $\{k{e}_{i\alpha}\}$ and readout weights $\{{d}_{\alpha 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 Appendix 1 for a discussion); (3) the state $\overrightarrow{x}$ of the dynamical system is observable; (4) the intrinsic time scales of the reference dynamics are much larger than the synaptic time scale and the time scale of the error feedback loop, and much smaller than the time scale of learning; (5) the feedforward and recurrent weights remain bounded; and (6) the input $\overrightarrow{u}$ and reference output $\overrightarrow{x}$ remain bounded.
The proof proceeds in three major steps: (1) using the autoencoder assumption to write the evolution equation of the lowdimensional output state variable in terms of the recurrent and feedforward weights; (2) showing that output follows the reference due to the error feedback loop; and (3) obtaining the evolution equation for the error and using it in the timederivative of a Lyapunov function $V$, to show that $\dot{V}\le 0$ for uniform stability, similar to proofs in adaptive control theory (Narendra and Annaswamy, 1989; Ioannou and Sun, 2012).
Role of network weights for lowdimensional output
Request a detailed protocolThe filtered lowdimensional output of the recurrent network is given by Equation (3) of Results and repeated here:
where ${d}_{\alpha j}$ are the readout weights. Since $\kappa $ is an exponential filter with time constant ${\tau}_{s}$, Equation (21) can also be written as
We convolve this equation with kernel $\kappa $, multiply by the error feedback weights, and sum over the output components $\alpha $
We would like to write Equation (23) in terms of the recurrent and feedforward weights in the network.
To do this, we exploit assumptions (1) and (4). Having shown the equivalence of the compressive and expansive descriptions of our autoencoder in the errorfeedback loop (Equations (15) and (17)), we formulate our nonlinear autoencoder as compressive: we start with a highdimensional set of inputs ${I}_{j}\equiv {J}_{j}{b}_{j}$ (where ${J}_{j}$ is the current into neuron $j$ with bias ${b}_{j}$, cf. Equations (5) and (6)); transform these inputs nonlinearly into filtered spike trains ${S}_{j}*\kappa $; decode these filtered spike trains into a lowdimensional representation $\overrightarrow{z}$ with components ${z}_{\alpha}={\sum}_{j}{d}_{\alpha j}({S}_{j}*\kappa )$; and increase the dimensionality back to the original one, via weights $k{e}_{i\alpha}$, to get inputs:
Using assumption (1) we expect that the final inputs ${I}_{i}^{\prime}$ are approximately $k$ times the initial inputs ${I}_{i}$:
This is valid only if the high$N$dimensional input ${I}_{i}$ lies in the low${N}_{d}$dimensional subspace spanned by $\{{e}_{i\alpha}\}$ (Equation (17)). We show that this requirement is fulfilled in the next major step of the proof (see text accompanying Equations (31)–(35)).
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 (25)) 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 ${I}_{i}^{\prime}$ on the left hand side of Equation (25) is varying only on the slow time scale. Therefore, we used firing rate equations to compute mean outputs given static input when precalculating the readout weights (earlier in Methods).
Inserting the approximate Equation (25) into Equation (23) we find
We replace ${I}_{i}\equiv {J}_{i}{b}_{i}$, using the current ${J}_{i}$ from Equation (6) for neuron $i$ of the recurrent network, to obtain
Thus, the change of the lowdimensional output ${\widehat{x}}_{\alpha}*\kappa $ depends on the network weights, which need to be learned. This finishes the first step of the proof.
Errorfeedback loop ensures that output follows reference
Request a detailed protocolBecause 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 $w}_{il}^{\mathrm{f}\mathrm{f}\ast$ and recurrent weights ${w}_{ij}^{*}$. We write an equation similar to Equation (27) for the output ${x}_{\alpha}^{*}$ of the target network:
where $({S}_{l}^{\mathrm{f}\mathrm{f}\ast}\ast \kappa )(t)$ and $({S}_{j}^{*}*\kappa )(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 ${S}_{j}^{*}$ gives the target output which we denote by ${\overrightarrow{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 input encoding 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 learned.
In view of potential generalization, we note that any nonlinear dynamical system is approximately realizable due to the expansion in a highdimensional nonlinear basis that is effectively performed by the recurrent network (see Appendix 1). Approximative weights (close to the ideal ones) could in principle also be calculated algorithmically (see Appendix 1). 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 ${\u03f5}_{\alpha}$ can be computed directly via a comparison of the true output $\overrightarrow{x}$ of the reference with the output $\overrightarrow{\widehat{x}}$ of the network: ${\u03f5}_{\alpha}={x}_{\alpha}{\widehat{x}}_{\alpha}.$ (In view of potential generalizations, we remark that the observable output need not be the state variables themselves, but could be a higherdimensional nonlinear function of the state variables, as shown for a general dynamical system in Appendix 1.)
As the second step of the proof, we now show that the error feedback loop enables the first network to follow the target network under assumptions (4  6). More precisely, we want to show that the readout and neural activities of the first network remain close to those of the target network at all times, that is ${\widehat{x}}_{\alpha}(t)\approx {x}_{\alpha}^{*}(t)$ for each component $\alpha $ and $({S}_{i}*\kappa )(t)\approx ({S}_{i}^{*}*\kappa )(t)$ 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, and (ii) the reference dynamics is slower than the synaptic and feedback loop time scales, so the reference output ${x}_{\alpha}$ can be assumed momentarily constant. Thus, we have a separation of time scales in Equation (27): for a given input (transmitted via the feedforward weights) and a given target value ${x}_{\alpha}^{*}$, the network dynamics settles on the fast time scale ${\tau}_{s}$ to a momentary fixed point ${\widehat{x}}^{\u2020}$ which we find by setting the derivative on the lefthand side of Equation (27) to zero:
We rewrite this equation in the form
We choose the feedback gain for the error much larger than 1 ($k\gg 1$), such that $k/(k+1)\approx 1$. We show below (in the text accompanying Equations (31)–(35)), that the factor in parentheses multiplying $1/(k+1)$ in the second term starts from zero and tends, with learning, towards $\sum _{\alpha}{e}_{i\alpha}({x}_{\alpha}^{\ast}\ast \kappa )$, which is the factor multiplying $k/(k+1)$ in the first term. Thus, the first term remains approximately $k$ times larger than the second during learning. To obtain $\hat{x}}_{\alpha}^{\u2020}\approx {x}_{\alpha}^{\ast$, we set $k\gg 1$.
To show that the momentary fixed point is stable at the fast synaptic time scale, we calculate the Jacobian $=[{}_{il}]$, for the dynamical system given by Equation (27). We introduce auxiliary variables ${y}_{i}\equiv \sum _{\alpha}{e}_{i\alpha}({\hat{x}}_{\alpha}\ast \kappa )$ to rewrite Equation (27) with the new variables in the form ${\dot{y}}_{i}={F}_{i}(\overrightarrow{y})$; and then we take derivative of its right hand side to obtain the elements of the Jacobian matrix at the fixed point $\sum _{\alpha}{e}_{i\alpha}{\hat{x}}_{\alpha}^{\u2020}$ (using ${\int}_{0}^{\mathrm{\infty}}\kappa (\tau )d\tau =1$):
where ${\delta}_{il}$ is the Kronecker delta function. We note that $\sum _{j}{w}_{ij}({S}_{j}\ast \kappa )$ is a spatially and temporally averaged measure of the population activity in the network with appropriate weighting factors ${w}_{ij}$. We assume that the population activity varies smoothly with input, which is equivalent to requiring that on the time scale ${\tau}_{s}$, the network fires asynchronously, i.e. there are no precisely timed population spikes. Then we can take the second term to be bounded, in absolute value, by say $B}_{1}/{\tau}_{s$. The Jacobian matrix $\mathcal{\mathcal{J}}$ is of the form $(k+1)\mathbb{I}/{\tau}_{s}+\mathrm{\Lambda}$, where $\mathbb{I}$ is the $N\times N$ identity matrix and $\mathrm{\Lambda}$ is a matrix with each element bounded in absolute value by $B}_{1}/{\tau}_{s$. If we set $k\gg N{B}_{1}$, then all eigenvalues of the Jacobian have negative real parts, applying the Gerschgorin circle theorem (the second term can perturb any eigenvalue from $(k+1)/{\tau}_{s}$ to within a circle of radius $N{B}_{1}/{\tau}_{s}$ at most), 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 system dynamics of the reference ${\overrightarrow{x}}^{*}$, the first network output follows the target network output at all times, ${\widehat{x}}_{\alpha}\approx {x}_{\alpha}^{*}$.
We now show that neurons are primarily driven by inputs close to those in the target network due to error feedback, and that these lie in the lowdimensional manifold spanned by $\{{e}_{i\alpha}\}$, as required for Equation (25). We compute the projected error using Equation (30):
and insert it into Equation (6) to obtain the current into a neuron in the recurrent network:
At the start of learning, if the feedforward and recurrent weights are small, then the neural input is dominated by the fedback error input that is the first term, making ${J}_{i}$ close to the ideal current
Thus, the neural input at the start of learning is of the form ${\sum}_{\alpha}{e}_{i\alpha}{x}_{\alpha}^{*}$ which lies in the lowdimensional subspace spanned by $\{{e}_{i\alpha}\}$ as required for Equation (25). Furthermore, over time, the feedforward and recurrent weights get modified so that their contribution tends towards ${\sum}_{\alpha}{e}_{i\alpha}({x}_{\alpha}^{*}*\kappa )$, such that the two terms of Equation (32) add to make ${J}_{i}$ even closer to the ideal current ${J}_{i}^{*}$ given by Equation (33). This is made clearer by considering the weight update rule Equation 10 as stochastic gradient descent on a loss function,
leading us to (for each recurrent weight ${w}_{ij}$, and similarly for $w}_{il}^{\mathrm{f}\mathrm{f}$):
which is identical to the FOLLOW learning rule for ${w}_{ij}$ in Equation (10) except for the timescale of filtering of the error current (see Discussion), and a factor involving $k$ that can be absorbed into the learning rate ${\eta}^{\prime}$. In the last step above, we used the projected error current from Equation (31) and the definition of ${I}_{i}^{\u03f5}$ in Equation (8). Thus, the feedforward and recurrent connections evolve to inject, after learning, the same ideal input within the lowdimensional manifold, as was provided by the error feedback during learning. Hence, the neural input remains in the lowdimensional manifold spanned by $\{{e}_{i\alpha}\}$ throughout learning, as required for Equation (25), making this major step and the previous one selfconsistent.
Since the driving neural currents are close to ideal throughout learning, the filtered spike trains of the recurrent neurons in the first network will also be approximately the same as those of the target network, so that $({S}_{i}*\kappa )(t)$ can be used instead of $({S}_{i}^{*}*\kappa )(t)$ in (Equation (28)). Moreover, the filtered spike trains $({S}_{l}^{\mathrm{f}\mathrm{f}}\ast \kappa )(t)$ of the command representation layer in the first network are the same as those in the target network, since they are driven by the same command input $\overrightarrow{u}$ and the command encoding weights are, by construction, the same for both networks. The similarity of the spike trains in the first and target networks will be used in the next major part of the proof.
Stability of learning via Lyapunov’s method
Request a detailed protocolWe now turn to the third step of the proof and consider the temporal evolution of the error ${\u03f5}_{\alpha}={x}_{\alpha}{\widehat{x}}_{\alpha}$. We exploit that the network dynamics is realized by the target network and insert Equations (27) and (28) so as to find
In the second line, we have replaced the reference output by the target network output; and in the third line we have used Equations (27) and (28), 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 ${\psi}_{ij}\equiv {w}_{ij}{w}_{ij}^{*}$ and $\varphi}_{il}\equiv {w}_{il}^{\mathrm{f}\mathrm{f}}{w}_{il}^{\mathrm{f}\mathrm{f}\ast$.
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:
where ${\stackrel{~}{\u03f5}}_{i}\equiv {\tau}_{s}{\sum}_{\alpha}{e}_{i\alpha}({\u03f5}_{\alpha}*\kappa )$ and ${\stackrel{~}{\eta}}_{1},{\stackrel{~}{\eta}}_{2}>0$ are positive constants. We use Lyapunov’s direct method to show the stability of learning. For this, we require the following properties for the Lyapunov function. (a) The Lyapunov function is positive semidefinite $V(\stackrel{~}{\u03f5},\psi ,\varphi )\ge 0$, with the equality to zero only at $(\stackrel{~}{\u03f5},\psi ,\varphi )=(0,0,0)$. (b) It has continuous firstorder partial derivatives. Furthermore, $V$ is (c) radially unbounded since
and (d) decrescent since
where ${(\stackrel{~}{\u03f5},\psi ,\varphi )}^{2}\equiv {\sum}_{i}{({\stackrel{~}{\u03f5}}_{i})}^{2}+{\sum}_{i,j}{({\psi}_{ij})}^{2}+{\sum}_{i,k}{({\varphi}_{il})}^{2}$ and $\mathrm{min}/\mathrm{max}$ take the minimum/maximum of their respective arguments.
Apart from the above conditions (a)(d), we need to show the key property $\dot{V}\le 0$ for uniform global stability (which implies that bounded orbits remain bounded, so the error remains bounded); or the stronger property $\dot{V}<0$ for asymptotic global stability (see for example [Narendra and Annaswamy, 1989; Ioannou and Sun, 2012]). Taking the time derivative of $V$, and replacing ${\dot{\stackrel{~}{\u03f5}}}_{i}$that is ${\tau}_{s}{\sum}_{\alpha}{e}_{i\alpha}({\dot{\u03f5}}_{\alpha}*\kappa )$ from (Equation (36)), we have:
If we enforce the first two terms above to be zero, we derive a learning rule
and then
requiring $k>1$, which is subsumed under $k\gg 1$ for the error feedback. Equation (39) with ${\eta}_{1}\equiv {\stackrel{~}{\eta}}_{1}{\tau}_{s}/k$ and ${\eta}_{2}\equiv {\stackrel{~}{\eta}}_{2}{\tau}_{s}/k$, and $\kappa $ replaced by a longer filtering kernel ${\kappa}^{\u03f5}$, is the learning rule used in the main text, Equation (10).
Thus, in the $(\stackrel{~}{\u03f5},\psi ,\varphi )$system given by Equations (36) and (39), we have proven the global uniform stability of the fixed point $(\stackrel{~}{\u03f5},\psi ,\varphi )=(0,0,0)$, which is effectively $(\u03f5,\psi ,\varphi )=(0,0,0)$, choosing ${\eta}_{1},{\eta}_{2}>0$ and $k\gg max(1,N{B}_{1})$, under assumptions (1  6), while simultaneously deriving the learning rule (Equation (39)).
This ends our proof. So far, we have shown that the system is Lyapunov stable, that is 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 subsection that the error $\overrightarrow{\u03f5}$ goes to zero asymptotically, so that after learning, even without error feedback, $\overrightarrow{\widehat{x}}$ reproduces the dynamics of $\overrightarrow{x}$.
A major caveat of this proof is that under assumption (2) the dynamics be 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 literature 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 nonrealizable 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 needed to stop weights from drifting. Various techniques to address such modelapproximation 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 Appendix 1.
To summarize, the FOLLOW learning rule (Equation (39)) on the feedforward or recurrent weights has two terms: (i) a filtered presynaptic firing trace $({S}_{l}^{\mathrm{f}\mathrm{f}}\ast \kappa )(t)$ or $({S}_{j}*\kappa )(t)$ that is available locally at each synapse; and (ii) a projected filtered error $\sum _{\alpha}k{e}_{i\alpha}({\u03f5}_{\alpha}\ast \kappa )(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 (6). Thus the learning rule can be classified as local. Moreover, it uses an error in the observable $\overrightarrow{x}$, not in its timederivative. While we have focused on spiking networks, the learning scheme can be easily used for nonlinear rate units by replacing the filtered spikes $({S}_{i}*\kappa )(t)$ by the output of the rate units $r(t)$. Our proof is valid for arbitrary dynamical transforms $\overrightarrow{h}(\overrightarrow{x},\overrightarrow{u})$ as long as they are realizable in a network. The proof shows uniform global stability using Lyapunov’s method.
Proof of error tending to zero asymptotically
Request a detailed protocolIn the above subsection, we showed uniform global stability using $\dot{V}=(k+1){\sum}_{i}{({\stackrel{~}{\u03f5}}_{i})}^{2}\le 0$, with $k\gg max(1,N{B}_{1})$ and ${\stackrel{~}{\u03f5}}_{i}\equiv {\tau}_{s}{\sum}_{\alpha}{e}_{j\alpha}({\u03f5}_{\alpha}*\kappa )$. 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 Barbălat’s lemma: if $f,\dot{f}\in {\mathcal{L}}_{\mathrm{\infty}}$ and $f\in {\mathcal{L}}_{p}$ for some $p\in [1,\mathrm{\infty})$, then $f(t)\to 0$ as $t\to \mathrm{\infty}$. Recall the definitions: function $f\in {\mathcal{L}}_{p}$ when ${x}_{p}\equiv {\left({\int}_{0}^{\mathrm{\infty}}{f(\tau )}^{p}d\tau \right)}^{1/p}$ exists (is finite); and similarly function $f\in {\mathcal{\mathcal{L}}}_{\mathrm{\infty}}$ when $x{}_{\mathrm{\infty}}\equiv \underset{t\ge 0}{sup}f(t)$ exists (is finite).
Since $V$ is positive semidefinite ($V\ge 0$) and is a nonincreasing function of time ($\dot{V}\le 0$), its ${lim}_{t\to \mathrm{\infty}}V={V}_{\mathrm{\infty}}$ exists and is finite. Using this, the following limit exists and is finite:
Since each term in the above sum ${\sum}_{i}$ is positive semidefinite, ${\int}_{0}^{\mathrm{\infty}}{({\stackrel{~}{\u03f5}}_{i}(\tau ))}^{2}d\tau $ also exists and is finite $\mathrm{\forall}i$, and thus ${\stackrel{~}{\u03f5}}_{i}\in {\mathcal{\mathcal{L}}}_{2}\phantom{\rule{thinmathspace}{0ex}}\mathrm{\forall}i$.
To show that ${\stackrel{~}{\u03f5}}_{i},{\dot{\stackrel{~}{\u03f5}}}_{i}\in {\mathcal{\mathcal{L}}}_{\mathrm{\infty}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{\forall}i$, consider Equation (36). We use assumption (6) that the input $\overrightarrow{u}(t)$ and the reference output $\overrightarrow{x}(t)$ are bounded. Since network output $\overrightarrow{\widehat{x}}$ is also bounded due to saturation of firing rates (as are the filtered spike trains), the error (each component) is bounded that is ${\stackrel{~}{\u03f5}}_{i}\in {\mathcal{\mathcal{L}}}_{\mathrm{\infty}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{\forall}i$. If we also bound the weights from diverging during learning (assumption (5)), then ${\psi}_{ij},{\varphi}_{il}\in {\mathcal{\mathcal{L}}}_{\mathrm{\infty}}\mathrm{\forall}i,j,l$. With these reasonable assumptions, all terms on the right hand side of the Equation (36) for ${\dot{\stackrel{~}{\u03f5}}}_{i}$ are bounded, hence ${\dot{\stackrel{~}{\u03f5}}}_{i}\in {\mathcal{\mathcal{L}}}_{\mathrm{\infty}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{\forall}i$.
Since ${\stackrel{~}{\u03f5}}_{i}\in {\mathcal{\mathcal{L}}}_{2}\phantom{\rule{thinmathspace}{0ex}}\mathrm{\forall}i$ and $\stackrel{~}{{\u03f5}_{i}},{\dot{\stackrel{~}{\u03f5}}}_{i}\in {\mathcal{\mathcal{L}}}_{\mathrm{\infty}}\phantom{\rule{thinmathspace}{0ex}}\mathrm{\forall}i$, invoking Barbălat’s lemma as above, we have ${\stackrel{~}{\u03f5}}_{i}\to 0\phantom{\rule{thinmathspace}{0ex}}\mathrm{\forall}i$ as $t\to \mathrm{\infty}$. 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 errorfeedback, and spiking shot noise (cf. Figure 5).
We do not further require the convergence of parameters to ideal ones for our purpose, since the error tending to zero, that is 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 (weights) 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 (Figure 5), but is not pursued here.
Appendix 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 $\overrightarrow{v}(\overrightarrow{u})$ corresponding to the $\overrightarrow{u}$ encoded in the command representation layer, from the spike trains ${S}_{l}^{\mathrm{f}\mathrm{f}}(t)$ of the neurons, by synaptically filtering and linearly weighting the trains with decoding weights ${d}_{\alpha l}^{(\overrightarrow{v})}$:
where $*$ denotes convolution $({S}_{l}^{\mathrm{f}\mathrm{f}}\ast \kappa )(t)\equiv {\int}_{\mathrm{\infty}}^{t}{S}_{l}^{\mathrm{f}\mathrm{f}}({t}^{\mathrm{\prime}})\kappa (t{t}^{\mathrm{\prime}})d{t}^{\mathrm{\prime}}={\int}_{0}^{\mathrm{\infty}}{S}_{l}^{\mathrm{f}\mathrm{f}}(t{t}^{\mathrm{\prime}})\kappa ({t}^{\mathrm{\prime}})d{t}^{\mathrm{\prime}}$, and $\kappa (t)\equiv \mathrm{exp}(t/{\tau}_{s})/{\tau}_{s}$ is a normalized filtering kernel.
We can obtain the decoders ${d}_{\alpha i}^{(\overrightarrow{v})}$ by minimizing the loss function
with respect to the decoders. The average ${\u27e8\cdot \u27e9}_{\overrightarrow{u}}$ over $\overrightarrow{u}$ guarantees that the same constant decoders are used over the whole range of constant inputs $\overrightarrow{u}$. The time average ${\u27e8\cdot \u27e9}_{t}$ denotes an analytic rate computed for each constant input for a LIF neuron. Linear regression with a finite set of constant inputs $\overrightarrow{u}$ was used to obtain the decoders (see Methods). With these decoders, if the input $\overrightarrow{u}$ varies slowly compared to the synaptic time constant ${\tau}_{s}$, we have ${\hat{v}}_{\alpha}=\sum _{l}{d}_{\alpha l}^{(\overrightarrow{v})}({S}_{l}^{\mathrm{f}\mathrm{f}}\ast \kappa )(t)\approx {v}_{\alpha}(\overrightarrow{u})$.
Any function of the input $\overrightarrow{v}(\overrightarrow{u})$ can be approximated with appropriate linear decoding weights ${d}_{\alpha l}^{(\overrightarrow{v})}$ from the highdimensional basis of nonlinear tuning curves of heterogeneous neurons with different biases, encoding weights and gains, schematized in Figure 1—figure supplement 1. 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 nonlinear basis functions (Funahashi, 1989; Hornik et al., 1989; Girosi and Poggio, 1990) successful usage in spiking neural network models by various groups (Seung et al., 2000; Eliasmith and Anderson, 2004; Eliasmith, 2005), and biological plausibility (Poggio, 1990; Burnod et al., 1992; Pouget and Sejnowski, 1997).
Here, the neurons that are active at any given time operate in the mean driven regime, that is 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 (41) with stationary input are good approximations even in the timedependent case, as long as the input varies on a time scale slower than the synaptic time constant.
Online learning based on a loss function and its shortcomings
Suppose that a dynamical system given by
is to be mimicked by our spiking network implementing a different dynamical system with an extra error feedback term as in Equation (27). This can be interpreted as:
Comparing with the reference Equation (42), after learning we want that ${\stackrel{\u02d8}{f}}_{\alpha}(\overrightarrow{\widehat{x}},\overrightarrow{u})+{\stackrel{\u02d8}{g}}_{\alpha}(\overrightarrow{u})$ should approximate ${\tau}_{s}{f}_{\alpha}(\overrightarrow{\widehat{x}})+{\widehat{x}}_{\alpha}+{\tau}_{s}{g}_{\alpha}(\overrightarrow{u})$. One way to achieve this (Eliasmith and Anderson, 2004) is to ensure that ${\stackrel{\u02d8}{f}}_{\alpha}(\overrightarrow{\widehat{x}},\overrightarrow{u})$ and ${\stackrel{\u02d8}{g}}_{\alpha}(\overrightarrow{u})$ approximate ${\stackrel{~}{f}}_{\alpha}(\overrightarrow{\widehat{x}})\equiv {\tau}_{s}{f}_{\alpha}(\overrightarrow{\widehat{x}})+{\widehat{x}}_{\alpha}$ and ${\stackrel{~}{g}}_{\alpha}(\overrightarrow{u})\equiv {\tau}_{s}{g}_{\alpha}(\overrightarrow{u})$ respectively, as used in the loss functions below. In our simulations, we usually start with zero feedforward and recurrent weights, so that initially $\stackrel{\u02d8}{f}(\overrightarrow{\widehat{x}},\overrightarrow{u})=0={\stackrel{\u02d8}{g}}_{\alpha}(\overrightarrow{u})$.
Assuming that the time scales of dynamics are slower than synaptic time scale ${\tau}_{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):
Using these loss functions, we can precalculate the weights required for any dynamical system numerically, similarly to the calculation of decoders in the subsection above.
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 (45) is
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 $\eta $. With requisite dynamics slower than synaptic ${\tau}_{s}$, and with large enough number of neurons, we have approximated ${\sum}_{i}{w}_{ji}{\u27e8{S}_{i}*\kappa \u27e9}_{t}(t)\approx {\sum}_{i}{w}_{ji}({S}_{i}*\kappa )(t)$. The third line defines an error in the projected $\overrightarrow{\stackrel{~}{f}}(\overrightarrow{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
where ${\u03f5}_{j}^{(\stackrel{~}{f})}\equiv \left({\sum}_{\alpha}{e}_{j\alpha}{\stackrel{~}{f}}_{\alpha}(\overrightarrow{x}){\sum}_{i}{w}_{ji}({S}_{i}*\kappa )(t)\right)$. This learning rule is the product of a multidimensional error ${\u03f5}_{j}^{(\stackrel{~}{f})}$ and the filtered presynaptic spike train $({S}_{i}*\kappa )(t)$. However, this error in the unobservable $\overrightarrow{\stackrel{~}{f}}$ is not available to the postsynaptic neuron, making the learning rule nonlocal. A similar issue arises in the feedforward case.
In mimicking a dynamical system, we want only the observable output of the dynamical system, that is $\overrightarrow{x}$ to be used in a supervisory signal, not a term involving the unknown $\overrightarrow{f}(\overrightarrow{x})$that appears in the derivative $\dot{\overrightarrow{x}}$. Even if this derivative is computed from the observable $\overrightarrow{x}$, it will be noisy. Furthermore, this derivative cannot be obtained by differentiating the observable versus time, if the observable is not directly the state variable, but an unknown nonlinear function of it, which however our FOLLOW learning can handle (see next subsection). Thus, this online rule, if using just the observable error, can learn only an integrator for which $f(x)\sim 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.
General dynamical system and transformed observable
General dynamical systems of the form
can be learned with the same network configuration (Figure 1B) used for systems of the form Equation 1. Here, the state variable is $\overrightarrow{x}$, but the observable which serves as the reference to the network is $\overrightarrow{y}$. The transformation equation of the observable (second equation) can be absorbed into the first equation as below.
Consider the transformation equation for the observable. The dimensionality of the relevant variables: (1) the state variables (say joint angles and velocities) $\overrightarrow{x}$; (2) the observables represented in the brain (say sensory representations of the joint angles and velocities) $\overrightarrow{y}$; and (3) the control input (motor command) $\overrightarrow{u}$, can be different from each other, but must be small compared to the number of neurons. Furthermore, we require the observable $\overrightarrow{y}$ to not lose information compared to $\overrightarrow{x}$, that is $\overrightarrow{K}$ must be invertible, so $\overrightarrow{y}$ will have at least the same dimension as $\overrightarrow{x}$.
The time evolution of the observable is
The last step follows since function $\overrightarrow{K}$ is invertible, so that $\overrightarrow{x}={\overrightarrow{K}}^{1}(\overrightarrow{y})$. So we essentially need to learn ${\dot{y}}_{\beta}={p}_{\beta}(\overrightarrow{y},\overrightarrow{u})$.
Having solved the observable transformation issue, we use $\overrightarrow{x}$ now for our observable instead of $\overrightarrow{y}$, consistent with the main text. The dynamical system to be learned is now ${\dot{x}}_{\beta}={p}_{\beta}(\overrightarrow{x},\overrightarrow{u})$. Since our learning network effectively evolves as Equation (43), it can approximate ${p}_{\beta}(\overrightarrow{x},\overrightarrow{u})$. Thus our network can learn general dynamical systems with observable transformations.
Approximation error causes drift in weights
A frozen noise term $\xi (\overrightarrow{x}(t))$ due to the approximate decoding from nonlinear tuning curves of neurons, by the feedforward weights, recurrent weights and output decoders, will appear additionally in Equation (36). If this frozen noise has a nonzero mean over time as $\overrightarrow{x}(t)$ varies, leading to a nonzero mean error, then it causes a drift in the weights due to the errorbased learning rules in Equation (10), 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 updating 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.
References

Building functional networks of spiking model neuronsNature Neuroscience 19:350–355.https://doi.org/10.1038/nn.4241

Learning longterm dependencies with gradient descent is difficultIEEE Transactions on Neural Networks 5:157–166.https://doi.org/10.1109/72.279181

Deep networks for motor control functionsFrontiers in Computational Neuroscience 9:32.https://doi.org/10.3389/fncom.2015.00032

BookEnforcing balance allows local supervised learning in spiking recurrent networksIn: Cortes C, Lawrence N. D, Lee D. D, Sugiyama M, Garnett R, Garnett R, editors. Advances in Neural Information Processing Systems, 28. Curran Associates, Inc. pp. 982–990.

Mirrored STDP implements autoencoder learning in a network of spiking neuronsPLoS Computational Biology 11:e1004566.https://doi.org/10.1371/journal.pcbi.1004566

Visuomotor transformations underlying arm movements toward visual targets: a neural network model of cerebral cortical operationsJournal of Neuroscience 12:1435–1453.

Modeling of continuous time dynamical systems with input by recurrent neural networksIEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications 47:575–578.https://doi.org/10.1109/81.841860

Every good regulator of a system must be a model of that systemInternational Journal of Systems Science 1:89–97.https://doi.org/10.1080/00207727008920220

A learningbased approach to artificial sensory feedback leads to optimal integrationNature Neuroscience 18:138–144.https://doi.org/10.1038/nn.3883

Widespread access to predictive models in the motor system: a short reviewJournal of Neural Engineering 2:S313–S319.https://doi.org/10.1088/17412560/2/3/S11

Efficient codes and balanced networksNature Neuroscience 19:375–382.https://doi.org/10.1038/nn.4243

A spiking neural model of adaptive arm controlProceedings of the Royal Society B: Biological Sciences 283:20162134.https://doi.org/10.1098/rspb.2016.2134

BookNeural Engineering: Computation, Representation, and Dynamics in Neurobiological SystemsMIT Press.

A unified approach to building and controlling spiking attractor networksNeural Computation 17:1276–1314.https://doi.org/10.1162/0899766053630332

Hierarchical models in the brainPLoS Computational Biology 4:e1000211.https://doi.org/10.1371/journal.pcbi.1000211

BookNeuronal Dynamics: From Single Neurons to Networks and Models of Cognition (1st edition)Cambridge, United Kingdom: Cambridge University Press.https://doi.org/10.1017/CBO9781107447615

Networks and the best approximation propertyBiological Cybernetics 63:169–176.https://doi.org/10.1007/BF00195855

The tempotron: a neuron that learns spike timingbased decisionsNature Neuroscience 9:420–428.https://doi.org/10.1038/nn1643

To spike, or when to spike?Current Opinion in Neurobiology 25:134–139.https://doi.org/10.1016/j.conb.2014.01.004

Gradient Flow in Recurrent Nets: The Difficulty of Learning LongTerm DependenciesGradient Flow in Recurrent Nets: The Difficulty of Learning LongTerm Dependencies.

BookAdaptive Control TutorialPhiladelphia, PA: SIAM, Society for Industrial and Applied Mathematics.https://doi.org/10.1137/1.9780898718652

A robust direct adaptive controllerIEEE Transactions on Automatic Control 31:1033–1043.https://doi.org/10.1109/TAC.1986.1104168

ReportThe ”Echo State” Approach to Analysing and Training Recurrent Neural NetworksTechnical report.

A Tutorial on Training Recurrent Neural Networks, Covering BPPT, RTRL, EKF and the "Echo StateNetwork" ApproachA Tutorial on Training Recurrent Neural Networks, Covering BPPT, RTRL, EKF and the "Echo StateNetwork" Approach.

Movement generation with circuits of spiking neuronsNeural Computation 17:1715–1738.https://doi.org/10.1162/0899766054026684

Neural basis of sensorimotor learning: modifying internal modelsCurrent Opinion in Neurobiology 18:573–581.https://doi.org/10.1016/j.conb.2008.11.003

Input prediction and autonomous movement analysis in recurrent circuits of spiking neuronsReviews in the Neurosciences 14:5–19.https://doi.org/10.1515/REVNEURO.2003.14.12.5

On the adaptive control of robot manipulatorsThe International Journal of Robotics Research 6:49–59.https://doi.org/10.1177/027836498700600303

BookOptimal Control for Biological Movement SystemsSan Diego: University of California.

Random synaptic feedback weights support error backpropagation for deep learningNature Communications 7:13276.https://doi.org/10.1038/ncomms13276

Deterministic nonperiodic flowJournal of the Atmospheric Sciences 20:130–141.https://doi.org/10.1175/15200469(1963)020<0130:DNF>2.0.CO;2

On the computational power of circuits of spiking neuronsJournal of Computer and System Sciences 69:593–616.https://doi.org/10.1016/j.jcss.2004.04.001

Explaining facial imitation: a theoretical modelEarly Development and Parenting 6:179–192.https://doi.org/10.1002/(SICI)10990917(199709/12)6:3/4<179::AIDEDP157>3.0.CO;2R

Span: spike pattern association neuron for learning spatiotemporal spike patternsInternational Journal of Neural Systems 22:1250012.https://doi.org/10.1142/S0129065712500128

Global stability of parameteradaptive control systemsIEEE Transactions on Automatic Control 25:433–439.https://doi.org/10.1109/TAC.1980.1102364

Stable adaptive controller design, part II: Proof of stabilityIEEE Transactions on Automatic Control 25:440–448.https://doi.org/10.1109/TAC.1980.1102362

Solving the problem of negative synaptic weights in cortical modelsNeural Computation 20:1473–1494.https://doi.org/10.1162/neco.2008.0706295

Gradient calculations for dynamic recurrent neural networks: a surveyIEEE Transactions on Neural Networks 6:1212–1228.https://doi.org/10.1109/72.410363

A theory of how the brain might workCold Spring Harbor Symposia on Quantitative Biology 55:899–910.https://doi.org/10.1101/SQB.1990.055.01.084

Spatial transformations in the parietal cortex using basis functionsJournal of Cognitive Neuroscience 9:222–237.https://doi.org/10.1162/jocn.1997.9.2.222

Computational approaches to sensorimotor transformationsNature Neuroscience 3:1192–1198.https://doi.org/10.1038/81469

Principles of neurodynamics. perceptrons and the theory of brain mechanismsTechnical report.

Parallel Distributed Processing: Explorations in the Microstructure of CognitionLearning Internal Representations by Error Propagation, Parallel Distributed Processing: Explorations in the Microstructure of Cognition, 1, Cambridge, MA, USA, MIT Press.

Gaussian networks for direct adaptive controlIEEE Transactions on Neural Networks 3:837–863.https://doi.org/10.1109/72.165588

The roles of vision and proprioception in the planning of reaching movementsAdvances in Experimental Medicine and Biology 629:317–335.https://doi.org/10.1007/9780387770642_16

BookAdaptive Control: Stability, Convergence, and RobustnessUpper Saddle River, NJ, USA: PrenticeHall, Inc.

Adaptive representation of dynamics during learning of a motor taskJournal of Neuroscience 14:3208–3224.

Adaptive sliding controller synthesis for nonlinear systemsInternational Journal of Control 43:1631–1651.https://doi.org/10.1080/00207178608933564

Python scripting in the Nengo simulatorFrontiers in Neuroinformatics 3:.https://doi.org/10.3389/neuro.11.007.2009

Generalization in reinforcement learning: Successful examples using sparse coarse codingAdvances in Neural Information Processing Systems 8:138–1044.

Learning universal computations with spikesPLoS Computational Biology 12:e1004895.https://doi.org/10.1371/journal.pcbi.1004895

Uncertainty of feedback and state estimation determines the speed of motor adaptationFrontiers in Computational Neuroscience 4:11.https://doi.org/10.3389/fncom.2010.00011

Computational principles of movement neuroscienceNature Neuroscience 3:1212–1217.https://doi.org/10.1038/81497

Forward models for physiological motor controlNeural Networks 9:1265–1279.https://doi.org/10.1016/S08936080(96)000354

Can proprioceptive training improve motor learning?Journal of Neurophysiology 108:3313–3321.https://doi.org/10.1152/jn.00122.2012

Fast adaptation of the internal model of gravity for manual interceptions: evidence for eventdependent learningJournal of Neurophysiology 93:1055–1068.https://doi.org/10.1152/jn.00833.2004

Visuomotor coordination and internal models for object interceptionExperimental Brain Research 192:571–604.https://doi.org/10.1007/s0022100816913

Stable adaptive control with recurrent neural networks for square MIMO nonlinear systemsEngineering Applications of Artificial Intelligence 22:702–717.https://doi.org/10.1016/j.engappai.2008.12.005
Decision letter

Peter LathamReviewing Editor; University College London, United Kingdom
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Predicting nonlinear dynamics by stable local learning in a recurrent spiking neural network" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors and the evaluation has been overseen by Timothy Behrens as the Senior Editor.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
This paper proposes a learning scheme that allows a recurrent network of spiking neurons to predict the dynamics of a nonlinear system by training feedforward and recurrent connections. The learning is achieved using a local learning rule. The error signal is fed back by fixed, random connections, with the feedback loop forming an autoencoder. After learning, which is done via motor "babbling," the network is able to reproduce the output in an openloop setting without feedback. The authors showcase their training algorithm in a number of linear and nonlinear applications, including the dynamics of a twolink arm. In addition, the authors prove analytically under mild assumptions that learning converges if the nonlinear dynamics can be realized by the network. The paper is mainly clearly written, and has the potential to be of significant interest to the community.
Essential revisions:
1) The firing rates are far higher than is biologically plausible. It is important that this is fixable, at least in principle. If not, then the scheme has nothing to do with the brain. We hesitate to ask for additional simulations, but we think you need to show that the network can work when firing rates are in a more biologically plausible range, as in less than about 20 Hz on average (still high, but closer to reality).
2) We did not follow a key step. The sixth equation in subsection “Network architecture and parameters” says:
sum_{beta, i} e_{j,beta} d_{beta,i} a_i(sum_alpha e_{i,alpha} epsilon_alpha)
\approx sum_alpha e_{j,alpha} epsilon_alpha
It seems that you effectively rewrote this as
sum_{beta, i} e_{j,beta} d_{beta,i} a_i(I_i) \approx I_i
(see Eq. (13)). Strictly speaking, this holds only if the current (whose ith component is I_i) lies in the low dimensional submanifold spanned by the e_{i,alpha}. However, because of the recurrent dynamic, the network is likely be chaotic, and so the current is likely to live on a highdimensional manifold.
This seems to indicate a problem with the derivation. Or are we missing something? In any case, this needs to be explained much more clearly.
3) It would be good to have a clear, succinct description of what new feature allowed you to train spiking networks with a local, biologically realistic learning rule, something that others have failed to do. Perhaps it's the autoencoder? (See next comment.)
4) The autoencoder is critical for the operation of the network. Can you provide any intuition about why it's there, and what it does? This may not be possible, but it would greatly strengthen the manuscript if you could do it.
5) A potentially important limitation of this paper is that it does not take into account feedback delays during training. This would be a significant limitation if the purpose of this work is to propose how forward models can be learned in the brain: one of the main reasons to have forward models in the first place is being able to deal with feedback delays (see e.g. Wolpert, Ghahramani & Jordan, 1995). Learning that ignores feedback delays is therefore not fully biologically realistic. If this limitation is indeed correct, it needs to be mentioned more prominently in the Discussion section (we only found it buried in the Methods section); if the network can in fact handle feedback delays during training (that are on a similar scale as the target dynamics), it should be demonstrated, as it would make the impact of the paper significantly stronger.
6) An alternative approach to training recurrent networks is to initialize the recurrent weights in some appropriate way and to just train the readout weights during learning (reservoir computing). How much more computationally powerful is the proposed method? A simple way to test this would be to compare the performance of the trained network to that of a network where the recurrent weights have been shuffled after training, the feedback weights set to zero, and the output weights retrained by linear regression. If such a network had comparable computational performance, it would suggest that only the correct recurrent weight distribution is required to learn the task.
7) Does learning still work if the dynamical system has small amounts of stochasticity, or if the error signal is noisy?
8) It's critical to have experimental predictions. Or at least a summary of the implications of this work for what's going on in the brain.
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Thank you for submitting your revised article "Predicting nonlinear dynamics by stable local learning in a recurrent spiking neural network" for consideration by eLife. The evaluation has been overseen by a Peter Latham as the Reviewing Editor and Timothy Behrens as the Senior Editor.
I want to congratulate you on a very nice paper, and thank you for your very thorough responses; you addressed almost all the points. The only problem is that in parts the manuscript is still not very clear (at least not to me). Please take a very hard look at the Methods section "Network architecture and parameters", and make sure it's crystal clear.
https://doi.org/10.7554/eLife.28295.022Author response
Summary:
This paper proposes a learning scheme that allows a recurrent network of spiking neurons to predict the dynamics of a nonlinear system by training feedforward and recurrent connections. The learning is achieved using a local learning rule. The error signal is fed back by fixed, random connections, with the feedback loop forming an autoencoder. After learning, which is done via motor "babbling," the network is able to reproduce the output in an openloop setting without feedback. The authors showcase their training algorithm in a number of linear and nonlinear applications, including the dynamics of a twolink arm. In addition, the authors prove analytically under mild assumptions that learning converges if the nonlinear dynamics can be realized by the network. The paper is mainly clearly written, and has the potential to be of significant interest to the community.
We are sincerely thankful to the reviewers and editors for their insightful comments. We have incorporated all the essential revisions and minor points in our revised manuscript, as detailed in the individual responses below. We hope we have addressed them all to your satisfaction.
Essential revisions:
1) The firing rates are far higher than is biologically plausible. It is important that this is fixable, at least in principle. If not, then the scheme has nothing to do with the brain. We hesitate to ask for additional simulations, but we think you need to show that the network can work when firing rates are in a more biologically plausible range, as in less than about 20 Hz on average (still high, but closer to reality).
We thank the reviewers for exhorting us to reduce the firing rates to a more biologically plausible range.
We have replaced Figure 2 with newer simulations run with modified neuronal parameters that yield 13 Hz mean firing rate (averaged over 16s) now compared to 37 Hz before. Third paragraph of subsection “Nonlinear oscillator” now reads:
“The distributions of firing rates averaged over a 0.25 s period with fairly constant output, and over a 16 s period with timevarying output, were longtailed, with the mean across neurons maintained at approximately 1213 Hz (Figure 2E, F).”
We have also modified the Methods section to include these different neuronal parameters in the second para of subsection “Network architecture and parameters”:
“For the low firing rate simulation in Figure 2, the gains $\nu_i$ were all set to 2, with the biases chosen uniformly from (−2, 2). Choosing a fixed gain for all neurons led to lower variance in firing rates compared to random gains.”
See also Figure 1—figure supplement 1 for representative neuronal tuning curves.
2) We did not follow a key step. The sixth equation in subsection “network architecture and parameters” says:
sum_{beta, i} e_{j,beta} d_{beta,i} a_i(sum_alpha e_{i,alpha} epsilon_alpha)
\approx sum_alpha e_{j,alpha} epsilon_alpha
It seems that you effectively rewrote this as
sum_{beta, i} e_{j,beta} d_{beta,i} a_i(I_i) \approx I_i
(see Eq. (13)). Strictly speaking, this holds only if the current (whose ith component is I_i) lies in the low dimensional submanifold spanned by the e_{i,alpha}. However, because of the recurrent dynamic, the network is likely be chaotic, and so the current is likely to live on a highdimensional manifold.
This seems to indicate a problem with the derivation. Or are we missing something? In any case, this needs to be explained much more clearly.
We profusely thank the reviewers for pointing out this issue.
The neural input actually remains always in the requisite lowdimensional manifold. We have now clarified and derived it in detail in the text as below.
Below equation (17) for the compressive autoencoder we write:
“Thus, the loss function of equation (17) is approximately proportional to the squarederror loss function of equation (15) used for the expansive autoencoder, showing that for an autoencoder embedded in a loop with fixed random encoding weights, the expansive and compressive descriptions are equivalent for those $N$dimensional inputs $I_i$ that lie in the $N_d$dimensional subspace spanned by $\{e_{i\alpha}\}$ i.e. $I_i$ is of the form $\sum_\alpha e_{i\alpha} \epsilon_\alpha$ where $\epsilon_\alpha$ lies in a finite domain (hypersphere).”
Below equation (25) [abovereferred previous equation (13)] in the first major step of the proof, we have added:
“This is valid only if the high$N$dimensional input $I_i$ lies in the low$N_d$dimensional subspace spanned by $\{e_{i\alpha}\}$ (equation (17)). We show that this requirement is fulfilled in the next major step of the proof (see text accompanying equations (31)(35)).”
The text in subsection “Errorfeedback loop ensures that output follows reference” around new equations (31)(35) in the second major step of the proof confirms this requirement. The text is too large to excerpt here.
In brief, the derivation shows that the fedback error current dominates over the feedforward and recurrent contributions, beginning with small feedforward and recurrent weights, at the start of learning. The fedback error input lies in the low dimensional manifold spanned by the error encoding weights, driving the neural activity close to ideal (also in a lowdimensional manifold) and this in turn ensures that the feedforward and recurrent weights evolve to contribute terms essentially in the lowdimensional manifold. Thus, at all times in the learning and testing, the neural currents lie in a low dimensional manifold.
3) It would be good to have a clear, succinct description of what new feature allowed you to train spiking networks with a local, biologically realistic learning rule, something that others have failed to do. Perhaps it's the autoencoder? (See next comment.)
We thank the reviewers for asking us to clarify the new features of our learning scheme.
We address both points 3 and 4 with this response. The key new features are the use of the errorfeedback loop along with the autoencoder.
We have introduced a subsection early in the Results section titled: “Negativefeedback of error via autoencoder enables local learning”.
This incorporates old and new text, explaining what new features enable learning and how.
In brief, the error feedback using the autoencoder of gain k, serves two purposes: (1) it makes the projected error available as a local current in every neuron in the recurrent network, and (2) it drives each neuron's activity in the recurrent network to be close to ideal. Thus, since the activities are already ideal, only the weights are responsible for the output error, and the weights can be trained with a simple perceptronlike learning rule, overcoming the creditassignment problem. In FORCE, for example, neural activities are also forced close to ideal, but by very fast weights swings initially, at a time scale faster than the reference dynamics (biologically implausible), and without the error being projected locally in each neuron (nonlocal).
In particular please see the paragraphs (much easier to read in the typeset pdf, hence not excerpting fully here):
Subsection “Negative error feedback via autoencoder enables local learning”:
"The combination of the autoencoder and the error feedback implies that the output stays close to the reference, as explained now..."
In the Methods section, we show that not just the lowdimensional output $\vec{\hat{x}}$, but also the spike trains $S_i(t)$, for $i=1,\ldots,N$, are entrained by the error feedback to be close to the ideal ones required to generate $\vec{x}$.
"During learning, the error feedback via the autoencoder in a loop serves two roles: (i) to make the error current available in each neuron, projected correctly, for a local synaptic plasticity rule, and (ii) to drive the spike trains to the target ones for producing the reference output. In other learning schemes for recurrent neural networks, where neural activities are not constrained by error feedback, it is difficult to assign credit or blame for the momentarily observed error, because neural activities from the past affect the present output in a recurrent network. In the FOLLOW scheme, the spike trains are constrained to closely follow the ideal time course throughout learning, so that the present error can be attributed directly to the weights, enabling us to change the weights with a simple perceptronlike learning rule \citep{rosenblatt_principles_1961} as in equation (10), bypassing the credit assignment problem. In the perceptron rule, the weight change $\Δ w \sim (\text{pre})\cdot \δ$ is proportional to the presynaptic input $(\text{pre})$ and the error $\δ$. In the FOLLOW learning rule of equation (10), we can identify $(S_i*\kappa)$ with $(\text{pre})$ and $(I_i^\epsilon * \kappa^\epsilon)$ with $\δ$. In Methods, we derive the learning rule of equation (10) in a principled way from a stability criterion."
We further mention how this is different from FORCE and some other similar methods (briefly here after introducing the rule and above intuition in the results, see also the detailed comparisons in the Discussion section):
"FORCE learning [Sussillo and Abbott, 2009; 2012; DePasquale et al., 2016; Thalmeier et al., 2016; Nicola and Clopath, 2016] also clamps the output and neural activities to be close to ideal during learning, by using weight changes that are faster than the time scale of the dynamics. In our FOLLOW scheme, clamping is achieved via negative error feedback using the autoencoder, which allows weight changes to be slow and makes the error current available locally in the postsynaptic neuron. Other methods used feedback based on adaptive control for learning in recurrent networks of spiking neurons, but were limited to linear systems [MacNeil and Eliasmith, 2011; Bourdoukan and Deneve, 2015], whereas the FOLLOW scheme was derived for nonlinear systems (see Methods section). Our learning rule of equation (10) uses an error $\epsilon_\alpha \equiv x_\alpha – \hat{x}_\alpha$ in the observable state, rather than an error involving the derivative $dx_\alpha/dt$ in equation (1), as in other schemes (see Appendix 1) [Eliasmith, 2005; MacNeil and Eliasmith, 2011]. The reader is referred to the Discussion section for detailed further comparisons. The FOLLOW learning rule is local since all quantities needed on the righthandside of equation (10) could be available at the location of the synapse in the postsynaptic neuron. For a potential implementation and prediction for errorbased synaptic plasticity, and for a critical evaluation of the notion of ‘local rule’, we refer to the Discussion section."
4) The autoencoder is critical for the operation of the network. Can you provide any intuition about why it's there, and what it does? This may not be possible, but it would greatly strengthen the manuscript if you could do it.
Thanks for asking us to clarify the autoencoder. It is key to ensuring that the error is fedback projected correctly to the neurons for a local learning rule, and to ensure that the neurons are driven to the ideal activities by the error current.
The text added for this clarification is incorporated in the response for the above revision.
See in particular Subsection “Negative error feedback via autoencoder enables local learning”.
“The combination of the autoencoder and the error feedback implies that the output stays close to the reference, as explained now. In open loop i.e. without connecting the output $\vec{\hat{x}}$ and the reference $\vec{x}$ to the error node, an input $\vec{\epsilon}$ to the network generates an output $\vec{\hat{x}} = k\vec{\epsilon}$ due to the autoencoder of gain $k$. In closed loop, i.e. with the output and reference connected to the error node (Figure 1B), the error input is $\vec{\epsilon} = \vec{x}\vec{\hat{x}}$, and the network output $\vec{\hat{x}}$ settles to:
\vec{\hat{x}} = k\vec{\epsilon} = k\left(\vec{x}\vec{\hat{x}}\right)
\implies \vec{\hat{x}} = \frac{k}{k^{+}1}\vec{x} \approx \vec{x},
i.e. approximately the reference $\vec{x}$ for large positive $k$. The fedback residual error $\vec{\epsilon}=\vec{x}/(k^{+}1)$ drives the neural activities and thence the network output. Thus, feedback of the error causes the output $\hat{x}_\alpha$ to approximately follow $x_\alpha$, for each component $\alpha$, as long as the error feedback time scale is fast compared to the reference dynamical system time scale, analogous to negative error feedback in adaptive control [Narendra and Annaswamy, 1989; Ioannou and Sun, 2012].”
5) A potentially important limitation of this paper is that it does not take into account feedback delays during training. This would be a significant limitation if the purpose of this work is to propose how forward models can be learned in the brain: one of the main reasons to have forward models in the first place is being able to deal with feedback delays (see e.g. Wolpert, Ghahramani & Jordan, 1995). Learning that ignores feedback delays is therefore not fully biologically realistic. If this limitation is indeed correct, it needs to be mentioned more prominently in the Discussion section (we only found it buried in the Methods section); if the network can in fact handle feedback delays during training (that are on a similar scale as the target dynamics), it should be demonstrated, as it would make the impact of the paper significantly stronger.
Thank you for pointing out this important use case of the forward model visavis our learning scheme.
We performed these simulations and incorporated the results (Figure 6EH). While our learning scheme cannot really learn a dynamical system followed by a delay taken together as the reference system, we propose a scheme in the Discussion section (Figure 7BC) using the Smith predictor configuration, adapted from suggestions by Miall & Wolpert, 1996 and Smith, 1957, whereby our forward model can learn with and compensate for sensory feedback delays in motor control. Our forward model is able to feedback a corrective prediction to the controller, before the sensory feedback arrives, thus improving control.
In the Results section, we have added panels EH in Figure 6 and added further text at the very end of Results section:
"We also asked if the network could handle sensory feedback delays in the reference signal. Due to the strong limit cycle attractor of the van der Pol oscillator, the effect of delay is less transparent than for the linear decaying oscillator (Figure 2—figure supplement 2), so we decided to focus on the latter. For the linear decaying oscillator, we found that learning degraded rapidly with a few milliseconds of delay in the reference, i.e. if $x(t − \Δ)$ was provided as reference instead of $x(t)$ (Figure 6EF). We compensated for the sensory feedback delay by delaying the motor command input by identical $\Δ$ (Figure 6G), which is equivalent to timetranslating the complete learning protocol, to which the learning is invariant, and thus the network would learn for arbitrary delay (Figure 6H). In the Discussion section, we suggest how a forward model learned with a compensatory delay (Figure 6G) could be used in control mode to compensate for sensory feedback delays."
In the Discussion section, we have added panels B and C in Figure 7 and the text below:
“One of the postulated uses of the forward predictive model is to compensate for delay in the sensory feedback during motor control [Miall and Wolpert, 1996; Wolpert et al., 1995] using the Smith predictor configuration [Smith, 1957]. We speculate that the switch from the closedloop learning of forward model with feedback gain k >> 1 to openloop motor prediction k = 0 could also be used to switch delay lines: the system can have either a delay before the forward model as required for learning (Figure 7B), or after the forward model as required for the Smith predictor (Figure 7C). We envisage that FOLLOW learning of the forward model occurs in closed loop mode (k >> 1) with a delay in the motor command path, as outlined earlier in Figure 6G and now embedded in the Smith predictor architecture in Figure 7B. After learning, the network is switched to motor control mode, with the forward predictive model in open loop (k = 0), implementing the Smith predictor (Figure 7C). In this motor control mode, the motor command is fed with zero delay to the forward model. This enables to rapidly feed the estimated state back to the motor controller so as to take corrective actions, even before sensory feedback arrives. In parallel, available sensory feedback is compared with a copy of the forward model that has passed through a compensatory delay after the forward model (Figure 7C).”
6) An alternative approach to training recurrent networks is to initialize the recurrent weights in some appropriate way and to just train the readout weights during learning (reservoir computing). How much more computationally powerful is the proposed method? A simple way to test this would be to compare the performance of the trained network to that of a network where the recurrent weights have been shuffled after training, the feedback weights set to zero, and the output weights retrained by linear regression. If such a network had comparable computational performance, it would suggest that only the correct recurrent weight distribution is required to learn the task.
Thank you for this interesting comparison to reservoir computing.
We performed these simulations. Our recurrent network does not behave like a reservoir. Rather as shown for revision 2, the neural activities remain lowdimensional as they are driven strongly by the error current.
At the end of subsection “Nonlinear oscillator” (van der Pol) oscillator, we added the text below:
“We also asked whether merely the distribution of the learned weights in the recurrent layer was sufficient to perform the task, or whether the specific learned weight matrix was required. This question was inspired from 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], where the recurrent weights are random, and only the readout weights are learned. To answer this question, we implemented a perceptron learning rule on the readout weights initialized at zero, with the learned network’s output as the target, after setting the feedforward and / or recurrent weights to either the learned weights as is or after shuffling them. The readout weights could be approximately learned only for the network having the learned weights and not the shuffled ones (Figure 2—figure supplement 3), supporting the view that the network does not behave like a reservoir (Methods section).”
We have also added the corresponding text in the Methods section.
7) Does learning still work if the dynamical system has small amounts of stochasticity, or if the error signal is noisy?
Thanks for asking about noise in the reference / error signals. Yes, it is robust.
We have added in the last subsection of Results section, panel B in Figure 6 with accompanying text:
“We added Gaussian white noise to each component of the error, which is equivalent to adding it to each component of the reference, and ran the van der Pol oscillator learning protocol for 10,000 s for different standard deviations of the noise (Figure 6B). The learning was robust to noise with standard deviation up to around 0.001, which must be compared with the error amplitude of the order of 0.1 at the start of learning, and orders of magnitude lower later.”
8) It's critical to have experimental predictions. Or at least a summary of the implications of this work for what's going on in the brain.
Thank you for pointing out this important lacuna. Our two predictions about errorcurrent driven synaptic plasticity and the presence of an autoencoder in a feedback loop have been highlighted now.
We have added in the Discussion section with Figure 7A, a prediction about errorcurrent driven synaptic plasticity:
“A possible implementation in a spatially extended neuron would be to imagine that the postsynaptic error current $I_i^\epsilon$ 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 feedforward and feedback input could be located, as depicted in Figure 7A. Consistent with the picture of a messenger molecule, we lowpass filtered the error current with an exponential filter $\kappa^\epsilon$ of time constant 80 ms or 200 ms, much longer than the synaptic time constant of 20 ms of the filter $\kappa$. Simultaneously, filtered information about presynaptic spike arrival $S_j*\kappa$ is available at each synapse, possibly in the form of glutamate bound to the postsynaptic receptor or by calcium triggered signalling 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.
The separation of the error current from the currents at feedforward or recurrent synapses could be spatial (such as suggested in Figure 7A) or chemical if the error current projects onto synapses that trigger a signalling 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. This leads us to a prediction regarding synaptic plasticity that, say in cortical pyramidal neurons, the plasticity of synapses that are driven by presynaptic input in the basal dendrites, should be modulated by currents injected in the apical dendrite or on stimulation of feedback connections.”
Further, the presence of an autoencoder in a feedback loop is a prediction also added to the Discussion section:
“The first assumption is that error encoding feedback weights and output decoding readout weights form an autoencoder. This requirement can be met if, at an early developmental stage, either both sets of weights are learned using say mirrored STDP [Burbank, 2015], or the output readout weights are learned, starting with random encoding weights, via a biological perceptronlike learning rule [D'Souza et al., 2010; Urbanczik and Senn, 2014]. A prelearned autoencoder in a highgain negative feedback loop is in fact a specific prediction of our learning scheme, to be tested in systemslevel experiments.”
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Thank you for submitting your revised article "Predicting nonlinear dynamics by stable local learning in a recurrent spiking neural network" for consideration by eLife. The evaluation has been overseen by a Peter Latham as the Reviewing Editor and Timothy Behrens as the Senior Editor.
I want to congratulate you on a very nice paper, and thank you for your very thorough responses; you addressed almost all the points. The only problem is that in parts the manuscript is still not very clear (at least not to me). Please take a very hard look at the Methods section "Network architecture and parameters", and make sure it's crystal clear.
We thank you for the compliments and for pointing out the lack of clarity in the Methods subsection: "Network architecture and parameters". We have rewritten the unclear subsection, and we hope that the same is now crystal clear. The changes are outlined below.
We have substantially rewritten the said Methods subsection. It is now renamed "Network parameters" as the network architecture was already described at the start of the Results section. We have partitioned this subsection into smaller subheadings. The key changes are, in brief, that the gains and biases of neurons and the encoding weights for each layer are first clearly defined. We then explain the procedure for initializing them randomly at the start of a simulation, after which they remain fixed. Figure —figure supplement 1 has also been suitably modified to show the different gain functions of the heterogeneous neurons. We have also improved the explanation of setting the readout weights to form an autoencoder with respect to the error encoding weights, as well as improved the discussion of the expansive versus compressive view of our autoencoder.
We have also made a few minor wordlevel corrections / improvements in other sections.
https://doi.org/10.7554/eLife.28295.023Article and author information
Author details
Funding
European Research Council (Multirules 268 689)
 Aditya Gilra
 Wulfram Gerstner
Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (CRSII2_147636)
 Aditya Gilra
 Wulfram Gerstner
Horizon 2020 Framework Programme (Human Brain Project 720270)
 Wulfram Gerstner
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Johanni Brea, Samuel Muscinelli and Laureline Logiaco for helpful discussions, and Samuel Muscinelli, Laureline Logiaco, Chris Stock, Tilo Schwalger, Olivia Gozel, Dane Corneil and Vasiliki Liakoni for comments on the manuscript. Financial support was provided by the European Research Council (Multirules, grant agreement no. 268689), the Swiss National Science Foundation (Sinergia, grant agreement no. CRSII2_147636), and the European Commission Horizon 2020 Framework Program (H2020) (Human Brain Project, grant agreement no. 720270).
Reviewing Editor
 Peter Latham, University College London, United Kingdom
Publication history
 Received: May 3, 2017
 Accepted: November 22, 2017
 Accepted Manuscript published: November 27, 2017 (version 1)
 Version of Record published: December 14, 2017 (version 2)
 Version of Record updated: February 23, 2018 (version 3)
Copyright
© 2017, Gilra et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 4,278
 Page views

 788
 Downloads

 34
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Further reading

 Neuroscience
Resolving trajectories of axonal pathways in the primate prefrontal cortex remains crucial to gain insights into higherorder processes of cognition and emotion, which requires a comprehensive map of axonal projections linking demarcated subdivisions of prefrontal cortex and the rest of brain. Here, we report a mesoscale excitatory projectome issued from the ventrolateral prefrontal cortex (vlPFC) to the entire macaque brain by using viralbased genetic axonal tracing in tandem with highthroughput serial twophoton tomography, which demonstrated prominent monosynaptic projections to other prefrontal areas, temporal, limbic, and subcortical areas, relatively weak projections to parietal and insular regions but no projections directly to the occipital lobe. In a common 3D space, we quantitatively validated an atlas of diffusion tractographyderived vlPFC connections with correlative green fluorescent proteinlabeled axonal tracing, and observed generally good agreement except a major difference in the posterior projections of inferior frontooccipital fasciculus. These findings raise an intriguing question as to how neural information passes along longrange association fiber bundles in macaque brains, and call for the caution of using diffusion tractography to map the wiring diagram of brain circuits.

 Medicine
 Neuroscience
Background: Deep Brain Stimulation (DBS) electrode implant trajectories are stereotactically defined using preoperative neuroimaging. To validate the correct trajectory, microelectrode recordings (MER) or local field potential recordings (LFP) can be used to extend neuroanatomical information (defined by magnetic resonance imaging) with neurophysiological activity patterns recorded from micro and macroelectrodes probing the surgical target site. Currently, these two sources of information (imaging vs. electrophysiology) are analyzed separately, while means to fuse both data streams have not been introduced.
Methods: Here we present a tool that integrates resources from stereotactic planning, neuroimaging, MER and highresolution atlas data to create a realtime visualization of the implant trajectory. We validate the tool based on a retrospective cohort of DBS patients (𝑁 = 52) offline and present single use cases of the realtime platform. Results: We establish an opensource software tool for multimodal data visualization and analysis during DBS surgery. We show a general correspondence between features derived from neuroimaging and electrophysiological recordings and present examples that demonstrate the functionality of the tool.
Conclusions: This novel software platform for multimodal data visualization and analysis bears translational potential to improve accuracy of DBS surgery. The toolbox is made openly available and is extendable to integrate with additional software packages.
Funding: Deutsche Forschungsgesellschaft (410169619, 424778381), Deutsches Zentrum für Luftund Raumfahrt (DynaSti), National Institutes of Health (2R01 MH113929), Foundation for OCD Research (FFOR).