Abstract
This investigation of the leech heartbeat neural network system led to the development of a low resources, real-time, biomimetic digital hardware for use in hybrid experiments. The leech heartbeat neural network is one of the simplest central pattern generators (CPG). In biology, CPG provide the rhythmic bursts of spikes that form the basis for all muscle contraction orders (heartbeat) and locomotion (walking, running, etc.). The leech neural network system was previously investigated and this CPG formalized in the Hodgkin–Huxley neural model (HH), the most complex devised to date. However, the resources required for a neural model are proportional to its complexity. In response to this issue, this article describes a biomimetic implementation of a network of 240 CPGs in an FPGA (Field Programmable Gate Array), using a simple model (Izhikevich) and proposes a new synapse model: activity-dependent depression synapse. The network implementation architecture operates on a single computation core. This digital system works in real-time, requires few resources, and has the same bursting activity behavior as the complex model. The implementation of this CPG was initially validated by comparing it with a simulation of the complex model. Its activity was then matched with pharmacological data from the rat spinal cord activity. This digital system opens the way for future hybrid experiments and represents an important step toward hybridization of biological tissue and artificial neural networks. This CPG network is also likely to be useful for mimicking the locomotion activity of various animals and developing hybrid experiments for neuroprosthesis development.
Introduction
Millions of people worldwide are affected by neurological disorders which disrupt connections between brain and body, causing paralysis or affecting cognitive capabilities. The number is likely to increase over the next few years and current assistive technology is still limited. In recent decades, extensive research has been devoted to Brain-Machine Interfaces (BMIs) and neuroprosthesis in general (Hochberg et al., 2006, ; Nicolelis and Lebedev, 2009), working toward effective treatment for these disabilities. The development of these devices has had and, hopefully, will continue to have a profound social impact on these patients' quality of life. These prostheses are designed on the basis of our knowledge of interactions with neuronal cell assemblies, taking into account the intrinsic spontaneous activity of neuronal networks and understanding how to stimulate them into a desired state or produce a specific behavior. The long-term goal of replacing damaged neural networks with artificial devices also requires the development of neural network models that match the recorded electrophysiological patterns and are capable of producing the correct stimulation patterns to restore the desired function. The hardware set-up used to interface the biological component is a Spiking Neural Network (SNN) system implementing biologically realistic neural network models, ranging from the electrophysiological properties of a single neuron to large-scale neural networks.
Our study describes the development of a neuromorphic hardware device containing a network of real-time biomimetic Central Pattern Generators (CPG). The main goal of this research is to create artificial CPGs that will be connected to ex vivo spinal cord of rats and guinea pigs, thus achieving one main objective of the Brainbow European project (Brainbow, ) toward hybridization. Hardware-based SNN systems were developed for hybrid experiments with biological neurons and the description of those pioneer platforms was reported in the literature (Jung et al., 2001; Le Masson et al., 2002; Vogelstein et al., 2006). The Brainbow project will go further by using a large-scale neural network instead of few neurons to substitute the functions of a biological sub-network. The final goal is the development of a new generation of neuro-prostheses capable to restore the lost communication between neuronal circuitries.
Locomotion is one of the most basic abilities of animals. Neurobiologists have established that locomotion results from the activity of half-center oscillators that provides alternating bursts. The first half-center oscillator was proposed by Brown (). Pools of interneurons control flexor and extensor motor neurons with reciprocal inhibitory connections. Most rhythmic movements are programmed by central pattern-generating networks consisting of neural oscillators (Marder and Bucher, 2001; Ijspeert, 2008). CPGs are neural networks capable of producing rhythmic patterned outputs without rhythmic sensory or central input. CPGs underlie the production of most rhythmic motor patterns and have been extensively studied as models of neural network function (Hooper, 2000). Half-center oscillators control swimming in xenopus, salamander (Ijspeert et al., 2007), and lamprey (Cohen et al., ), as well as leech heartbeat (Cymbalyuk et al., ), as described in numerous publications. One key article on modeling the leech heartbeat system is Hill et al. (), where the Hodgkin–Huxley formalism is used to reproduce the CPG.
The main novelty of this research was to implement the leech heartbeat system neural network with minimum resources while maintaining its biomimetic activity. Indeed, the final application is a hybrid experiment that requires spike detection, spike sorting, and micro-electrode stimulation. All of these modules are implemented in the same digital board. To achieve this, the Hill et al. () model and results were reproduced using a simpler model (Izhikevich, 2004), implemented in an FPGA (Field Programmable Gate Array) board. This digital board made it possible to design a frugal, real-time network of several CPGs (in this case, a network of 240 CPGs implemented on a Spartan6 FPGA board). For instance, this CPG network is capable of mimicking the activity of a salamander, which requires 40 CPGs (Ijspeert, 2001), or developing hybrid experiments (Le Masson et al., 2002) for neuroprosthesis development (Brainbow, ).
The first part of this article describes the biological leech heartbeat system, based on one segmental CPG. The next section focuses on choosing a frugal neuron model to match the same biological behavior. The following section explains the topology of a single neuron and its implementation in the hardware, followed by its extension to a neuron computation core for increasing the size of the neural network. The next stage was to develop a new synaptic model reproducing activity-dependent depression phenomena to fit the biological activity of a leech heartbeat. The architecture of this digital system is then described in full, including the various blocks. Finally, the system was used to design a CPG network, validated by comparing our measurements with ex vivo rat spinal cord locomotion results following pharmacological stimulation.
Materials and methods
Description of the leech biological heartbeat system
All leech heartbeat studies agree that the CPG (Figure 1C) responsible for this activity (Figure 1A) requires few neurons, making it an ideal candidate system for elucidating the various biomechanisms governing CPG behavior.
Figure 1
Modeling studies indicate that the burst duration of a leech heart interneuron in an elemental oscillator is regulated by the interneuron itself and by the opposite interneuron (see L3 and R3 in Figure 1B) (Calabrese,
Nadim et al. (1995) and Olsen et al. (1995) developed a biophysical model of a pair of reciprocally inhibitory interneurons in the leech heartbeat system. This model included synaptic ionic currents based on voltage-clamp data. Synaptic transmissions between the interneurons consist of spike-mediated and graded synaptic currents. The Hill et al. (
System modeling for hardware implementation
State of art
Some previous studies used silicon neurons (Indiveri et al., 2011) to simulate the leech heartbeat system (Simoni et al., 2004; Simoni and DeWeerth, 2007). Sorensen et al. (2004) created a hybrid system of a heart interneuron and a silicon neuron. The silicon neuron provides real-time operation and implements a version of the Hodgkin–Huxley formalism (Hodgkin and Huxley, 1952). However, due to the complexity of the model, it was only possible to use a small number of silicon neurons and, therefore, only one CPG. This study describes the same results using a large CPG network (240 CPGs on a Spartan6 FPGA board), in preparation for future hybrid experiments with different CPGs. For instance, in the salamander model (Ijspeert, 2001), the body CPG consists of 40 interconnected segmental networks.
When a silicon neuron and heart interneuron are connected with reciprocal inhibitory synapses of appropriate strength, they form a hybrid elemental oscillator that produces oscillations remarkably similar to those seen in the living system. Olypher et al. (2006) described the control of burst duration in heart interneurons using a hybrid system, where a living, pharmacologically-isolated, heart interneuron was connected with artificial synapses to a model heart interneuron running in real-time (software). Using an FPGA board will make it possible to operate in real time using a large number of neurons, together with customized systems for various applications (hybrid experiments).
A few studies (Torres-Huitzil and Girau, 2008; Rice et al., 2009; Serrano-Gotarredona et al., 2009; Barron-Zambrano et al.,
Choice and presentation of the izhikevich model
In designing a SNN, the first step is the choice of a biologically realistic model. Indeed, a mathematical model based differential equations is capable of reproducing a behavior quite similar to that of a biological cell. The choice of model was based on two criteria: the family of neurons able to be reproduced and the number of equations. These criteria were used to compare several models, including the Leaky Integrate and Fire model (LIF) (Indiveri, 2007), the Hodgkin–Huxley model (HH), and the Izhikevich model (IZH).
Hill et al. (
The HH model reproduces all types of neurons with good accuracy (spike timing and shape). Its main drawbacks are the large number of parameters and the equations required. In the heartbeat network, the main focus is on excitatory neurons, like RS. The HH model required 32 parameters for an RS and 26 for a fast-spiking neuron (FS) (Grassia et al.,
The IZH represents a good solution, as it is based on two equations and is capable of reproducing many different families of neurons by changing four parameters. Furthermore, according to Izhikevich (2004), this model is resource-frugal, a key advantage when the aim is to design a large CPG network embedded in the same board as other modules required for hybrid experiments (spike detection, spike sorting, stimulation, etc.).
The IZH model depends on four parameters, which make it possible to reproduce the spiking and bursting behavior of specific types of cortical neurons. From a mathematical standpoint, the model is described by a two-dimensional system of ordinary differential equations (Izhikevich, 2003): with the after-spike resetting conditions: In equation (3), v is the membrane potential of the neuron, u is a membrane recovery variable, which takes into account the activation of potassium and inactivation of sodium channels, and IIzh describes the input current from other neurons.
The IZH model was chosen to emulate the behavior of the excitatory cells for its simplicity and its capacity to implement various families of neurons. The next step was to determine the network system topology. The next section describes the design of one neuron and its extension to a neuron computation core, then the different synapse models implemented, and, finally, the topology of the network.
System topology
Topology of one neuron core: architecture and implementation
In order to make the Izhikevich neural network more biomimetic, the IIzh current from equation (1) was split into three: Ibias, Iexc, and Iinh. Ibias is the biasing current, Iexc is the positive contribution due to excitatory synapses, and Iinh is the negative contribution of inhibitory synapses. Those currents will be detailed in Synapse Model. As suggested in Cassidy and Andreou (
Moreover, and, as the time step of the IZH model is equal to one millisecond (Δt = 1): One neuron was implemented on the FPGA board according to these equations and specifications. This neuron was then extended into a neuron computation core that updated the u and v values of all neurons in the network. Consequently, the neuron implementation became a neuron computation core. For instance, around 2000 independent neurons could be implemented on our digital board. In this system, the type of neuron is defined by the four Izhikevich parameters: a, b, c, and d from equations (2) and (3). Moreover, the state of a neuron is defined by values u and v, and the three current values. Those 9 values were saved in a RAM for use in the next millisecond in the step computation. By extension, the same process can be used for every neuron in the network.
Each u and v computation step is run in parallel, using two pipelines based on the architecture presented in [9]. The topology is presented in Figure 2. All parameters from equations (2), (3), and (4), as well as the u and v values used in the computation are synchronized in one cycle before going through the pipelines (not shown in Figure 2).
Figure 2

Architecture of the “u” and “v” pipelines in the neural computation core. The computation cycles are separated by dotted lines.
To resume, each neuron is represented by one “v” and “u” value, four Izhikevich coefficients (a, b, c, and d), and three currents (Ibias, Iexc, and Iinh).
Iexc, Iinh, and Ibias are added in two cycles at the beginning of the “v” pipeline, while the “u” pipeline is still inactive (steps 1 and 2). The current sum is added to the constant 109.375 and at the same time as the first multiplication (step 3). By multiplexing operands, the same multiplier is used for the following multiplications in different computation cycles. In step 4, v2 is obtained by another multiplication. A simple two-bit shift makes it possible to obtain 4v and add it to v. At the same time, u is used in two subtractions. Step 5 consists of a 5-bit shift to obtain (1/32)v2, an addition, and the last multiplication. In step 6, the computation of both u and v is completed. In the next step, the v value is tested against the threshold to determine whether the neuron has emitted a spike or not. This test gives the next u and v values for this neuron to be stored in the RAM.
An RS neuron with a = 0.002, b = 0.2, c = −65, and d = 8 was used to implement the CPG.
Once the neuron computation core was implemented, the synaptic model was chosen and implemented.
Synapse model
A network is defined by a group of neurons and a group of synapses. Once the neuron model had been chosen, it was obviously necessary to choose a synapse model. Like the neuron model, this model had to be biomimetic but frugal in its use of resources. In biology, synapses are described as links between neurons that transmit different types of synaptic currents to each other to either excite or inhibit neuron activity. In our implementation, a synaptic weight (Wsyn) was added to the synaptic current. When Wsyn was positive, it was added to Iexc (excitatory synaptic current) and when Wsyn was negative, it was added to Iinh (inhibitory synaptic current).
Thanks to AMPA and GABA effects, all synaptic current excitations or inhibitions, respectively decay exponentially (Ben-Ari et al.,
The synaptic current Isyn was implemented with a time constant τsyn for the exponential decay, as follows: When computation step T equals one millisecond and τsyn is in ms: Adding the synaptic weight to the synaptic current, the new equation is: The synaptic computation core implementation is based on the same principle as the neuron computation core. However, this model is not adequate to fit biological data. It was, therefore, decided to implement an activity-dependent depression, where the new synaptic weight, Ws, was dependent on Wsyn.
Activity-dependent depression
As the synaptic behavior described in Hill et al. (
This phenomenon has been reported in neurobiology literature but no model had been devised. This paper proposes a model of this activity-dependent depression that was implemented in digital hardware to improving our CPG network.
As previously explained, each time a neuron emits a spike; the synapse adds a synaptic weight (Wsyn) to the synaptic current. At the same time, the factor (δsyn) indicating the level of depression on a synaptic weight increases. Furthermore, δsyn regulates Wsyn. The value of δsyn is between 0 and 1. Consequently, when δsyn equals zero there is no depression on Wsyn and when δsyn equals one there is maximum depression on Wsyn and the synapse is exhausted.
Ws was used instead of Wsyn as the synaptic weight for each synapse. Then, according to the activity-dependent depression effect, when there is a spike, Ws is added to the synaptic current: The other effect of activity-dependent depression is to increase δsyn after each spike, thanks to the percentage dissipation (P).
The regeneration or reloading of synaptic vesicles is represented by δsyn decreasing to zero. Thus, δsyn decays exponentially when no spike is emitted. So, using the method described in Synapse Model: To summarize, all synapses are now represented by (12), (13), (14) and: The main parameters are: synaptic weight, Wsyn; level of depression, δsyn; and percentage dissipation, P. All these parameters are stored in the RAM on the digital board. Furthermore, this computation required greater precision due to the sensitivity of the parameters. The 26-bit signed fixed representation chosen had 1-bit for the sign, 9-bits for the whole numbers, and 16 for the decimals.
Once the neuron and synapse models had been designed, it was possible to develop the neural network topology.
Network topology
Three elementary blocks. The architecture was based on three main blocks: the neuron implemented (or neuron computation core), a synapse, and the RAM. The connectivity between those blocks is shown in Figure 3.
Figure 3

Global architecture of the spiking neural network.
So far, the neuron computation core can update the state (“u” and “v” variables) of each neuron. In the digital network, the role of the synapse is to update all synaptic currents and weights related to the activity of all neurons, so the synapse block exhibits two behaviors (spiking or not). These two behaviors are summarized in Table 1.
Table 1
| When a spike is emitted | When no spike is emitted |
|---|---|
| Synaptic current | Synaptic current |
| Iexc [n + 1] = Iexc [n] + Ws [n] or Iinh [n + 1] = Iinh [n] + Ws [n] | Iexc [n + 1] = Iexc[n] and Iinh [n + 1] = Iinh [n] |
| Activity-dependent depression | Activity-dependent depression |
| δsyn [n + 1] = δsyn [n] + P(1 − δsyn [n]) |
Description of the equations for synaptic currents and activity-dependent depression.
The IZH model has a time step of one millisecond, so the other computation was synchronized with this time step. The new values of u and v, the exponential decay of Isyn, and the new values of each synaptic current are computed in the same millisecond.
Moreover, a biological neural network is composed of Nn neurons and Ns synapses. To define which neuron is connected to which and with which kind of synapse (excitatory or inhibitory), the network is described using two matrixes: connectivity and synaptic weight (see Figure 3). To save RAM, both matrixes are implemented as sparse matrixes with Nn lines. The ith line in the connectivity matrix corresponds to the connectivity of presynaptic neuron Ni to the other neurons. The synapses are identified by the postsynaptic neuron addresses. For example, the connection to neuron Nj is identified by the number j on the ith line. In the worst case, each neuron is connected to itself and all the others, giving Nn columns. Each matrix line ends with a virtual neuron (address Nn + 1). This implementation is not optimum for the worst case, but the gain is significant for biologically plausible networks, where the total number of synapses is at least four times smaller. Marom and Shahaf (2002) and Garofalo et al. (
There is a direct link between the matrixes: the synaptic weight matrix is the same size as the connectivity matrix, i.e., the same number of lines and columns, with the virtual neurons in the same position (Figure 4). The connectivity between two neurons described by the coordinates (k, l) in the connectivity matrix has the weight shown in box (k, l) in the synaptic weight matrix. A third matrix based on the same principle completes the system: the percentage efficiency matrix, which gives the percentage dissipation, P, of each synapse in a network, as defined in the previous section on activity-dependent depression. We will describe now the state machine of the neural network.
Figure 4

Example of matrix design depending on the neural network. Neuron 2 (N2) is connected by an inhibitory synapse to neuron 1 (N1) and by an excitatory synapse to neuron 3 (N3). Then, on line 2 in the connectivity matrix, N2 is connected to N1, N3, and a virtual neuron (VN), indicating the end of the connection. In the synaptic weight matrix, the synapses for neuron 2 (S2) have a negative weight for the inhibitory synapse and a positive weight for the excitatory synapse. Note the correspondence of its position in both matrixes.
Network machine states. The synaptic current is computed in three successive steps:
– EXT state: for closed-loop experiments, we implement this state in which external feedback can interact with the artificial neural network. This first state consists of using the synaptic block to update the synaptic current. In this case, presynaptic spikes are external events (see Figure 3), such as stimulation from biological neurons in the case of neuroprosthesis. This state makes it possible to stimulate each neuron.
– NEUR state: during this step, the neuron membrane (“u” and “v” from Figure 2) and all exponential decay values are computed in parallel.
– SYN state: the last step consists of updating the synaptic current to reflect the presynaptic spikes computed in the NEUR state. These updated current values are used in the EXT state during the next cycle.
The EXT, NEUR, and SYN states must be completed within a one millisecond time step. If the computation of all three states is completed in less than 1 ms, an IDLE state is implemented until the end of the cycle. Moreover, the blocks (neuron computation core and synapse computation core) described in Figure 3 are multiplexed in time to reduce the implementation area in large-scale neural networks.
Our architecture has two main limits: the number of available cycles (Nc) in one millisecond and the size of the RAM used to save all parameters. Two equations derived from these limits determine the maximum size of the implementable neural network, in terms of number of neuron (Nn) and synapses (Ns).
In the EXT state, all synaptic currents are updated in 10 cycles for each neuron, i.e., 10· Nn cycles. Each neuron requires 11 cycles to compute the NEUR state, i.e., 11· Nn cycles. The synaptic current update during state SYN requires 10 cycles per synapse, i.e., 10· Ns. Figure 2 describes 7 cycles for the neuron computation core, but 4 more cycles are required to read and save the various parameters in the RAM.
This leads to the following equation for computing the maximum number of neurons that may be implemented, depending on the number of cycles available: Having built all the component parts of this real-time, biomimetic digital system, it was possible to validate it by several experiments, presented in the following section.
Results
A CPG is defined by the number of neurons and the families of neurons and synapses. The leech heartbeat neural network was simulated by an appropriate CPG configuration.
Hill et al. (
Biological CPGs provide specifications concerning their behavior. Indeed, their activity is characterized by periodic long bursts (lasting many seconds). Each burst begins by a quick rise in spike frequency to a maximum and ends with a low final spike frequency.
Comparison of biological/digital elemental oscillator
The first example of a CPG was the elemental oscillator (with only two neurons). To reproduce activity accurately, it was necessary to obtain the following values: τampa (time constant of the inhibitory synaptic current exponential decay), τreg (time constant of the recovery of synaptic vesicles), and P (percentage dissipation). These values will be the same for each synapse. The following values were chosen to match biological behavior: τcurrent = 100 ms and τreg = 4444 ms (so 1/τcurrent = 0.01 and 1/τreg = 0.0002). The Ibias current was equal to 8 for both neurons. The synaptic weights are −5.1 and the percentages of dissipation are 1.49.
This model was validated by comparing its implementation with the complex model in Hill et al. (
Figure 5

Comparison between elemental oscillator (Figure 1B) bursting activity in the complex model simulated by scilab, as described in Hill et al. (
In this case, the activity of one neuron inhibits the second neuron. Due to activity-dependent depression and the GABA effect, the inhibition ends and lets the second neuron fire again. In both cases (biological modeling system and digital system), the bursting activity was similar in terms of period and duty cycle, thus validating the simplified elemental oscillator with the complex one. The next step was to validate the segmental oscillator and compare its implementation with biological data.
Comparison of biological/digital segmental oscillator
Keeping the time constant, the biological behavior of the eight-neuron network was duplicated using the following parameters. This time, an eight-neuron CPG was implemented using the same values for τcurrent and τreg than as those used for the elemental oscillator. The use of 8 neurons made it possible to maintain the period without variation (see Table 2) by slowing down the two pairs of oscillators with coordination neurons (De Schutter,
Table 2
| Biological system (Hill et al., | Elemental oscillator (digital) | Segmental oscillator (digital) | |
|---|---|---|---|
| Mean period | 10–12 s | 12.6 ± 1.4 s | 11.2 ± 1 s |
| Mean duty cycle | 57.2 ± 2.9% | 54.7 ± 6% | 46.1 ± 6% |
| Mean spike frequency | 11.9 ± 2.1 Hz | 12.1 ± 1 Hz | 11.2 ± 1 Hz |
| Initial spike frequency | 4.3 ± 0.7 Hz | 8.5 ± 0.2 Hz | 8.6 ± 0.4 Hz |
| Peak spike frequency | 17.5 ± 3.2 Hz | 13 ± 0 Hz | 12.5 ± 0 Hz |
| Final spike frequency | 5.8 ± 1.0 Hz | 8.1 ± 0.2 Hz | 9.3 ± 3 Hz |
Comparison of burst characteristics in the two digital implementations and the biological system.
In Figures 1C, 6, L3/R3 and L4/R4 correspond to the two elemental oscillators and are coupled to the L1/R1 and L2/R2 coordination neurons. The connectivity between each neuron is following Figure 1C. The synaptic weights are −7 and the percentage of dissipation is 2.65.
Figure 6

Logic analyzer measurements of the digital eight-neuron CPG. L3 and R3 show the activity of the first oscillator. L4 and R4 show the activity of the second oscillator. L1/R1 and L2/R2 are the coordination neurons.
The mean period, duty cycle, and variations in spike frequency depending on their position in the burst were measured to quantify the overlap of bursting activity (Table 2). The mean period of this digital implementation was similar to biological values. Note that the segmental oscillator exhibited less variation than the elemental system, thanks to its coordination neurons.
Also, in general, the spike frequency of our implementation was similar to that of the biological system. Due to our synapse model, the frequency reached a maximum in each spike burst but remained on a plateau instead of decreasing to the minimum frequency immediately. In the biological system, the behavior described is due to the enhancement and attenuation of variations in conductance. However, the IZH model does not include conductance, so it cannot be as biomimetic as the HH model. This highlights a weak point of the implementation presented here, but even the HH model, Hill et al. (
These experiments validated the implementation of our elemental and segmental oscillators. This table also confirms that designing a biomimetic system was a good choice. Indeed, the variations of the duty cycle and the period for the bursting activity could not be reproduced by bio-inspired oscillators. The next step was to identify one parameter that would modify the bursting activity period, which would be useful in closed-loop applications.
Variation in the mean period depending on one parameter
A CPG is defined here by the number of neurons and the type synapses involved, the static currents of each neuron, the percentage dissipation, and the synaptic efficiency time constant.
Changing the synaptic efficiency time constant τreg modifies the period of each spike burst (Table 3). The variation in τreg affects the period and duration of each burst, as well as the duty cycle and the variability of these parameters: the greater the value of 1/τreg, the longer the mean period of bursting activity.
Table 3
| 1/τreg(ms−1) | Mean period (s) |
|---|---|
| 0.09 | 4.4 ± 1.6 |
| 0.15 | 7.2 ± 1.2 |
| 0.20 | 11.2 ± 1 |
| 0.22 | 12.9 ± 1.1 |
Variation in the mean period depending on the τreg parameter.
The possibility of modifying the period using a single parameter is very useful and was applied in a closed-loop hybrid experiment concerning locomotion behaviors.
FPGA resources
Originally, a CPG consisted of 8 neurons and 12 synapses but 2 additional synapses per CPG were required to create a network of CPG, by connecting CPG to another one. Thus, each CPG consisted of 8 neurons and 14 synapses. In terms of cycles and available memory, this implementation was capable of running 240 CPGs on a Spartan 6 digital board [see equation (15) and Table 4]. The power consumption of one CPG is 8 mW and for CPGs is 20 mW. We could reduce it in the future by designing a custom ASIC. For neuroprosthesis application, the power consumption should be lower than 80 mW/cm2 chronic heat dissipation level considered to prevent tissue damage (Zumsteg et al., 2005).
Table 4
| Resources | Total available | Used for one CPG | Used for 240 CPGs |
|---|---|---|---|
| Slice FF's | 184304 | 1,093 (0.6%) | 1,459 (0.8%) |
| Slice LUT | 92152 | 1,037 (1.2%) | 1,756 (1.9%) |
| DSP48A1 | 180 | 10 (5.6%) | 10 (5.6%) |
| RAMB16BWER | 268 | 1 | 42 (756 kb) |
| Total RAM | 4824 kb | 9 kb (0.2%) | 765 kb (16%) |
Resources required for one CPG on a Spartan 6 digital board.
Comparison with ex-vivo rat spinal cord results using pharmacological stimulation
The final validation of this system consisted of comparing the CPG output with ex vivo physiological data obtained from the spinal cord of newborn rat [postnatal day (P)1–2]. Bursting locomotor-like activity was induced by bath-application of aCSF (artificial cerebrospinal fluid) mixed with N-methyl-DL-aspartate (NMA; 10 μM), serotonin (5HT; 5 μM), and dopamine (DA; 50 μm) (all purchased from Sigma-Aldrich, France).
For the elemental oscillator. Neuron N1 (corresponds to neuron L3 in Figure 1B) is connected to neuron N2 (corresponds to neuron R3 in Figure 1B) by an inhibitory synapse with a synaptic weight of −7 and a percentage dissipation of 12%. The Ibias current is equal to 7 for both neurons.
Figure 7 shows that the digital system fits the biological recordings of the newborn rat spinal cord. The period and duty cycle of the bursting activity are the same, confirming that the digital system was suitable for hybrid experiments. Instead of using pharmacological stimulation, the digital board will be used in the near future to create a hybrid experiment involving the ex vivo spinal cord and the digital CPGs. A closed-loop is also possible thanks to the possibility of changing the mean period of bursting activity by modifying a single parameter (τreg).
Figure 7

Comparison of pharmacological in-vitro spinal cord with digital CPG.
Discussion
One key step in designing a neuroprosthesis is to produce a large, resource-frugal biomimetic SNN. A biologically realistic CPG (i.e., the leech heartbeat system neural network) was implemented with a minimum resource cost in terms of neuron model, while maintaining its biomimetic activity, as shown in the Results. The first step was to model the biological leech heartbeat system using a single, segmental CPG. The next stage was to choose an efficient neuron model that required few resources for its digital implementation but remained biorealistic enough to match the behavior of biological cells. The topology and hardware implementation of a single neuron were then extended to form a neuron computation core built into a large-scale neural network: 240 CPGs on a Spartan6 FPGA board. Furthermore, the new synaptic model proposed reproduced the activity-dependent depression phenomenon, which had only previously been described in biology literature. The architecture of the entire real-time systemwas described in detail. Finally, the system was validated by several experiments comparing both elemental and segmental oscillators with biological data, and comparing the segmental oscillator with ex vivo rat spinal cord stimulated by pharmacological solutions.
The short-term prospect of this work is to improve the stability of the system using another neuron model. Currently our work is focused on the quartic model (Touboul, 2009), which is more stable than the Izhikevich one and also requires few resources. As described in Table 2, this system is subject to variations in duty cycle and mean period, likely to be reduced by using the new model. However, these variations also exist in biology, so it is necessary to study the actual effect of these variations in the biological system to determine whether they should be eliminated or not.
In the medium term, this system will be included in a hybrid experiment using an ex vivo rat spinal cord. The experiment board includes several modules, including an MEA (Micro-Electrode Array) and spike detection block, to detect and record neural activity in the spinal cord. All these modules, together with the CPG network, will be implemented in the same FPGA. Our neurophysiologist colleagues will identify the best spinal cord sites to stimulate and record bursting activity. These sites will be hybridized to the output of the artificial CPG described in this paper and, in turn, its activity will drive the various ventral root outputs of the spinal cord into full locomotor-like activity. These future experiments aim to demonstrate that hybrid artificial/biological networks provide possible solutions for rehabilitating lost central nervous system function.
Our CPG network could be also used to study the locomotion of different animals. Indeed, according to Ijspeert (2001), the locomotion activity of a salamander requires 40 CPGs, so the 240 CPGs implemented on the Spartan 6 digital board would be suitable for studying more complex locomotion. Our system will be used in a closed-loop system with different sensors and actuators.
Conflict of interest statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Statements
Acknowledgments
This work is supported by the European Union's Seventh Framework Programme (ICT-FET FP7/2007-2013, FET Young Explorers scheme) under grant agreement n° 284772 BRAIN BOW (www.brainbowproject.eu).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
References
1
AmariS. (1972). Characteristic of the random nets of analog neuron-like elements. IEEE Trans. Syst.Man Cybern. 2, 643–657. 10.1109/TSMC.1972.4309193
2
AngstadtJ. D.CalabreseR. L. (1991). Calcium currents and graded synaptic transmission between heart interneurons of the leech. J. Neurosci. 11, 746–759.
3
Barron-ZambranoJ. H.Torres-HuitzilC. (2013). FPGA implementation of a configurable neuromorphic CPG-based locomotion controller. Neural Netw. 45, 50–61. 10.1016/j.neunet.2013.04.005
4
Barron-ZambranoJ. H.Torres-HuitzilC.GirauB. (2010). FPGA-based circuit for central pattern generator in quadruped locomotion. Aust. J. Intell. Inform. Process. Syst. 12, 24–29.
5
Ben-AriY.KhazipovR.LeinekugelX.CaillardO.GaiarsaJ. L. (1997). GABAA, NMDA and AMPA receptors: a developmentally regulated ‘ménage à trois’. Trends Neurosci. 20, 523–529. 10.1016/S0166-2236(97)01147-8
6
Brainbow. (2012). Brainbow Project European Union's Seventh Framework Programme (ICT-FET FP7/2007-2013, FET Young Explorers scheme) Under Grant Agreement n° 284772. Available online at: www.brainbowproject.eu
7
BrownT. (1914). On the nature of the fundamental activity of the nervous centres; together with an analysis of the conditioning of rhythmic activity in progression and a theory of the evolution of function in the nervous system. J. Physiol. 48, 18–46.
8
CalabreseR. L. (1995). Half-center oscillators underlying rhythmic movements, in The Handbook of Brain Theory And Neural Networks, ed ArbibM. A. (Cambridge, MA: MIT Press), 444–447.
9
CalabreseR. L.AngstadtJ.ArbasE. (1989). A neural oscillator based on reciprocal inhibition. Perspect. Neural Syst. Behav. 10, 33–50.
10
CassidyA.AndreouA. G. (2008). Dynamical digital silicon neurons, in IEEE Biomedical Circuits and Systems Conference, BioCAS 2008, 289–292. 10.1109/BIOCAS.2008.4696931
11
CohenA. H.ErmentroutG. B.KiemelT.KopelN.SigvardtK. A.WilliamsT. L. (1992). Modelling of intersegmental coordination in the lamprey central pattern generator for locomotion. Trends Neurosci. 15, 434–438. 10.1016/0166-2236(92)90006-T
12
CymbalyukG. S.GaudryQ.MasinoM. A.CalabreseR. L. (2002). Bursting in leech heart interneurons: cell-autonomous and network-based mechanisms. J. Neurosci. 22, 10580–10592.
13
De SchutterE. (ed.) (2000). Computational Neuroscience: Realistic Modeling for Experimentalists. Boca Raton, FL: CRC Press. 10.1201/9781420039290
14
GarofaloM.NieusT.MassobrioP.MartinoiaS. (2009). Evaluation of the performance of information theory-based methods and cross-correlation to estimate the functional connectivity in cortical networks. PLoS ONE4:e6482. 10.1371/journal.pone.0006482
15
GrassiaF.BuhryL.LeviT.TomasJ.DestexheA.SaighiS. (2011). Tunable neuromimetic integrated system for emulating cortical neuron models. Front. Neurosci. 5:134. 10.3389/fnins.2011.00134
16
Guerrero-RiberasR.MorrisonA.DiesmannM.PearceT. (2006). Programmable logic construction kits for hyper-real-time neuronal modeling. Neural Comput. 18, 2651–2679. 10.1162/neco.2006.18.11.2651
17
HillA. A.LuJ.MasinoM. A.OlsenO. H.CalabreseR. L. (2001). A model of a segmental oscillator in the leech heartbeat neuronal network. J. Comput. Neurosci. 10, 281–302. 10.1023/A:1011216131638
18
HochbergL. R.BacherD.JarosiewiczB.MasseN. Y.SimeralJ. D.VogelJ.et al. (2012). Reach and grasp by people with tetraplegia using a neurally controlled robotic arm. Nature485, 372–375. 10.1038/nature11076
19
HochbergL. R.SerruyaM. D.FriehsG. M.MukandJ. A.SalehM.CaplanA. H.et al. (2006). Neuronal ensemble control of prosthetic devices by a human with tetraplegia. Nature442, 164–171. 10.1038/nature04970
20
HodgkinA. L.HuxleyA. F. (1952). A quantitative description of membrane current and its applications to conduction and excitation in nerve. J. Physiol. 117, 500–544.
21
HooperS. (2000). Central pattern generators. Curr. Biol. 10, 176–177. 10.1016/S0960-9822(00)00367-5
22
IjspeertA. (2001). A connectionist central pattern generator for the aquatic and terrestrial gaits of a simulated salamander. J. Biol. Cybern. 84, 331–348. 10.1007/s004220000211
23
IjspeertA. (2008). Central pattern generators for locomotion control in animals and robots: a review. J. Neural Netw. 21, 642–653. 10.1016/j.neunet.2008.03.014
24
IjspeertA.CrespiA.RyczkoD.CabelguenJ. (2007). From swimming to walking with a salamander robot driven by a spinal cord model. Science315, 1416–1420. 10.1126/science.1138353
25
IndiveriG. (2007). Synaptic plasticity and spike-based computation in VLSI networks of integrate-and-fire neurons. Neural Inform. Process. Lett. Rev. 11, 135–146.
26
IndiveriG.Linares-BarrancoB.HamiltonT.Van SchaikA.Etienne-CummingsR.DelbruckT. (2011). Neuromorphic silicon neuron circuits. Front. Neurosci. 5:73. 10.3389/fnins.2011.00073
27
IzhikevichE. M. (2003). Simple model of spiking neurons. IEEE Trans. Neural Netw. 14, 1569–1572. 10.1109/TNN.2003.820440
28
IzhikevichE. M. (2004). Which model to use for cortical spiking neurons. IEEE Trans. Neural Netw. 15, 1063–1070. 10.1109/TNN.2004.832719
29
JezziniS.HillA. A.KuzykP.CalabreseR. L. (2004). Detailed model of intersegmental coordination in the timing network of the leech heartbeat central pattern generator. J. Neurophysiol. 91, 958–977. 10.1152/jn.00656.2003
30
JungR.BrauerE. J.AbbasJ. J. (2001) Real-time interaction between a neuromorphic electronic circuit and the spinal cord. IEEE Trans. Neural Syst. Rehabil. Eng. 9, 319–326. 10.1109/7333.948461
31
KrahlB.Zerbst-BoroffkaI. (1983). Blood pressure in the leech. J. Exp. Biol. 107, 163–168.
32
Le MassonG.Renaud-Le MassonS.DebayD.BalT. (2002). Feedback inhibition controls spike transfer in hybrid thalamic circuits. Nature417, 854–858. 10.1038/nature00825
33
LeeY. J.LeeJ.KimK. K.KimY. B.AyersJ. (2007). Low power cmos electronic central pattern generator design for a biomimetic underwater robot. Neurocomputing71, 284–296. 10.1016/j.neucom.2006.12.013
34
LewisM. A.HartmannM. J.Etienne-CummingsR.CohenA. H. (2001). Control of a robot leg with an adaptive VLSI CPG chip. Neurocomputing38–40, 1409–1421. 10.1016/S0925-2312(01)00506-9
35
Linares-BarrancoB.Sánchez-SinencioE.Rodríguez-VázquezA.HuertasJ. L. (1993). CMOS Analog Neural Network Systems Based on Oscillatory Neurons, in Silicon Implementation of Pulse Coded Neural Networks, eds ZaghloulM.MeadorJ.NewcombR. (Boston: Kluwer Academic Publishers), 199–247.
36
MarderE.BucherD. (2001). Central pattern generators and the control of rhythmic movements. Curr. Biol. 11, 986–996. 10.1016/S0960-9822(01)00581-4
37
MaromS.ShahafG. (2002) Development, learning and memory in large random networks of cortical neurons: lessons beyond anatomy. Q. Rev. Biophys. 35: 63–87. 10.1017/S0033583501003742
38
MatsuokaK. (1987). Mechanism of frequency and pattern control in the neural rhythm generators. Biol. Cybern. 56, 345–353. 10.1007/BF00319514
39
NadimF.OlsenO. H.De SchutterE.CalabreseR. L. (1995). Modeling the leech heartbeat elemental oscillator. J. Comput. Neurosci. 2, 215–235. 10.1007/BF00961435
40
NakadaK. (2003). An analog cmos central pattern generator for interlimb coordination in quadruped locomotion. IEEE Tran. Neural Netw. 14, 1356–1365. 10.1109/TNN.2003.816381
41
NicolelisM. A. L.LebedevM. A. (2009). Principles of neural ensemble physiology underlying the operation of brain-machine interfaces. Nat. Rev. Neurosci. 10, 530–540. 10.1038/nrn2653
42
NorrisB.WeaverA.WenningA.GarciaP.CalabreseR. L. (2007). A central pattern generator producing alternative outputs: phase relations of leech hear motor neurons with respect of premotor synaptic input. J. Neurophysiol. 98, 2983–2991. 10.1152/jn.00407.2007
43
OlsenO. H.CalabreseR. L. (1996). Activation of intrinsic and synaptic currents in leech heart interneurons by realistic waveforms. J. Neurosci. 16, 4958–4970.
44
OlsenO. H.NadimF.CalabreseR. L. (1995). Modeling the leech heartbeat elemental oscillator: II. Exploring the parameter space. J. Comput. Neurosci. 2, 237–257. 10.1007/BF00961436
45
OlypherA.CymbalyukG.CalabreseR. L. (2006). Hybrid systems analysis of the control of burst duration by low-voltage-activated calcium current in Leech heart interneurons. J. Neurophysiol. 96, 2857–2867. 10.1152/jn.00582.2006
46
RiceK. L.BhuiyanP. A.TahaT. M.VutsinasC. N.SmithM. C. (2009). FPGA Implementation of Izhikevich Spiking Neural Network for Character Recognition, in International Conference on Reconfigurable Computing and FPGAs, (Cancun), 451–456.
47
Serrano-GotarredonaR.OsterM.LichtsteinerP.Linares-BarrancoA.Paz-VicenteR.Gómez-RodríguezF.et al. (2009). CAVIAR: A 45k-Neuron, 5M-Synapse, 12G-connects/sec AER Hardware Sensory-Processing-Learning-Actuating System for High Speed Visual Object Recognition and Tracking. IEEE Trans. Neural Netw. 20, 1417–1438. 10.1109/TNN.2009.2023653
48
SimoniM.DeWeerthS. (2007). Sensory feedback in a half-center oscillator model. IEEE Trans. Biomed. Eng. 54, 193–204. 10.1109/TBME.2006.886868
49
SimoniM.CymbalyukG.SorensenM., R.CalabreseR. L.DeWeerthS. (2004). A multi-conductance silicon neuron with biologically matched conductances. IEEE Trans. Biomed. Eng. 51, 342–354. 10.1109/TBME.2003.820390
50
SorensenM.DeWeerthS.CymbalyukG.CalabreseR. L. (2004). Using a hybrid neural system to reveal regulation of neuronal network activity by an intrinsic current. J. Neurosci. 24, 5427–5438. 10.1523/JNEUROSCI.4449-03.2004
51
StillS.TildenM. W. (1998). Controller for a four legged walking machine, in Neuromorphic Systems: Engineering Silicon from Neurobiology, eds SmithL. S.HamiltonA. (World Scientific Publishing Co Pte Ltd), 138–14810.1142/9789812816535_0012
52
StillS.HeppK.DouglasR. J. (2006). Neuromorphic walking gait control. IEEE Trans. Neural Netw. 17, 496–508. 10.1109/TNN.2005.863454
53
TabakJ.SennW.O'DonovanM.RinzelJ. (2000). Modeling of spontaneous activity in developing spinal cord using activity-dependent depression in an excitatory network. J. Neurosci. 20, 3041–3056.
54
Torres-HuitzilC.GirauB. (2008). Implementation of central pattern generator in an FPGA-based embedded system. 18th International Conference on Artificial Neural Networks, Vol. 5164, 179–187. 10.1007/978-3-540-87559-8_19
55
TouboulJ. (2009). Importance of the cutoff value in the quadratic adaptive integrate-and-fire model. Neural Comput. 21, 2114–2122. 10.1162/neco.2009.09-08-853
56
Van Der PolB. (1928). The heartbeat considered as a relaxation oscillation, and an electrical model of the heart. Philos. Mag. 6, 763–775.
57
VogelsteinR. J.TenoreF.Etienne-CummingsR.LewisM. A.CohenA. H. (2006). Dynamic control of the central pattern generator for locomotion. Biol. Cybern. 95, 555–566. 10.1007/s00422-006-0119-z
58
WijekoonJ.DudekP. (2008). Compact silicon neuron circuit with spiking and bursting behavior. Neural Netw. 21, 524–534. 10.1016/j.neunet.2007.12.037
59
ZumstegZ.KemereC.O'DriscollS.SanthanamG.AhmedR. E.ShenoyK. V.et al. (2005). Power feasibility of implantable digital spike sorting circuits for neural prosthetic systems. IEEE Trans. Neural Syst. Rehabil. Eng. 13, 272–279. 10.1109/TNSRE.2005.854307
Summary
Keywords
central pattern generator, biomimetic, neuron model, spiking neural networks, digital hardware, FPGA
Citation
Ambroise M, Levi T, Joucla S, Yvert B and Saïghi S (2013) Real-time biomimetic Central Pattern Generators in an FPGA for hybrid experiments. Front. Neurosci. 7:215. doi: 10.3389/fnins.2013.00215
Received
06 August 2013
Accepted
29 October 2013
Published
21 November 2013
Volume
7 - 2013
Edited by
André Van Schaik, The University of Western Sydney, Australia
Reviewed by
Jorg Conradt, Technische Universität München, Germany; Runchun M. Wang, University of Western Sydney, Australia; M. Anthony Lewis, Qualcomm, QTI, USA
Copyright
© 2013 Ambroise, Levi, Joucla, Yvert and Saïghi.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Timothée Levi, Laboratoire IMS, UMR Centre National De La Recherche Scientifique 5218, Université de Bordeaux, 351 Cours de la Libération, 33405 Talence, France e-mail: timothee.levi@ims-bordeaux.fr
This article was submitted to Neuromorphic Engineering, a section of the journal Frontiers in Neuroscience.
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.