A spiking central pattern generator for the control of a simulated lamprey robot running on SpiNNaker and Loihi neuromorphic boards

Central pattern generator (CPG) models have long been used to investigate both the neural mechanisms that underlie animal locomotion, as well as for robotic research. In this work we propose a spiking central pattern generator (SCPG) neural network and its implementation on neuromorphic hardware as a means to control a simulated lamprey model. To construct our SCPG model, we employ the naturally emerging dynamical systems that arise through the use of recurrent neural populations in the neural engineering framework (NEF). We define the mathematical formulation behind our model, which consists of a system of coupled abstract oscillators modulated by high-level signals, capable of producing a variety of output gaits. We show that with this mathematical formulation of the CPG model, the model can be turned into a spiking neural network (SNN) that can be easily simulated with Nengo, an SNN simulator. The SCPG model is then used to produce the swimming gaits of a simulated lamprey robot model in various scenarios. We show that by modifying the input to the network, which can be provided by sensory information, the robot can be controlled dynamically in direction and pace. The proposed methodology can be generalized to other types of CPGs suitable for both engineering applications and scientific research. We test our system on two neuromorphic platforms, SpiNNaker and Loihi. Finally, we show that this category of spiking algorithms displays a promising potential to exploit the theoretical advantages of neuromorphic hardware in terms of energy efficiency and computational speed.


Introduction
Our work can be placed in the emerging field of Neurorobotics, a field that combines knowledge acquired from different scientific fields and applies them to the study and the control of animal models and robots.Within the context of Neurorobotics, an artificial brain, either biologically or AI inspired, is interacting with a robot model in physical or virtual experiments [1].This enables the testing of hypotheses on virtual embodiment, a concept which encompasses the idea that a brain is not a system isolated from the outer world, but one that constantly receives and processes stimuli and acts according to them.Neurorobotics problems can fall into various categories, for example robotic control based on cerebellar models [2,3], dynamic vision systems based on event-based cameras [4,5], visual perception [6], motor control and locomotion tasks [7,8] and action selection [9].
A major limitation of existing neuronal models that are often used as artificial brains is that they are both energy and computationally demanding, since they are usually running on conventional CPUs.Even though spiking neural network (SNN) models are computationally sparse by definition [10], this characteristic is not taken into account when running them on conventional hardware.Thus specialized hardware that is optimized to run these models has been researched and developed, among others Intel Loihi [11], IBM TrueNorth [12], SpiNNaker [13] and BrainScale [14], the latter two developed within the context of the Human Brain Project.Our work makes use of a SpiNNaker and a Loihi chip that runs the spiking neural network that we developed.
Many fields of robotics have taken inspiration from biological systems, and particularly from the locomotor system.Locomotion of animals is hypothesized to be controlled to a large extent by functional units in the central nervous system (CNS) called called Central Pattern Generators (CPGs) [15,16], which are usually described as neuronal systems that create rhythmic activity patterns with minimal sensory feedback.In vertebrates, these locomotor circuits are located mostly in the spinal cord, and receive stimulation from the brainstem and other areas of the brain such as the motor cortex, the cerebellum and the basal ganglia [17].One interesting finding is that these networks are capable of producing rhythmic output in the absence of feedback with minimal stimulation, even if the spinal cord has been completely isolated from the body [18].The investigation of CPG based locomotion control is motivated by the insight that it can give on animals locomotion systems and by the fact that these kind of bio-inspired controllers present good capabilities in terms of autonomy and modulation [19].So far the CPG approach has been largely validated for the locomotion of snake-like robots [20,21,22,23].On an implementation level there exist several CPG models which are formulated as SNNs, and and these spiking CPGs (SCPGs) are often running on specialized or generic neuromorphic hardware.It was shown that such SCPGs running on Neuromorphic hardware such as FPGAs, SpiNNaker or VLSI are providing a robust and efficient way to control a complex movement [24] including sensory feedback, namely for bipedal walking [25,26], for the movement of an arm [27,28] or to control a six-legged robot [29,30].
The mathematical modelling of CPGs can be categorized into roughly 3 approaches.The first treats the neural circuitry to the abstraction level of biophysical models and incorporates information about ion pumps and ion channels located in the neural cells membrane and their influence on membrane potentials and the generation of action potentials, frequently modelled by Hodgkin-Huxley neuron models.The second approach uses simpler leaky integrate-and-fire neurons as the basis of computation, abstracting away low-level biological information.The third category which is also our starting point is deprived of lower level biological information and treats CPGs as systems of nonlinear coupled oscillators, where one oscillator models the activity of a whole oscillatory neural network at an abstract level.Although conceptually the latter is a phenomenological approach based on the observation of the emerging locomotor patterns, it still offers many explanations of the underlying mechanisms of rhythmic pattern generation.One of the first successful attempts to use a high-level mathematical formulation of a CPG and model it as a dynamical system which can be simulated with spiking neurons was the work of Eliasmith and Anderson [31].Many of the described models are accompanied with neuromechanical simulations that close the loop between body and brain.For an extensive review on CPGs in robotics and biology we refer to [16].
In this article, we present a high-level SCPG for a lamprey robot that was trained to replicate the dynamics of a system of coupled Hopf-like oscillators.This model is able to produce a set of travelling waves with high-level modulation which correspond to a continuous space of swimming gaits.It can run directly on the neuromorphic SpiNNaker and Loihi boards.It builds on the core Neurorobotics idea of interaction between a virtual robot or animal agent and a virtual brain that runs on neuromorphic hardware and achieves a complex locomotion task.In Section 2, we present the underlying mathematical formulation of the system of coupled Hopf-like oscillators as a first step of the modeling, in Section 2.3 we present the spiking version of the CPG and its performance on the two boards.We provide simulations of both the isolated spiking CPG model as well as neuromechanical simulations under different scenarios in 3. We then present our future work (4.1) and a conclusion (4).

Overall model architecture
Locomotor CPGs are modulated by higher level control centers of the brain with low-dimensional control signals, a property which makes CPG models good candidates for robotic control problems.This property of CPGs gives them a role similar to a feed-forward controller inside a control framework, of producing oscillatory signals that are modulated by external stimulation.To test whether our CPG model can successfully control a lamprey robot we implemented a neuromechanical simulation for which we employed an accurate 3D model of a lamprey robot that is composed of nine body parts similar to the Amphibot robot in [32].These parts are bound together by eight joints that have one degree of freedom: the rotation around the vertical axis.To produce the swimming patterns, the angular positions of these joints oscillate with amplitudes, frequencies and phases prescribed by the CPG model.The complete controller architecture can then be divided in three components (see Figure 1): 1. the mesencephalic locomotor region (MLR), that emits high level signals on each side of the spinal cord: the drives; 2. the central pattern generator (CPG), that generates travelling waves for each joint corresponding to the received drives; 3. the proportional derivative (PD) controller, that controls the torques applied to the joints to reach the timevarying target angle positions. .

Oscillatory signals generation based on coupled abstract Hopf-like oscillators
In order to explain the synchronization phenomena between the different oscillatory centers in the vertebrate spinal cord, Ijspeert [7] proposed a model of nonlinear coupled oscillators, and used this model to control a salamander robot.This model proposes a coupling between different oscillatory centers based on coupling weights that dictate the phase difference and frequency of the oscillatory centers.The oscillators can be chained either in a single or double chain.In the double chain model, the one that we employ here, the activity of the one side of the spinal cord is in antiphase with the activity of the other side, a phenomenon which is also observed in measurements of muscle activity of lampreys.Providing different stimuli, coming from the high-level control centers, between the oscillators found on each side can lead to a shift of the overall oscillatory patterns, which when applied to a robot model induces turning due to the change of the overall curvature of the robot.This dynamical system can be described by the following differential equations which describe a system of phase oscillators with controlled amplitude.The oscillators are described first in phase space, which gives an intuition of how the coupling is induced, and then rewritten in Cartesian space which as we explain is a form suitable for modelling with an SNN: In this system the θ i , v i are the phase and the preferred frequency of the i-th oscillator, r i , the amplitude, x i is the output of the i-th oscillator which represents motoneuron activity, and Ψ i is the output of the model that is applied to the robot and combines the activity of the oscillators of left and the right side of the double chained model.From equation 1 one can observe that the first derivative with respect to time of the phase of each oscillator, is modulated by the coupling weights w ij and the amplitude of the oscillators it is connected to.It is interesting to note that when the phase differences Φ ij are reached between the coupled oscillators the term θ j -θ i -Φ ij becomes zero, and thus the oscillator oscillates with the preferred frequency 2πν i .This is indeed the case when the steady state is reached, which takes place when certain convergence criteria are met.Equation 2 describes how the amplitude of each oscillator converges to the preferred amplitude R i , with parameter a i dictating the speed of convergence.This ensures smooth transitions of the amplitude when abrupt changes of the high-level drive occur.Even though this system fully describes a CPG in phase space, it is not suitable for approximation with an SNN, as integrating equation 1 in time, leads to a constantly increasing phase.This constantly increasing value quickly saturates the representational capabilities of neural populations, as they excel in approximating values within a subset of a larger space.The solution for this problem is to reformulate the problem in Cartesian space as follows [33]: where x i , y i denote the x and y-coordinates of a point in 2-D space moving in a circle through time, with frequency controlled by equation 7. The parameter a dictates the speed of convergence of the amplitude to the steady state, and r i it the norm of the [x,y] vector.This formulation is close to the standard form of coupled Hopf oscillators with coupling to other oscillators.This equation has the advantage that the x,y values stay within a limit cycle, whose radius is dictated by the amplitude of the oscillation, solving the problem of continuously increasing phase when one attempts to use the phase representation.
To incorporate the drive corresponding to the high-level stimulation we use two piece-wise linear functions, which saturate when the stimulation is outside of a certain range.These two functions control the target frequency and the target amplitude of each oscillator according to the relations: These two equations replicate biological observations that the frequency and amplitude of muscle contraction increase together with increased stimulation, hence leading to faster locomotion.They complement the CPG with high-level modulation, and with them we have a complete mathematical formulation of the control framework, which we implement in an SNN.
2.3 Implementation of the coupled oscillators system in a spiking network

Architecture of the spiking CPG neural network
The model that we introduced in the previous section is a mathematical formulation of a system of coupled abstract Hopf-like oscillators, modulated in frequency and amplitude by high-level stimulation.We show that such a system can be easily simulated with an SNN simulator.To do so we designed a modular SNN architecture where one oscillatory center is represented by one population of spiking neurons and computes the equations described in (5 -7).This population at the same time encodes equation 9.For the coupling between the neural oscillators we introduce an intermediate population which receives the x,y values from neighbor oscillators, and computes the coupling term of equation 7.This intermediate population facilitates the exchange of data between the neural oscillators, and it's presence is dictated purely by the framework that we chose to implement the SNN.The overall architecture of the model can be seen in Figure 2. At the same time each of the oscillatory centers is receiving input from the high-level drive through equations 8 -9.

Choice of the neural simulator
In order to replicate the system of modulated oscillators with a spiking neural network the choice of a framework that can perform such numerical computations was necessary.A characteristic shared by most neural simulators is that they allow the simulation of simple leaky integrate-and-fire neuron models (LIF).According to this model [34] the neuron spikes when its membrane potential reaches a certain threshold.Each neuron is excited by the neurons that are connected to it either in an excitatory or inhibitory fashion, increasing or decreasing the membrane potential respectively.After a period of inactivity the membrane potential is reset -leaks-to a base value.A neuron is usually connected with multiple other neurons via junctions called synapses.The information flow from one neuron to the other is dictated among other factors by the level of present in the synapse neurotransmitters and whose release is regulated by dedicated proteins.The overall strength of the connection between neurons is dictated by the synaptic weight.
From a computational perspective, the adaptation of the synaptic weights through synaptic plasticity mechanisms is the process which allows these networks of neurons to learn a representation.Synaptic plasticity mechanisms can be either biologically accurate, i.e.STDP [35], or variations of some machine learning inspired approach such as the ones making use of backpropagation algorithms [36], or biologically plausible mechanisms such as the e-prop algorithm [37].Most computational models of spiking neurons employ the simple Leaky integrate-and-fire neuron model.We use these types of neurons for our study as well.Several simulation platforms were suitable for the task of simulating such neurons, but Nengo [38] was chosen for two reasons.First, it has built-in methods for generating neural networks that approximate differential equations.This approach is described in section 2.3.3.Second, it can generate versions of these networks that can run on dedicated neuromorphic hardware, as we discuss in section 2.5.

Nengo and the Neural Engineering Framework
In this section we give an overview of the Neural Engineering Framework (NEF), which is a general methodology for creating neural networks that approximate differential equations [39].Importantly, it generalizes to any neuron model, including LIF spiking neurons, and takes into account the timing of synapses.
To understand the NEF, we start with the standard observation that a normal feed-forward neural network is a function approximator.That is, if we have some input x and some output y, we can train a neural network produce the desired output y = f (x).While this training can be done using any neural network learning algorithm, here we just use the simple method of having a network with a single hidden layer of LIF neurons (no non-linearities at the input or output), randomly generate the first layer of weights, and use least-squares minimization to solve for the second layer of weights.This method works for a large range of functions and is robust to spiking neuron models [39].
However, to generate the CPG model described here, we need networks that approximate differential equations.Here, the NEF applies the following method.Suppose we want the differential equation ẋ = f (x, u).We build a feed-forward network where the inputs are x and u and the output approximates τ f (x, u) + x.We introduce the variable τ here, which will be used as the time constant of a simple exponential low-pass filter synapse that will connect the neurons.Now to generate the recurrent neural network, we simply connect the output of that network back to itself, and scale the u input by τ .The resulting network will approximate ẋ = f (x, u).See [39] for a full proof, which is based on the observation that the Laplace transform of the low-pass filter is F (s) = 1/(1 + sτ ).Similar transformations can be done for more complex synaptic filters, but we do not use those here.
As an example of this process, Figure 4 shows an NEF model of a single Hopf-style oscillator.This was formed by creating a feed-forward single-hidden-layer neural network with three inputs (x, y, and ω) and two outputs (τ (a(R 2 − r 2 )x − ωy) + x and τ (a(R 2 − r 2 )y + ωx) + y).The weights for this network were found by randomly sampling the inputs (x, y, and ω), computing the desired outputs for each input, and then training the network given this data.Afterwards, the resulting input and output connection weights were multiplied together to create the recurrent neural network shown.
The Nengo software toolkit [38], which is the software implementation of the more general Neural Engineering Framework, provides high-level tools for creating such networks for a variety of neuron models.Crucially, it also provides facilities for linking networks together, so that large systems can be built out of these components.Futhermore, the resulting systems can be automatically compiled to run on CPUs, GPUs, or a variety of neuromorphic hardware.

The Nengo model
Based on the third principle of the NEF we employ the dynamical systems that emerge through the use of recursive neurons to implement the oscillators in our model.It is worth noting that recurrent neural populations can implement various dynamical systems, such as integrators, oscillators, even chaotic systems such as Lorenz attractors.The network computes each function from equations(5-9) according to the NEF principles.By doing so the decoded spiking activity of each neural population can be seen as a real-valued vector with the appropriate dimensions.For the populations that encode the oscillators (depicted with theta i in Figure 2) this 4-dimensional vector represents the values [ ẋ, ẏ, ω, R].For the intermediate neuron populations that compute the coupling part of equation 7 the 4-dimensional vector represented is [ ẋi , ẏi , ẋj , ẏj ].The high-level drive is approximated by the decoded activity of a neuronal population dedicated in receiving the drive and translating it to neural activity.A dedicated readout output node (non-spiking) can be used to read the decoded output of the system, that corresponds to the x-coordinate of the Hopf-like oscillator.The complete system with input and output for 4 oscillatory centers can be seen in Figure 3.As will be shown the system can scale to a larger number of oscillatory centers but the scaling can be limited by the capabilities of the neuromorphic hardware that it is running on.
As mentioned in 2.3.3 the Neural Engineering Framework can be used to approximate any linear or non-linear function with spiking activity by computing the connection weights between the different components of a spiking neural network, acting as a neural compiler.This alleviates the need for explicit training of the SNN, as in the NEF the information that needs to be provided is limited to the properties of the neurons (i.e.membrane threshold potential, neuron types), the values that the neural populations need to represent and the functions that they compute, and the NEF solves for the connection weights that will compute the desired functions.This enables specifying the high-level mathematical functions that are encoded by the SNN and that works both for feed-forward as well as for recurrent connections.The latter is particularly relevant for our work as it enables dynamical systems such as the oscillator system that we employ to emerge from the neuronal activity.In order for the connection weights to be computed by the NEF, during the initialization phase of the simulation a random selection of sampling points to be used as inputs to the function to approximate is selected.These points are based on the input space that the neuronal population approximates, f.e. points in the space [0,1] for a population that encodes 1-D values.Then these points are used to generate training data from the functions, by providing the points as inputs to the desired functions and collecting the output.Subsequently a least-squares optimization computes the weights that best fit the decoded neuronal activity to the training data.For a more detailed technical overview of this method we refer the viewer to [40].

Perturbations and robustness of the CPG model
Animal CPGs have been documented to adapt to various perturbations (i.e.external application of a force), by reacting smoothly and exhibiting stable limit cycle behavior, i.e. recovering the gait patterns without losing synchronization.Furthermore different degrees of stimulation of the oscillatory centers on the spinal cord can lead to different gaits.Simple asymmetrical stimulation between the right and left side drive of the spinal cord can induce a shift of the gait patterns to the left or to the right, and can induce turning.We show that these characteristics are exhibited by our model under the following scenarios: 1. Perturbation of a single oscillatory center by external stimulation 2. Asymmetrical stimulation of the spinal cord from left to right side of the spinal cord These scenarios show the CPG model's ability to quickly recover under external perturbations as well as to modulate swimming gaits.

Neuromechanical simulation in the Neurorobotics Platform
To test the output and the high-level adaptation of the control signals we performed a closed-loop neuromechanical simulation of our model with a robot model as a body.The motivation behind simulating our model within a physical simulation framework comes from the fact that neural circuits and control algorithms cannot be separated from their natural habitat, the body.Only within an embodied simulation can we test whether the system that we propose can successfully control a robot.For such a full closed-loop robot-brain interaction simulation we made use of a framework built exactly for this purpose, the Neurorobotics Platform.The Neurorobotics Platform (NRP) is a software simulator developed within the Human Brain Project [41] that enables the synchronization and exchange of data between modelled brains and virtual robots within a physical simulation environment.The Robotic Operating System [42] is the middleware which enables the communication between the different software components, which is also supported by a multitude of physical robots.Within the NRP there is no need for an explicit synchronization mechanism between the physical world and the modelled brain, as such a mechanism is built into the framework.The physical simulation is provided by Gazebo [43], which interfaces with multiple physics engines.It supports directly many different brain simulators such as NEST [44], Nengo and SpiNNaker, and through Nengo one can run models on Loihi.We used this framework to connect the Nengo model presented in section 2.3.4 with the lamprey robot (Figure 1).
To complement the simulation with a simplified fluid dynamics model, we implemented a drag model, which is computing the forces produced by the swimming motion, forcing the robot to move forward.The drag model is the one presented in [45], and computes the forces applied on each robot link based on the formulas: and the coefficients λ can be computed by where υ i and υ i ⊥ are the velocity components of each link relative to the water in the parallel and perpendicular directions.The parameter λ depends on the fluid density ρ and the parameter S i is the surface of the link perpendicular to the link movement.This drag model is only a simple approximation of the fluid forces applied on the robot, but offers simplicity and computational speed compared to the 3D Navier-Stokes equations.

The neuromechanical simulation scenarios
We tested the arising swimming gaits under different simulation scenarios.Firstly we show that the spiking CPG can produce swimming even with a low number of neurons.Secondly we show unperturbed swimming with no high-level modulation.Thirdly, we present modulation of the swimming by the high-level drive with control of direction and speed.To show the ability of the controller to incorporate sensory feedback from the simulation dynamically we add a water speed barrier to the simulation.This speed barrier forces the robot to move to the side without adaptation of the high-level drive, but with modulation the robot manages to overcome it.The water speed barrier is implemented in the form of a global fluid velocity vector opposite to the forward direction.A summary of the scenarios: 1. Unperturbed swimming, effect of varying number of neurons per neural population 2. Unperturbed swimming, no high-level modulation 3. Unperturbed swimming, control of the speed and direction of the robot 4. Presence of water speed barrier, no high-level modulation

Presence of water speed barrier, high-level modulation
The method that we used to modulate the high-level drive of the robot in the presence of a speed barrier consists of a high-level feedback loop that modulates the turning commands (i.e. the left-right asymmetry of drive signals) towards a desired target angle (e.g.similarly to a fish aiming to swim towards a particular far away target).This is implemented through a linear minimization of the error between a target global angle around the z-axis of the robot's head and the actual angle of the robot's head around the z-axis.Thus, when the robot turns i.e. left, the error between the target angle and the measured angle increases and the right drive increases linearly to compensate for the deviation from the target angle.The equations that we used for this strategy: Where the left drive is increased when the error is positive, and the right when negative.u target is the target lateral velocity, R z is the recorded rotation around the z-axis of the robot's head, CF is the correction factor that linearly multiplies the error, and d right0 and d lef t0 provide the baseline of the drive stimulation.This simple error correction strategy proves to be enough to correct the deviation of the robot from a target angle by modulating the CPG with the high-level drive.

Nengo on SpiNNaker-3 and Loihi boards
As stated in [46], the computational limitations for running spiking models on conventional CPUs are originating in the von Neumann architecture.Conventional computers are built and optimized to perform Boolean algebra operations and arithmetic on the data stored in memory.Hence, this data needs to be transferred back and forth between the memory and the CPUs, which can be time consuming.Neuromorphic hardware on the other hand is specialized in running spiking neural networks.The computation takes place in many small calculators that have access to a small amount of local data.This strategy reveals itself to be more time and energy efficient for neuron oriented computations.For this reason, we tested our Nengo model on a SpiNNaker-3 [13] and a Loihi board [11].Due to the direct connection of SpiNNaker and Loihi boards to Nengo with a software interface our model remained high-level but could be run directly on the boards.
It should also be emphasized that, for efficiency reasons, the actual neuron model running on conventional CPUs, SpiNNaker-3, and Loihi, are all slightly different.They can all implement Leaky Integrate-and-Fire neurons (and other neuron models), but they all make slightly different approximations (e.g.fixed-point rounding).This means that the optimal neural network connection weights for these different hardware platforms will all be slightly different.However, because we specify our model in Nengo using only the mathematical function to be approximated, this means that Nengo can take the hardware details into account when solving for the connection weights, and the user does not have to modify their model to adjust for different hardware platforms.
That said, there are still some areas where the Nengo-SpiNNaker and Nengo-Loihi interfaces have room for improvement.
In particular, the software support for automatically splitting a group of neurons to run across multiple hardware cores is lacking, effectively giving an upper limit on the size of a single group of neurons that is hardware-dependent.We also encountered hardware limitations on the amount of data that could be probed (i.e.recorded) during the running of the simulation, as discussed in Section 3.2.3.

Running the isolated CPG model
The first test that we performed on the isolated (i.e.no time-varying external modulation) spinal cord model, shows that our system can produce oscillations and traveling waves from random initial conditions meaning that it exhibits limit cycle behavior.For such a scenario there is a clear periodic activation of the spiking neurons inside the oscillatory populations as can be seen in 6.In order to provide benchmarks for the neuromorphic platforms vs the CPU as well as to show the adaptive capabilities of our model we ran the model with different numbers of neurons and different numbers of oscillatory centers.An interesting finding is that oscillatory patterns are generated even with low numbers of neurons as can be seen in Figure 8.
Furthermore, perturbing the model by providing explicit stimuli on specific oscillatory centers, can lead to some interesting behaviours which show the stability of the circuit.As can be seen in Figure 7 a single external perturbation on one of the oscillatory centers leads to a temporary disruption of the signals, localized around the neighbouring oscillatory centers.Upon removal of the perturbation the oscillators quickly recover and stabilize.This is the limit cycle property of the high-level mathematical model that is captured well by the spiking network, and exhibits the robustness of the model, a property which is of particular importance for robotics problems.
The high-level modulation and control of the signals when varying the input to the network under the scenario described in 2.3.5 can be seen in Figure 5.In this scenario a simple asymmetrical variation of the input signals between the left and the right side of the spinal cord leads to a formulation of different travelling wave patterns, which can induce different swimming behaviours.A variation between the left and right side of the spinal cord leads according to equation 4 to a shift of the center of the signals towards positive or negative angles, which in turn induces a shift of the joints angles towards one side, causing the robot's curvature to change, inducing a change of direction.

Unperturbed swimming
As mentioned in section 3.1 swimming patterns arise even with a smaller number of neurons for every neural population in the spiking neural network, albeit the fewer neurons the less precise the approximation is.A comparison of the three simulation scenarios with consecutively larger numbers of neurons can be seen in videos 3 (500 neurons), 4 (1000 neurons), 5 (2000 neurons).The robot configurations in the scenario of the 2000 neurons can be seen in Figure 9.The videos correspond to Figure 8, and as can be observed the less neurons, the less smooth the swimming is.Nevertheless, even the 280 neurons per neural population are enough to provide a swimming pattern.
Asymmetry of the driving signals between left and right induces turning as can be seen in video 6 , and providing such drives is a simple way to navigate the robot towards one direction.Using a closed loop control method such as the one described in 2.4.1 such asymmetries can be computed and provided automatically to the control loop.

Presence of water speed barrier
As described in section 2.4.1, to demonstrate the controllability of the robot with a closed loop controller we examine the behaviour of the robot with the presence of a speed barrier, first without adaptation of the high-level signal 7 and then with high-level adaptation 8 .In the first video, the speed barrier causes the robot to follow a trajectory towards the side, by applying higher drag forces to the robot in the lateral direction.In this scenario the robot does not manage to compensate for the presence of the speed barrier as the unmodulated oscillatory signals do not induce a correction of the direction of the robot.In the second video on the other hand, the error correction mechanism described in 2.4.1 is activated, causing the trajectory of the robot to be corrected to compensate for the speed barrier, and eventually it manages to orient itself and swim forward.We can observe that the model adapts well when the high-level tonic drive signal is regulated by the error correction mechanism, which conceptually corresponds to the adaptation that a decision making center of the brain would perform in order to follow a certain trajectory.

Energy and computational speed metrics on SpiNNaker-3 and Loihi boards
For robotics applications it is important that the control signals are generated in real-time.In order to be able to control a robot with the two neuromorphic boards that we examined, the quality of the generated signals has to be similar to the one coming from the CPU.Such comparison of the quality for a simulation of 10 secs can be seen in Figures 11 and  10.As can be observed, the signals are of better quality than the CPU for a low number of neurons.The quality of the produced signals depends heavily on the number of neurons that are used to represent them.Due to limitations arising from the architecture of the two neuromorphic boards we tested, the total number of neurons that we could run on a SpiNNaker board is limited to 30000, for a Loihi board the limitations are reached at a similar number of neurons when no probes for measuring the networks output are used.With probes the limit on Loihi is reached at approximately 22000 neurons.The concept of a probe corresponds to a software construct that can be used to collect simulation data from the neuron activity, energy consumption etc.They are used to record the decoded output value of the neural population representing the oscillatory centres.
A more detailed comparison of the runtime performance for the different platforms can be see in figure 12.What we observed during the execution on the neuromorphic chips is that most of the time is spent during phases other than the network execution, mostly during the initialization phase where the network configuration is being setup, and during input-output(I/O) operations such as the transfer of spikes between the neuromorphic board and the host computer.This is especially true for the Loihi board, as can be observed in figure 13, where the actual execution of the network is around 1 second for 10 seconds of simulation time, almost 10 times faster than real-time, slightly increasing as the network's size increases.In contrast, most of the time during execution is spent on other operations such as the exchange of spikes.It is clear, that this is the main bottleneck of Loihi's execution time.SpiNNaker on the other hand, and especially the execution of spiking networks on SpiNNaker through Nengo, is already optimized for real-time execution.This is the reason why the total operation of SpiNNaker including I/O operations and network execution is staying almost real-time.It should be noted that this time also includes waiting times induced by Nengo to make sure the simulation runs in real-time.The network itself is executed on SpiNNaker at around 2 seconds, marking a slightly slower execution time than Loihi.
A more detailed analysis of the time spent during the execution of the network on Loihi during larger simulation times is provided in figure 14.To explain the observations it is useful to separate the operation of the board in three distinct phases.The first would be the initialization and setup phase which includes software overhead, overhead to boot the board, setup of the host server, compilation of neurons and synapses on the board and which is performed only once.
The second phase would be the loading of the spikes into the neuromorphic board which can be done in parallel with the execution of the network, or before the execution of the simulation.The third phase corresponds to the actual execution on the board.From these findings we can conclude that as soon as the execution of the network is separated from the setup it can perform much faster than real-time.It should be noted that these metrics are relevant for this specific neural network and do not provide an accurate metric for other types of models.
Due to software limitations it was not possible to provide accurate energy benchmarks for the SpiNNaker board.However, a comparison of the energy consumption between a CPU and Loihi is provided in figure 15.On Loihi the energy consumption was measured with the built in time and energy probes.For measuring the energy consumption on the CPU, the RAPL interface was used.RAPL is an Intel processor feature that provides the ability of monitoring and controlling the SoC power consumption [47].As the power measurement control domain we used the PACKAGE domain which includes the energy consumption of all cores, integrated graphics and other uncore components like caches and memory controllers.For the actual measurement, a framework developed by [48] was used.
As a result, in figure 15 you can see that the energy consumption of the Loihi chip is by three orders of magnitude lower than executing the same network with Nengo CPU.This shows neuromorphic hardware can deliver significant energy reductions for executing spiking neural networks when compared to traditional CPU architectures.

Conclusions
In this paper we presented a Spiking Central Pattern Generator based on a high-level system of abstract coupled Hopf-like oscillators that can run on both software and neuromorphic hardware.The method which we used can be generalized to any type of similar CPG controller.Our model is highly parametrizable, and is an excellent candidate for optimization methods.With different parametrizations it can provide a vast number of possible synchronized gaits, f.e.travelling and standing waves.Our method enables us to smoothly control a lamprey robot that with regulation of the high-level drive adapts to various simulation scenarios.We presented a closed-loop neurorobotics simulation within the Neurorobotics Platform achieving multiple locomotor tasks.Lastly, we showed that running the controller on neuromorphic hardware can achieve real-time operation and has potential advantages in terms of energy efficiency and computational speed.
Our work is related to other works in the field that attempt to provide insight on the performance of neuromorphic hardware.In particular, SpiNNaker was benchmarked for its performance in terms of energy efficiency and computational speed with similar accuracy, to an HPC system running a full-scale microcircuit of the human cortex model [49].It was shown that for such complex models the energy consumption per synaptic event, which provides an estimate of the energy efficiency is 5.9 µJ, close to the 5.8 µJ consumed by the HPC system.However for simpler models, closer in terms of synaptic connections and number of neurons to the model that we employ, the cost per synaptic event can be as low as 8 nJ [50].Similarly, in [12] they compared the performance of an IBM TrueNorth neuromorphic chip running a set of computer vision neural networks with the performance of a dual 2.4 GHz E5-2440 processor x86 system, as well as a Blue Gene/Q system with up to 32 compute cards and found two to three orders of execution time speedup and five orders of magnitude less energy consumption compared to the non-neuromorphic systems.Blouw et al. [51] showed that the energy performance of Intel's Loihi chip compared to the Movidius Neural Compute Stick, Nvidia's Jetson TX1, a CPU, and a GPU was significantly lower (5.3x,20.5x, 23.2x, 109.1xtimes respectively), for a keyword spotting task.However it should be noted that generating precise energy consumption benchmarks is a cumbersome task, and often the claims about the theoretical energy efficiency of neuromorphic hardware are not accompanied with the corresponding metrics.

Future work
In order to study the challenges presented in animal swimming locomotion, a realistic simulation framework that can model all the different aspects of the physical world is necessary.The dynamics of the system, the control part, and their communication and synchronization is already solved in the Neurorobotics Platform, but a realistic fluid simulation is still missing.We are planning to address this problem and present a unified framework in our future works.This would allow providing realistic force feedback in the control loop, thus enabling the generation of more complex computational models.
Furthermore, our CPG model can be enriched with various form of environmental or sensory feedback, which can be incorporated into the model itself.Sensory data such as stretch receptors, high-level cognitive controllers that regulate the tonic drive are examples of this type of feedback.
One natural continuation of our work would be the transfer of the control framework on a real robot, such as the Amphibot.This is currently limited by the size of the SpiNNaker board that would prevent it from being fitted on the robot.However Loihi comes with a USB stick that is more compact in size and would potentially fit on the robot.One important consideration would be waterproofing the neuromorphic boards, as well as making sure that the changes induced in the dynamics of the system by the extra weight would be negligible.Figure 10: The output of the network for different number of neurons per oscillatory population when executed on SpiNNaker.On SpiNNaker the output of the network is relatively accurate and better than the CPU even for a small number of neurons.The weights were trained with a random seed of 0. Note that high-frequency filtering is applied by default on the output signals.
Figure 11: The output of the network for different number of neurons per oscillatory population when executed on Loihi.The results have similar accuracy as SpiNNaker and perform better than the CPU for a low number of neurons.The weights were trained using the random seed 0. Note that high-frequency filtering is applied by default on the output signals.Python timings refer to the execution of the network from an application developer's point of view and include all the software and IO induced times.The Executing series shows the actual execution time on the chip and is linearly increasing as the number of neurons increase.The Executor series includes both the execution and the transferring of spikes between the board and the CPU.It should be noted that these two processes can be performed in parallel.The times spent during the setup and initialization phases (Host server up, encoding axons/synapses, booting the board, configuring registers) are performed only once and their relative duration is less significant if the simulation time increases, see also 14 Figure 14: Nengo Loihi execution times when the simulation time increases.All the benchmarks were performed with a network with 450 neurons per oscillatory center.In this figure it is evident that the initialization and setup times play an increasingly less significant role as the simulation time increases, making it possible to execute the network in real-time after roughly 35 secs of simulation time.This is important from the perspective of the application developer as it is taking into account all the software and IO bottlenecks, which usually treats the chips as black boxes and optimizes on the software and network layer.From the figure we can observe that the times spent during the operation of the chip are on the transfer of spikes and on the actual execution, which increase linearly in time, whereas all the other times remain relatively stable.

Figure 1 :
Figure 1: The control framework.The brainstem component is abstracting the brain areas that are stimulating the spinal cord, separated into two stimulations, one for each side of the spinal cord.The CPG component, comprised of coupled oscillatory centers organised in a double chain, produces the swimming gaits modulated by the high-level brainstem control.A PD controller is receiving the output of the CPG network and applies it to the robot, controlling the angular rotation of each joint.

Figure 2 :
Figure 2: Architecture of the spiking CPG model.Each oscillatory center, noted with theta i is coupled with its neighbours through an intermediate population, depicted with C ij .The intermediate population is computing the coupling term of equation 7. The x-y diagrams corresponding to each oscillator show the trajectory of a point traversing the limit circle through time for the ideal mathematical model.As can be observed, the oscillators in each side of the spinal cord have an antiphase relationship between them, whereas the ones upwards or downwards have a fixed phase difference of 4π/N umOsc.

Figure 3 :
Figure 3: (Left) The Nengo simulated model where 4 oscillatory centers are shown.In this simulation the high-level stimulation is driving the oscillations.(Right) The output of each oscillator that corresponds to the decoded spiking activity, when 2000 neurons per oscillatory center are used, is depicted.

Figure 4 :
Figure4: The behavior of a single Hopf-like oscillator implemented in spiking neurons using Nengo and the Neural Engineering Framework (NEF).The model consists of an all-to-all recurrently connected layer of LIF neurons with exponential synapses with 100ms time constants.Their spiking activity is shown in the middle row, sorted by similarity.A single input (ω) is provided, and the two outputs show that it functions as a controlled oscillator.The input weights, recurrent weights, and output weights are found using the NEF such that the network approximates ẋ = a(R 2 − r 2 )x − ωy and ẏ = a(R 2 − r 2 )y + ωx.

Figure 5 :
Figure 5: The output of the CPG network for 16 oscillatory centers, where each oscillator is depicted with θ i .An asymmetric drive is provided to the network after 5 seconds of simulation, increasing the drive on the right side of the spinal cord, and decreasing it on the left.As can be observed the amplitude of the oscillations on the right side increases, whereas on the left side decreases.

Figure 6 :
Figure 6: Spike train of the first 50 neurons of an oscillatory population with 2000 neurons for 4 secs.The activity of the neurons shows clears signs of periodicity.The neurons are continuously alternating between high and low firing rates.

Figure 7 :
Figure 7: The output of the network when the 5th oscillator is perturbed by an external signal.The perturbation lasting from 4.8 to 5 secs causes disturbance of the neighbouring oscillators' θ 2 , θ 5 , θ 6 wave patterns.The model quickly recovers when the perturbation is removed.

Figure 8 :
Figure 8: The output of the network for different number of neurons per oscillatory population.Even with 500 neurons the network can produce an oscillatory output, of lower quality as some of the oscillators' waves are not smooth and there is more high-frequency noise.With 100 neurons there is an improvement of the quality of the signals, whereas with 2000 neurons the signals are smooth and without high-frequency noise.Even with a low number of neurons the patterns are capable of producing simulated swimming.The network was trained in Nengo with a random seed of 0.

Figure 9 :
Figure 9: Swimming with the simulated robot, with snapshots at 160ms intervals for the unperturbed non-adaptive scenario.The network consists of 2000 neurons per neural population.The travelling wave is propagated along the robot's body from head to tail.

Figure 12 :
Figure12: Runtime of a 10 seconds experiment for various number of neurons per platform.The total execution time in SpiNNaker is referring to the complete execution cycle from the moment the simulation is launched to the moment the execution data is collected, likewise in Loihi.It is important to note that these values represent the execution of Nengo on the neuromorphic hardware from the perspective of an application developer, treating the hardware as a black box.The SpiNNaker on-chip execution time measures only the time spent on the board for the execution of the network.The Loihi execution measures the execution time reported by Loihi and represents the actual time spent executing the network.The execution + spike transfer represents the execution time plus the time spent during the exchange of spikes between the Loihi board and the CPU.The reasoning behind these benchmarks is to demonstrate that the times spent on the chip are very low compared to real-time and the rest of the times is spent on IO operations or other operations induced by the software.For a more detailed breakdown of the execution times in Loihi see also Figure13.It can be observed that the actual execution time on the boards is much faster than real-time, showing that neuromorphic hardware is a great candidate for running the CPG model in real-time.

Figure 13 :
Figure13: Breakdown of total execution time on the Loihi chip into different parts for 10 seconds of simulation time and increasing neurons.Python timings refer to the execution of the network from an application developer's point of view and include all the software and IO induced times.The Executing series shows the actual execution time on the chip and is linearly increasing as the number of neurons increase.The Executor series includes both the execution and the transferring of spikes between the board and the CPU.It should be noted that these two processes can be performed in parallel.The times spent during the setup and initialization phases (Host server up, encoding axons/synapses, booting the board, configuring registers) are performed only once and their relative duration is less significant if the simulation time increases, see also14

Figure 15 :
Figure 15: Energy Benchmark of the CPG with Nengo Loihi and Nengo CPU, measured with built-in energy probes in Loihi and with the RAPL interface on the CPU.Is it clear that the energy consumption on the chip is orders of magnitude smaller that the consumption on the CPU.