Efficient hardware implementation of the subthalamic nucleus–external globus pallidus oscillation system and its dynamics investigation

Accepted Manuscript Efficient hardware implementation of the subthalamic nucleus-external globus pallidus oscillation system and its dynamics investigation Shuangming Yang, Xile Wei, Jiang Wang, Bin Deng, Chen Liu, Haitao Yu, Huiyan Li

PII: DOI: Reference:

S0893-6080(17)30164-8 http://dx.doi.org/10.1016/j.neunet.2017.07.012 NN 3793

To appear in:

Neural Networks

Received date : 23 February 2017 Revised date : 26 May 2017 Accepted date : 13 July 2017 Please cite this article as: Yang, S., Wei, X., Wang, J., Deng, B., Liu, C., Yu, H., et al., Efficient hardware implementation of the subthalamic nucleus-external globus pallidus oscillation system and its dynamics investigation. Neural Networks (2017), http://dx.doi.org/10.1016/j.neunet.2017.07.012 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Manuscript Click here to view linked References

Efficient hardware implementation of the subthalamic nucleus-external globus pallidus oscillation system and its dynamics investigation Shuangming Yang1, Xile Wei1, Jiang Wang1, Bin Deng1, Chen Liu1, Haitao Yu1* and Huiyan Li2 1

School of Electrical Engineering and Automation, Tianjin University, 300072, P. R. China. 2

School of Automation and Electrical Engineering, Tianjin University of Technology and Educations, 300222, P.R. China. (*Corresponding author: [email protected], Tel: 86-22-27402293)

HIGHLIGHTS We propose a novel method to implement real-time neural system on a reconfigurable device. We engineer digital subthalamic nucleus-globus pallidus system with biologically plausible dynamics. A novel digital implementation approach is proposed with lowered hardware cost. The oscillation dynamics are highlighted and investigated using the proposed hardware platform. The obtained results show that the proposed method is efficient and feasible in the field of neuro-robotic control. ABSTRACT Modeling and implementation of the nonlinear neural system with physiologically plausible dynamic behaviors are considerably meaningful in the field of

computational neuroscience. This study introduces a novel hardware platform to investigate the dynamical behaviors within the nonlinear subthalamic nucleus-external globus pallidus system. In order to reduce the implementation complexities, a hardware-oriented conductance-based subthalamic nucleus (STN) model is presented, which can reproduce accurately the dynamical characteristics of biological conductance-based STN cells. The accuracy of the presented design is ensured by the investigation of the dynamical properties including bifurcation analysis and phase portraits. Hardware implementation on a field-programmable gate array (FPGA) demonstrate that the proposed digital system can mimic the relevant biological characteristics with higher performance, which means the resource cost is cut down and the computational efficiency is improved by introducing the multiplier-less techniques including novel "shift MUL" approach and piecewise linear approximation. The central pattern generator (CPG) coupled by the presented system is also investigated, which can be applied as an embedded intelligent system in the field of neuro-robotic engineering. Keywords: Field-programmable gate array (FPGA), biological oscillation system, efficient digital implementation, modified neuron model, piecewise linear approximation, central pattern generator (CPG). 1. Introduction In recent decades a large amount of knowledge about the functions of the brain have been accumulated, but the mechanism and the functional dynamics of the brain are still poorly understood (Bakken et al., 2016; Eliasmith et al., 2012). The brain is a

dynamic system that is nonlinear at multiple levels of analysis, and the characterization of its nonlinear dynamics is essential for the understanding of brain function. The dynamical oscillation system in brain has been developed which is composed of spiking neurons. Information transmission and coding in nonlinear systems reproduce the information transmission and coding among biological neurons via the precise timing of spikes or a sequence of spikes. Addition of the temporal dimension for information encoding in the nonlinear system presents new insight into the dynamics of the human brain and could lead to compact representations of neural networks (Ghosh-Dastidar & Adeli, 2009; Yu et al., 2013; Cassidy, Georgioub & Andreou, 2013). A large domain of applications would benefit from the dynamical characteristics of the nonlinear systems including studying brain mechanisms, biomedical applications and various neural engineering applications. The subthalamic nucleus-external globus pallidus (STN-GPe) system plays a vital role in the functions and dysfunction of the basal ganglia (BG). Normal information processing within the STN-GPe oscillation system is featured with complicated spatiotemporal patterns of firing activities (Parent & Hazrati, 1995; Nambu, Tokuno, Inase & Takada, 1997; Afsharpour, 1985). Previous studies have revealed that the oscillations in the dopamine-depleted BG may occur in the network consisting of two nuclei which are the STN and the globus pallidus pars externa (GPe) nucleus (Plenz & Kital, 1999). In this field of research, it has been proposed that the variations in the firing patterns of neurons in the STN-GPe system are of vital importance with or without changes in their overall firing rates (Alexander &

Crutcher, 1990; Smith et al., 1998; Matsumura, Kojima, Gardiner & Hikosaka, 1992). Indeed, the STN-GPe system displays synchronous patterns of firing activities, which may be a reflection of emergent, low-frequency oscillation activity within the network (Magill, Bolam & Bevan, 2000; Limousin et al., 1998; Bergman, Wichmann & Delong, 1990). Pathological low-frequency oscillatory activity in the STN-GPe system is associated with resting tremor and may also contribute to other symptoms of Parkinson's disease (PD) (Hamani et al., 2004; Rubchinsky, Park & Worth, 2012). Furthermore, investigating the STN-GPe circuit has attracted increasingly attentions in neural biophysics and neurophysiology (Bevan et al., 2002; Mallet et al., 2008; Bogacz & Gurney, 2007). The modeling of STN-GPe system demonstrates the oscillatory behavior under normal and various diseased states (Holgado, Terry & Bogacz, 2010). The system coupled by STN and GPe nuclei will be responsible for synchronized oscillatory activities in the normal and pathological BG network (Brown, 2003). In order to explore the biological and pathological mechanisms of the STN-GPe system based on the physiological experiments, a number of computational models have been proposed. Albin et al. presented a model and explains the symptoms of Parkinson's disease (PD) based on changes in mean rate of the BG output (Albin, Young & Penney, 1989). A few models of the STN and GPe neuron such as the Albin et al. model are static models, which means that they represent the inputs and outputs of the nuclei as firing rates (Terman et al., 2002). These rate models of the BG are deficient in the investigation for the dynamical interaction between the STN and GPe

nuclei which can reproduce the pathological changes of tremors. Besides, a computational model of the STN-GPE system has been developed based on electrophysiological and anatomical experimental results (Holgado, Terry & Bogacz, 2010). However, the revelation of ionic dynamics underlying the neural biological behaviors cannot be obtained using this model. Recently, Terman and Rubin presented conductance-based models of the STN and GPe nuclei to explore the dynamic interactions (Terman et al., 2002; Rubin & Terman, 2004). The dynamical mechanisms underlying the biological behaviors can be investigated and understood based on the intrinsic characteristics. As a result, this model will be meaningful for exploring the main mechanisms of the dynamic subthalamopallidal system underlying the emergence of correlated rhythmic activities. Previous theoretical STN neuron models can reproduce the biological dynamics in the electrophysiological experiments in a high level of precision. However, the model complexity will limit their implementations for large-scale spiking neural networks. Thus, a new conveniently feasible STN neuron model is significantly essential for the application in the field of neural information coding and network dynamics. The implementation of the nonlinear dynamic neural system is of vital importance and has been studied increasingly in the field of neuromorphic engineering. During the past few decades, neuromorphic engineers have focused on building hardware integrated circuits that can reproduce and mimic the biological behaviors as well as tasks including vision, audition and locomotion. The hardware implementation is superior to the software simulation due to its low power

consumption, compactness and real-time calculation. A number of neuron arrays are combined with neurons which have biological dynamics similar to that investigated in human brain. As a result, there have been a number of studies for hardware neural system, learning mechanisms in hardware and neural networks (Grossi & Pedersini, 2008; Beyeler, Oros, Dutt & Krichmar, 2015; Basu et al., 2010). Neuromorphic VLSI studies have been proposed by Indiveri et al. and Liu et al. to build cognitive or intelligent neuromorphic machines (Indiveri, Chicca & Douglas, 2009; Indiveri et al., 2011; Levi et al., 2008). Although the analog-based implementations are efficient and more convenient for the straight realization of nonlinear function, they are constrained by the high cost as well as significant development cycle. Reconfigurable digital platforms have been applied to implement neuron models using field-programmable gate arrays (FPGAs), which are generic and programmable digital devices to deal with complex logical calculations (Bonabi et al., 2014; Cassidy, Georgioub & Andreou, 2013; Nanami and Kohno, 2016; Nazari, Faez, Amiri & Karami, 2015; Yang et al., 2016, Yang et al., 2015). Although FPGA-based implementation will cost more silicon area as well as power in comparison with the analog-based design, its shorter development time, high density of logic gates and fast configuration speed makes is significantly powerful for the hardware implementations of biological spiking neurons and neural networks. Hardware resource cost, computational speed and biological validity are the major challenges in terms of the digital neuron implementation. The constrained number of available fast multipliers on a FPGA device will limit the design and realization of a complex conductance-based neuron

model as well as the large-scale spiking neural network composed of the nonlinear neuron models. Besides, multipliers are extravagant and high cost under the consideration of delay, power consumption and hardware resource. There are a large number of nonlinear functions as well as the multiplication calculations in the original STN-GPe system, which will limit the cost-efficient implementation of neuron and the corresponding network. In this paper, a STN-GPe system coupled by the STN and GPe neurons is implemented on a reconfigurable digital device with low hardware resource cost, high computational efficiency and low power consumption. It should be noted that a multiplier is a high cost computational block in terms of the hardware area and power consumption. The number of available multipliers on a single chip is also limited, which constrains the application in the implementation of the large-scale neural networks. Therefore, cost-efficient and multiplier-less designs are essential and necessary. A hardware-oriented subthalamic nucleus (HOSTN) digital circuit model is proposed to reproduce the subthalamic dynamical behaviors. Besides, aiming at high-speed realization, low hardware overhead and considerably efficient design of the STN-GPe system model as a major block for complex algorithm and large-scale considerations, a set of piecewise linear approximation (PLA) functions are presented to replace the nonlinear functions in the original model to reduce the hardware expense and improve the hardware feasibility. Since previous study has revealed that the precision analysis is important for the dynamics of the FPGA-based neurons, a precision analysis has also been presented. Experimental results and hardware

synthesis have shown that the proposed modified models can mimic the dynamical characteristics of the biological STN neuron with a high level of precision. The central pattern generator (CPG) is implemented with the coupling STN-GPe system on a FPGA device, which can mimic the bipedal primary gaits. This lays the foundation for the exploration of the neural oscillations that are essential for the neural information coding in the BG and the mechanism investigation of the PD, and is meaningful for the applications in neuro-robotic engineering. The rest of the paper is organized as follows. The next section introduces the proposed model of the biological STN-GPe system. Section 3 investigates the dynamical properties of the STN-GPe system, including both the coupling dynamics and the single neuronal dynamics. The detailed digital implementation approach of the STN-GPe system is presented in Section 4. In this section, a novel CPG network is established, which can mimic the basic bipedal traits. Implementation results are presented in Section 5. Discussions are given in Section 6 and Section 7 finally concludes the entire work. 2. Dynamic subthalamic nucleus-external globus pallidus system The STN-GPe system is regarded as a central pacemaker of the BG and has far-reaching implications for the functional dynamics of the BG. Firstly, in the culture system, the GPe nucleus was essential for the temporal organization of subthalamocortical activities, and the excitation of the GPe nucleus was dictated by the STN. Secondly, elctrophysiological and in situ hybridization experimental results from models of PD have demonstrated that the temporal dynamics of STN neurons

and synchronization of biological activities among STN nuclei is essential for the dysfunction of BG induced by striatal dopamine deficiency. Besides, the STN cell is one of the major parts of the BG, where it receives the cortical input current and interacts with the GPe nucleus. Fig. 1(a) shows the physiological location of the STN nucleus and the pacemaker in the human brain. Fig. 1(b) schematically summarizes the current view of the complex connectivity of the BG under both the normal and the pathological states. Reciprocally connected glutamatergic cells of STN and GPBAergic neurons of the GPe form a key pacemaker for the BG network to modulate its rhythms. The equivalent circuit representation of a patch of a single subthalamic membrane is shown in Fig. 1(c). Fig. 1 is here. 2.1 Dynamic model of the biological STN-GPe system In the STN-GPe network, each STN neuron sends excitation to one or more GPe neurons. Each GPe neuron sends inhibition to the entire GPe population or some part of it, as well as to one or more STN neurons. The synaptic current IGPe→STN that represents synaptic input from the GPe to STN neuron is modeled as j IGPeSTN  gGPeSTN VSTN  EGPeSTN   sGPe

(1)

where gGPe → STN=0.9 nS/μm2 and EGPe → STN=-85 mV are the maximal synaptic conductance and the synaptic reversal potential from the GPe neuron to STN neuron respectively. The variable VSTN represents the membrane potential of the STN neuron and sGPe is the synaptic variable of the GPe neuron. The summation is taken over the presynaptic GPe neurons. Each synaptic variable sj obeys a first order differential

equation as: j j j sGPe '   1  sGPe  H  VGPej   g    sGPe

(2)

j where VGPe is the membrane potential of the GPe neuron j. The synaptic time courses

are controlled by α and β, which are not simply instantaneous switches. The function H∞ (V) is the smooth approximation of the Heaviside step function, which is defined as



H  V   1 1  exp   V   gH   gH 



(3)

Two kinds of synaptic currents are included in the GPe model. They are modeled by the same type of expression as used for the synaptic current IGPe → STN with appropriately renamed synaptic parameters. The parameter values in the synapse models are given in Table 1. The STN cells contain multiple ionic channels with different biological dynamics. The dynamical characteristics are determined by the intrinsic features in the ionic channels and their mutual interaction. In order to reproduce the dynamics of the STN neurons, the kinetic equations are introduced by Terman et al. to describe the activation and inactivation variables of ionic channels (Terman et al., 2002). The proposed model is featured with spike-producing currents (IK and INa), a high-threshold Ca2+ current (ICa) and a low-threshold T-type Ca2+ current (IT), a leak current (IL) and a Ca2+-activated, voltage-independent K+ current (IAHP). This model is meaningful to investigate the mechanisms underlying the biological behaviors in terms of the intrinsic properties and the dynamics. The equation of the STN neurons is of the form:

 dV  I app  I L  I K  I Na  IT  I Ca  I GPe STN  I AHP Cm dt   dn n  n   n n  dt h  h  dh   h h  dt  dr r r   r  r  dt  d Ca      ( I Ca  IT  kCa Ca )  dt

(4)

where the membrane potential V represents the dynamical evolution of the STN neurons. The model features spike-producing currents (IK and INa), a low-threshold T-type Ca2+ current (IT) and a high-threshold Ca2+ current (ICa), a Ca2+-activated, voltage-dependent "afterhyperpolarization" K+ current (IAHP), and a leak current (IL). This model is introduced in Terman et al. (2002). The leak current and the other voltage-dependent currents are described by the Hodgkin-Huxley formalism and given in Table 1.

Table 1. The currents and functions for the conductance-based STN neuron model. Functions

Values

Functions

Values





IL

g L (V  EL )

r V 

1 1  exp   V   r   r 

IK

g K n4 (V  EK )

m (V )

1 1  exp   V  m   m 

INa

g Na m3h(V  ENa )

a (V )

1 1  exp   V   a   a 

IT

gT a3b 2 (V  ECa )

b (r )

1 1  exp   r  b   b   1 1  exp  b  b  

ICa

gCa s 2 (V  ECa )

s (V )

1 1  exp   V   s   s 

 r (V )

 r0   r1 1  exp  V  r   r

IAHP

g AHP (V  EK )

Ca Ca  15









































 n (V )

 n0   n1 1  exp  V  n   n





 h (V )

 h0   h1 1  exp  V  h   h

h V 

1 1  exp   V   h   h 

n V 

1 1  exp   V   n   n 

The Izhikevich neuron model is claimed to be one of the simplest possible neuron models that is able to reproduce main characteristics of the dynamical behavior of a single biological neuron including bursting and spiking activities. Since it can cut down the hardware resource cost considerably in comparison with the conductance-based neuron model while maintaining the critical dynamical characteristics, we use the Izhikevich model for the modeling of GPe neuron. It has been used in several previous studies for the investigation of the coupled nonlinear dynamic system with biological behaviors (Ambroise et al, 2013; Nazari et al., 2015; Reato et al., 2012). In order to develop a bio-inspired STN-GPe model and investigate STN dependent regulation roles in the dynamics of BG, the STN neuron is interacted with the GPe neuron. The phenomenological manner is preceded and the pacemaker composed of STN and GPe neurons is investigated. Previous studies have used the Izhikevich model to reproduce the dynamical behaviors of the GPe model (Thibeault & Srinivasa, 2013). The membrane potential equation of GPe neuron is described as:  dV 2  dt  0.04V  5V  140  u  I   du  a  bV  u   dt

(5)

with the auxiliary after-spike resetting rule V  c if V  30 mV, then  u  u  d

(6)

where V and u represent the membrane potential and the membrane recovery variable

respectively, which are responsible for the activation of potassium ionic currents and the inactivation of sodium ionic currents. The variable u provides negative feedback to V and I is current input composed of synaptic current and external applied current. Two kinds of synaptic currents are contained in the GPe model. They are the excitatory input from the STN neuron ISTN→GPe, and the inhibitory synaptic current from other GPe neurons IGPe→GPe. They are modeled by the same type of expressions used for IGPe→STN in equation (1-3), with the parameters given in Table 2.

Table 2. The parameter values for the synaptic current of the GPe neuron model Parameters

Values

Parameters

Values

α

2.0 msec-1

θg

20.0

β

0.08 msec-1

-57.0

E STN→GPe

0.0 mV

2.0

EGPe→GPe

-100.0 mV

g STN→GPe

0.3 nS/μm2

gGPe→GPe

0.1 nS/μm2

2.2 The hardware-oriented BG network A general method to implement the nonlinear functions is based on the look-up table (LUT). However, a lot of hardware resource is cost by LUTs. For a LUT with a 10-bit address and 30-bit data, the resource cost is up to 30720 bits. As a result, the number of the LUTs demands to be cut down. The cost-efficient design of the conductance-based STN neuron model is an essential task for the hardware-based

implementation in the field of neuromorphic engineering. A mass of nonlinear functions in the original model are big challenges in implementation of neural network and low-cost design. The major purpose of the modifications is to reduce the hardware resource cost. As a result, a novel method to deal with the challenges is presented, aiming at the hardware resource reduction and the computational efficiency improvement. In the proposed method the nonlinear functions in the neuron models are replaced by the piecewise linear functions implemented by shift and add/subtraction operations. The nonlinear parts of the ionic current models are replaced by the piecewise linear functions and the ionic current models can be described as:  I Na  g Na  f1 V   h   I K  g K  g1  n   V  EK    IT  gT  f 2 V   g 2  n    I Ca  gCa  f 3 V  I  AHP  g AHP V  EK   Ca 

(7)

where the functions f1 (V), f2 (V) and f3 (V) are expressed as







f V  V E Na  1  exp   V   m   m   1    3  f 2 V   V  ECa  1  exp   V   a   a   2  f3 V   V  ECa  1  exp   V   s   s  

 

The four gating variables can be expressed as:

 

3

(8)

 dn  dt  n   h1 V   n   g3 V    dh    h V   h  h V   4 h  2  dt  dr   r   h3 V   r   h5 V      dt  d Ca      ( I Ca  IT  kCa Ca )  dt

(9)

The second type of the nonlinear functions is as follows:  4  g1  n   n 2        1 1     g2  r         r     b   1  exp  1  exp  b      b   b      g V   1  0   1 1+exp V       n n n n  3





(10)





The third type of the nonlinear functions is the sigmoid functions, which have s-shaped curves. These types of functions are as follows:   h1 V   1   h2 V   1    h3 V   1   h4 V   1    h5 V   1

1+exp   V       1+exp   V       1+exp   V          1+exp   V         1+exp   V      n

n

h

h

r

r

0 h

1 h





h

h

0 r

1 r





r

r

(11)

The function h5 (V) is approximately a constant 0.041. The PLA method for the first type of the nonlinear functions are plotted in Fig. 2(a1-a3) Four segments of lines are used to fit the second type of the original function, which are shown as Fig. 2(b1-b3). Five segments of lines are used to fit the third type of the original functions, which are shown as Fig. 2(c1-c5). The piecewise linear functions can be described as follows:

b2  b1  k1V  b1 when V  k  k 1 2   b b b2  b1 V  3 2 k2V  b2 when k1  k2 k 2  k3 f PWL V    ...  bn  bn 1  knV  bn when V  k  k n 1 n 

(12)

where ki and bi are the slope and intercept of the piecewise linear functions separately (i=1, 2,…, n), which can be implemented using shift and addition/subtraction operations. Since the multiplication operation is replaced by the “adder” and “shifter” in our study, the value of coefficient ki in the modified functions defined in equation (12) should be a power of 2 (for example: 1, 2, 4 or 0.5, 0.25, etc.). An exhaustive search algorithm is used to find the coefficient values in the piecewise linear functions. The determination for the values of the coefficients in the linear function is of vital importance. The intersection of the adjacent two lines is denoted by the coefficients. Three factors needs to be considered in this method, which are high fitting degree of the original function, high fitting degree of the ionic currents and accurate neural dynamics of the modified model. A cost function for the evaluation of the approximating error is defined as 1 CFRE  n

n

 f i   f i 

i 1

f ori  i 



ori

mod 2

2

(13)

where n is the total sampling points. An exhaustive search algorithm is used in the proposed approach to get a set of coefficients values in the modified functions. The value ranges of the segment points are limited so that the amount of computations can be reduced. If the modification of the nonlinear function cannot meet the standards aforementioned, the segment number of the piecewise linear function will be added

until the three criteria can be guaranteed. Fig. 2 is here. The coefficient values of the modified functions are demonstrated in Table 3. The coefficient values of ki should be a power of 2 in order to implement multiplication operation by using shift and addition/subtraction operations.

Table 3. The coefficient values in the piecewise linear functions. The double dash represents the absence of the values in the corresponding function. a1

a2

a3

a4

a5

a6

a7

a8

f1(V)

-3/256

-3/32

-11/32

-19/32

-1

-0.25

0.25

0.75

f2(V)

1.5

6.125

75/32

0

-31/32

--

--

--

f3(V)

3/256

23/16

39/8

13/8

-9/32

-31/32

--

--

g1(n)

0

3/64

3/16

51/128

--

--

--

--

g2(r)

0

1/16

7/16

20/16

--

--

--

--

g3(V)

3/4096

5/2048

3/512

9/1024

--

--

--

--

h1(V)

1/512

3/256

7/256

9/1024

1/4096

--

--

--

h2(V)

-1/4096

-1/64

-9/128

-1/64

0

--

--

--

h3(V)

-3/128

-13/128

-3/128

-6/4096 0

--

--

--

h4(V)

1/16384 7/512

9/128

1/64

1/8192

--

--

--

h5(V)

0

--

--

--

--

--

--

--

b1

b2

b3

b4

b5

b6

b7

b8

-0.8856

-5.709

-17.4

-27.3

-41.2

-37.35

-39.3

-46.18

f1(V)

f2(V)

112.3

411.5

245.95

161.8

138.56

--

--

--

f3(V)

0.9

81.1

238.2

164.3

138.785

138.56 --

--

g1(n)

0

-0.007

-0.05

-0.137

--

--

--

--

g2(r)

0

-0.001871 -0.0394 -0.165

--

--

--

--

g3(V)

0.076

0.148

0.19

0.153

--

--

--

--

h1(V)

0.147

0.690

1.37

1.02

0.991

--

--

--

h2(V)

0.985

0.18

-2.24

-0.4

0

--

--

--

h3(V)

-0.7757

-6.31

-1.374

-0.078

0

--

--

--

h4(V)

0.0064

0.709

3.2

1.4

0.999

--

--

--

h5(V)

21/512

--

--

--

--

--

--

--

3. Dynamic investigation of the STN-GPe system. The STN-GPe system is essential for the functions and dysfunction of the BG. Normal information processing within the STN-GPe oscillation system is featured with complicated spatiotemporal patterns of firing activities. The STN cell in the STN-GPe system plays a vital role in the oscillation dynamics and the synchronization of the nonlinear system. In this study, we propose a novel HOSTN model to reproduce the dynamical behaviors of the STN neuron instead of the original model. Thus, we first explore the biological dynamics of the HOSTN model to guarantee the precise dynamic behaviors, which includes the bifurcation analysis and phase portraits plot. Then we further investigate the dynamics of the STN-GPe system with relevant biological characteristics.

3.1 Biological dynamics of the HOSTN model In order to validate the reproduction of the neural dynamics of the HOSTN model, the neural dynamical analysis of the proposed model are given in this section by comparing the dynamical behaviors between the original and modified model. The major dynamical characteristics include the ionic current dynamics, the dynamical response under external stimuli, the neural biological behaviors, the bifurcation properties and the phase trajectory of the HOSTN circuit. The V-n, V-h, V-r and V-[Ca] phase portraits of the HOSTN model are given in Fig. 3(a-d). It shows the biological behaviors of the STN neuron. The phase trajectory of the membrane potential and its deviation in the HOSTN model is depicted in Fig. 3(e), which shows a highly consistence of the dynamics between the HOSTN and the original models. Fig. 3(f) shows that the HOSTN model reproduces the consistent response dynamics in comparison with the original model. Fig. 3 is here. The ionic current plays essential roles in the neural dynamics such as the generation of oscillatory activities (Bose et al., 2014). In this paper the comparison of the ionic current between the original and modified neuron models is proposed. Fig. 4(a) and (b) show the ionic dynamics where the ionic current is considered as a function of voltage. The slow variables are set to the limiting values for fixed voltages to compute the steady-state currents. It shows that the ionic dynamics of the modified models are consistent with the original one. The dynamic responses of the original and HOSTN model under the external applied stimuli are depicted in Fig. 4(c) and Fig.

4(d). It reveals that there exists a saddle-node on invariant circle bifurcation, which is composed of two trajectories called heteroclinic tractories connecting the saddle and the node. Since the ionic currents and bifurcation dynamics strongly depend on the nonlinear portions of the neuron model, it is supposed to obtain some deviation between the original and the modified models in which all of the nonlinear portions have been replaced by the linear functions. However, the HOSTN model can still maintain high precision of the neural dynamics. Thus, slight deviation between the original model and the modified ones can be negligible. Fig. 4 is here. Fig. 5(a) and Fig. 5(b) show the spiking response of the STN neuron to external applied currents and its relevant dynamical mechanism. It can be found that under the external applied current stimulus, the STN neuron can exhibit either quiescent state or repetitive spiking state with the variation of the stimulation current. It shows the firing rates of the HOSTN model for the dynamical response under different levels of stimulus are consistent with the original model. Fig. 5(c) and Fig. 5(d) show the dynamical responses of the original and HOSTN model with the normal and pathological synaptic inputs. It shows that the HOSTN model can reproduce both the normal and the pathological dynamics, which are consistent with the original model. Fig. 5 is here. 3.2 Dynamic behaviors and synchronization of the STN-GPe system Bidirectional communication between STN and GPe is important for information processing, synaptic transmission and essential in the functioning of the BG system.

By the advantages of the extensive innervation of the GABAergic neurons of the BG output nuclei, the STN and GPe neurons powerfully influence the communication of the BG with the rest of the brain as shown in Fig. 6. It shows that different phase relationship can be reproduced with different combination of the coupling strengths. The first and second lines show the biological activities of the coupled STN-GPe system using the original model, and the third and fourth lines depict the biological behaviors of the STN-GPe system with the proposed HOSTN neuron model. It shows that the STN-GPe system with the HOSTN model can reproduce the accurate dynamics in comparison with the original model. Fig. 6 is here. Fig. 7(a) and (b) depicts the synchronization behaviors of the STN-GPe system. The synchronization situation of the system is impressed by the parameter values of the synapses gSTN→GPe, gGPe→STN, gGPe→GPe. One standard deviation method to evaluate the synchronization situation between variables A and B is described by:   A, B  

1 N  Ai   A  Bi  B     N  1 i 1   A   B 

(14)

where μA and σA are the mean and standard deviation of variable A, and μB and σB are the mean and standard deviation of variable B respectively. Each variable has N scalar observations. It gives a proper value reference in Fig. 7(a) and Fig. 7(b) to obtain the required synchronicity. It also shows that a low gGPe→STN induces a low synchronicity and a high gGPe→GPe guarantees a high synchronicity. Fig. 7 is here. Although weak beta oscillations exist under the normal state and strong beta

oscillations exist in the pathological state, the upper beta rhythms always exist in the cortical input. In order to investigate the roles of the intrinsic excitation and inhibition of the STN-GPe system with the normal and pathological beta oscillations, the dependency of the firing frequency on the different combinations of the coupling strength gGPe→STN, gGPe→GPe, and gSTN→GPe is shown in Fig. 8. It shows that different combinations of the coupling strength can induce different beta rhythmic activities in the GPe and STN neurons. When the inhibition of the GPe neuron is low (i.e. the value of gGPe→STN is small) and the excitation of the STN neuron is high (i.e. the value of gSTN→GPe is large), the excitatory-inhibitory dynamic system is active. Also, when the basal ganglia oscillate, the firing frequency increases and power decreases as gGPe→STN becomes larger with gSTN→GPe unchanged, and vice versa. It also shows that increased self-inhibition within the GPe neuron contributes to increased upper beta oscillations (21-35Hz) driven by the cortical rhythm, while decrease in the self-inhibition within the GPe facilitates an enhancement of the lower beta oscillations (12-20Hz) induced by the increased excitability of the BG. Fig. 8 is here. 4. Hardware implementation of the STN-GPe system. In this section, the architecture for the digital implementation of the STN-GPe system model is proposed in detail. The dynamical activities of the STN-GPe system can be reproduced by the proposed digital implementation. A digital topology of the CPG network is also proposed in order to generate the desired bipedal traits based on the synchronization dynamics of the STN-GPe system.

4.1 Equation Discretization The hardware pipelines include five blocks to calculate the variables as shown in the corresponding equations. In the proposed system, Euler method is utilized to implement the HOSTN model, which can reduce the hardware resource cost. In the HOSTN model, there are five blocks to calculate (dV)/(dt), (dn)/(dt), (dh)/(dt), (dr)/(dt) and (dCa)/(dt) for neuron model and one block to calculate (dsa)/(dt) for synapse model as shown in the corresponding equations. The derivatives of the variables are discretized into the following form: t  V  n  1  V  n   C  I app  I L  I K  I Na  IT  I Ca  I GPe STN  I AHP  m   n V  n   n  n  n  n  1  n  n   t  n  n   h n  1  h n  t    h V  n   h  n     h   h   r V  n   r  n  r  n  1  r  n   t  r  r  Ca n  1  Ca n  t     I  I  k Ca n     Ca T Ca       sa  n  1  sa  n   t  Aa 1  sa  n  H  V  n    a   Ba sa  n 



(15)



where n represents the number of iterative steps and ∆t stands for the time step in the Euler method. 4.2 Digital implementation of the HOSTN model The pipelining design is essential in the implementation of the high-performance digital circuit that can improve the throughput of the system. Since five neuron variables and one synapse variable exist in the HOSTN model, its hardware topology has six pipelines and six buffers as shown in Fig. 9. The six pipelines are performed with V_p, n_p, h_p, r_p, Ca_p and s_p stages, where V_buf, n_buf, h_buf, r_buf,

Ca_buf and s_buf are the buffer registers for the V, n, h, r, Ca and s values. The buffers are shifted to obtain new values with every rising edge of the clock. The bit number can be calculated based on the usage and precision. The buffer outputs are applied to the related arithmetic blocks to generate the repetitive states based on the variable equations. Synchronization should be guaranteed within the variable pipelines. Thus, if the numbers of stages are not consistent, the delay of each stage is scheduled in order to enable the synchronization of the equations. This step is provided at the latest step of the pipeline stages. The variable pipelines of the HOSTN model are depicted in Fig. 9(b-f), which are used in the pipelining structure of the HOSTN model. In order to determine the bit width of the neural variables, two major points should be considered. They are the parameter value bounds and the spans of the logic shifts. In the STN neuron model, the spans of V is -70 to 40 mV and the minimum bit length for implementing the membrane potential are 9 bits. The bit-width determination is essential because if the maximum logic shifts are not considered the overflow will occur. In order to eliminate the overflow and guarantee the required precision, a bit width of 30 bits is considered which consists of 10 bits for integer parts and 20 bits for the fractional part. Fig. 9 is here. 4.3 Design details of the digital HOSTN implementation The method to realize the "shift MUL" module is shown in Fig. 10(a). The "shift MUL" block is used to implement the multiplication between the two variables with

powers of two using the barrel shifter, which enables multiplierless calculation. The inputs of the "shift MUL" module are the variables "a[n]" and "b[n]". It should be noted that the value of a[n] is expected to be between 0 and 1 for the exact calculation. The bus splitter is employed to split a number into single-bit data and the output ports are numbered from least significant to most significant bit. The "MUX" block is used to realize the multiplexer operation, which is used to choose the input dataflow. In the first data line of the "MUX" block, the value of the input number is zero. When the output of the bus splitter is 1, the second data line is chosen and the variable "b[n]" is shifted rightwards according to the bit number of the output from the bus splitter. The single-bit outputs are then used as the distances of the barrel shifters. There are two inputs of a barrel shifter, i.e., the distance and the input data, which are represented by the black and orange arrows respectively. The value of variable "b[n]" is inputted into the data port of each barrel shifter, and the outputs of the barrel shifters are then added together using a parallel adder to get the final calculation result. Using this block and method, the multiplication between two variables can be implemented in a multiplierless form, which enables a significant reduction of the hardware resource cost. Fig. 10(b) shows the digital structure of the "f1_block" module that is used to implement the piecewise linear function f1(V). The hardware structures of other piecewise linear functions are identical with the "f1_block" module with different coefficient values and segment number. The number “nseg” is the segment number of the piecewise linear functions. The bus builder block is employed to construct outputs from the single-bit input numbers.

Fig. 10 is here. 4.4 Digital implementation of the STN-GPe system and its network model The proposed architecture for the STN-GPe system is presented in Fig. 11(a), which depicts a system block diagram. The "HOSTN neuron" block can be implemented using the proposed digital HOSTN neuronal circuit, and the presented HOSTN neuron model can be used in this case. Three kinds of synaptic current computational blocks are implemented to realize the coupling between the STN and GPe nucleus. There are two pipelines to implement the synaptic current computational blocks, which are pipelines of "sSTN" and "sGPe". The "GPe neuron" block is implemented using the Izhikevich neuron model, which is shown in Fig. 11(b) in detail. In order to describe the subthalamic nucleus–external globus pallidus network under both the normal and the Parkinsonian states, we propose a digital random sparsely connected network model based on the previous study as shown in Fig. 11(c). In this approach, the HOSTN model is used for the reproduction of biological dynamics. As a result, it will be revealed that the biological dynamics under physiological states can be mimicked. The synaptic current computational blocks of the IGPe→GPe, ISTN→GPe and IGPe→STN are implemented using the topology as shown in Fig. 11(a). Fig. 11 is here. 4.5 Central pattern generator with the coupled STN-GPe systems The central pattern generator (CPG) is the neural circuit in nervous system that

can generate rhythmic patterns automatically without any rhythmic inputs. It is generally believed that animal locomotion including invertebrate and vertebrate is controlled by CPGs. CPG is made up of pacemakers, so previous studies have investigated the pacemaker implemented by the neural system. The STN-GPe system is composed of spiking neuron models, and we have witnessed an increment of emphasis in the artificial neural network to biologically plausible SNN recently. However, according to the author's knowledge, the utilization of the subthalamic nucleus–external globus pallidus pacemaker is rare. We introduce the STN-GPe system to implement a CPG, which can give an explanation of the action potential generation in the basal ganglia system in terms of time- and voltage-dependent sodium and potassium conductance. The constructed CPG network based on the proposed STN-GPe system can be applied into practice for robotic walking control. The mathematical model of the CPG network proposed in Fig. 12 is a network of coupled STN-GPe system. The STN neuron is coupled among the different pacemakers. Thus, the membrane potential in the STN model is described by the following equation: Cm

dV  I app  I L  I K  I Na  IT  I Ca  I GPeSTN  I AHP   wij  x j   ge dt j 1

(16)

j i

where i is the index of neurons and wij is the synaptic coupling strength which indicates the reciprocity between neurons and more specifically, the act on neuron i imposed by neuron j. The variable e is the sensory feedback and g is the gain. The sensory feedback is a series of pulse with the phase Φe(i). The parameter wij=0 means that neuron i and j are not coupled with each other; wij>0 means that neuron j has an

excitatory effect on neuron i. Contrarily wij<0 indicates that neuron j has an inhibitory effect on neuron i. The synaptic coupling function we presented is modeled by Ψ(xj)=max(xj, 0). A function of Ψ(xj) is the Heaviside step function in Ref. 46. The coupling weight matrix has the following form: 0 k W  1  k3   k2

k1 0 k2 k3

k2 k3 0 k4

k3  k2  . k4   0

(17)

Fig. 12 is here. where k1, k2, k3 and k4 are the corresponding parameters of the coupling weight. 5. Implementation results The proposed digital circuit on FPGA is implemented on Altera Cyclone-Ⅳ EP4CE115 development system, which utilizes 60 nm process technologies. In the proposed design, the calculation of the fixed-point number is realized in binary form. The oscilloscope photograph of the membrane potential of a single STN neuron using the HOSTN model is shown in Fig. 13. It reveals that the STN neuron using the proposed neuron models can reproduce exactly the relevant biological dynamics of the STN cell. Fig. 13 is here. The resource utilization for the hardware implementation of the original and HOSTN neuron models is summarized in Table 4. Implementation results on the FPGA device show that the HOSTN model can eliminate the use of the embedded multiplier 9-bit elements, while the original model will cost 96% of the total available resource. The proposed digital implementation of the STN-GPe system can improve

the computational speed and reduce the power consumption compared to the original model due to their remove of the nonlinear functions. Table 4. Device utilization of the Altera Cyclone-IV development system. Abbreviations: Logic elements (LEs), combinational functions (CFs), dedicated logic circuit (DLC), total memory bits (TMBs), embedded multiplier 9-bit elements (EM9Es), phase locked loops (PLLs). Original model

HOSTN model

STN-GPe system

Resource

Avaliable

Total LEs

114480

16978(15%)

12590(11%)

18469(16%)

Total CFs

114480

12499(11%)

11216(10%)

16223(14%)

DLC

114480

7888(7%)

8433(7%)

12387(11%)

TMBs

3981312

289886(7%)

2632(<1%)

47507(<1%)

EM9Es

532

512(96%)

0(0%)

0(0%)

1(25%)

1(25%)

1(25%))

68.05MHz

125.38MHz

96.03MHz

Cost function

Original model

HOSTN model

Improvement

Neuron cost NC

1.1867

0.5983

49.58%

Delay D

2.1014

1.1405

45.73%

Resource cost H

27.2%

5.6%

79.41%

PLLs Frequency

4 143MHz

Resource Utilization

The results of hardware implementation in Table 4 show that the HOSTN model is low-cost in comparison with the implementation of the original model. Since the implementation for the large-scale spiking neural networks is based on the

performance including the resource cost and the operating frequency, a cost function is presented to consider both the hardware overhead and working frequency, which is described as:

NC =

a1  H  a2  D a1  a2

(18)

where we consider the same impact of both parameters as a1=a2=1. The hardware overhead H for each neuron is described by

1  H

LE Used CFUsed DLC Used TMBs Used EM9Es Used  2   3   4   5  LE Available CFAvailable DLCAvailable TMBs Available EM9Es Available 1  2  3  4  5 (19)

where we consider the same weights for all the factors as λ1=λ2=λ3=λ4=λ5=1. The delay D is described as D

1  f

1 fUsed

(20)

f Available

where f is the frequency of the neuron model. The improvement of the neuron cost NC, delay D and resource cost H are listed in Table 4. The cost reduction of the HOSTN neuron model is 49.58%, which reveals that the proposed HOSTN is useful in the large-scale implementation of spiking neural networks. Besides, it should be noted that the number of neurons is limited to the available number of embedded multipliers for an FPGA device, which means the network scale will be limited by the EM9Es cost without multiplier-less design. Thus, the digital multiplier-less implementation of the conductance-based neuron model can make significant contributions to the implementation of large-scale neural networks.

The proposed multiplier-less HOSTN model can not only cut down the hardware resource cost, but also enhance the working frequency. It is implemented with an increase in the 40.9% speed-up in comparison with the implementation of the original model on a Cyclone-IV FPGA that uses the straightforward LUT-based approach. It is meaningful because if the operating frequency is enhanced, more virtual neurons can be realized on a fixed number of physical neurons with the same sampling time (dt=1/64+1/256). The proposed CPG network can be applied as an embedded system for the neuro-robotic control. The modular network developed in this study is based on previous studies, which generates motor patterns of the leg rhythms (Pinto, 2007; Golubitsky et al., 1998). The network is symmetrical and each type of bond corresponds to a fixed type of coupling related to particular strength parameters as shown in Fig. 14(a). The previous studies have introduced four primary bipedal gaits with four population networks, which are two-legged hop, walk, two legged jump and run. The phase relation of each population is given in Table 5. The main parameters of the proposed CPG network to reproduce different bipedal primary gaits are given in Table 6. We determine the coupling weights and inner parameters of each population in the following procedure. Firstly, the CPG network is constructed using the subthalamic nucleus–external globus pallidus pacemaker based on equation (23). Secondly, the format of coupling weights matrix in view of the phase relationship between pacemakers is determined. If pacemaker i and j are anti-phase, set wij=wij<0; if neuron i and j are in-phase, set wij=wij>0; if Φij=Φmn and Φi=Φm, set wij=wmn. There

is no rule in other parameter values, however, they should be small in case they may affect the whole network behavior. Thirdly, adjust the parameter value g and Φe(i), and adjust the initial values of the model variables (V1, V2, V3, V4) if needed. Using the above procedure, we can get all the four primary gaits. The hardware implementation results of all the four primary gaits are shown in Fig. 14(b). Fig. 14 is here. Table 5. Four bipedal primary gaits. Suppose (x1(t), x2(t), x3(t), x4(t)) to be a periodic solution with period normalized to 1. Gait

Left leg

Right leg

Two-legged hop

x1(t), x1(t)

x1(t), x1(t)

Walk

x1(t), x1(t+1/2)

x1(t+1/2), x1(t)

Two-legged jump

x1(t), x1(t+1/2)

x1(t), x1(t+1/2)

Run

x1(t), x1(t)

x1(t+1/2), x1(t+1/2)

Table 6. Four bipedal primary gaits. Suppose (x1(t), x2(t), x3(t), x4(t)) to be a periodic solution with period normalized to 1. Gait

Paramters

Two-legged hop

V1, V2, V3, V4 are not changed g=30

Walk

V1=V4=51 V2=V3=-20 g=30

Coupling weight matrix W 0.15 0.15 0.15  0  0.15 0 0.15 0.15   0.15 0.15 0 0.15   0   0.15 0.15 0.15 0.1 0.1 0.12  0  0.1 0 0.12 0.1   0.1 0.12 0 0.1   0.12  0.1  0.1 0  

Two-legged jump

0.12 0.1 0.1  0 0.12 0 0.1 0.1   0.1 0.1 0 0.12    0.1  0.1 0.12 0  

V1=V3=51 V2=V4=-20 g=30

Run

0.1 0.12 0.1  0  0.1 0 0.1 0.12  0.12 0.1 0 0.1   0   0.1 0.1 0.12

V1=V4=51 V2=V3=-20 g=30

To validate the FPGA implementation, the bit-level fixed-point simulation was compared with the numerical simulation on MATLAB. For bit-level simulation, all hardware components were designed and implemented based on the very-high-speed integrated circuit hardware description language (VHDL). Fig. 15 shows the four outputs of the proposed CPG in the gait mode of "run" with FPGA implementation (red lines) and MATLAB simulation (black lines). As shown in Fig. 15 the proposed CPG generated periodic bursting patterns and the difference between the outputs on hardware and software is small. The major cause of the difference is the fixed-point calculation of FPGA implementation. However, the small difference can be reduced by increasing the number of bits in the proposed system. Fig. 15 is here. The cost function CFNMAE to quantify the relative error has been defined in equation (18). Another criterion for the error assessment is defined as

CFNMAE

1 n  ERRABS n i 1   f max  f min 

(21)

which describes a cost function for the normalized mean absolute error. In this normalized cost function, fmax is the maximum value of the modified function and fmin is the minimum value. The function |ERRABS| represents the absolute error of the modified function. The parameter n=1400, which represents 1400 sampling points to get the value of CFNMAE. The values of the error evaluation criteria CFRE and CFNMAE are presented in Fig. 16(a) and (b). Fig. 16(c) depicted the relative error of the variable x in the "shift MUL" block with the fractional bit width increasing. As we can see from Fig. 16(c), when the bit width of the fractional part was less than 10, the relative error was always great so that it’s unfit for the calculation. A compromise was needed to reach between the required computing accuracy and the available hardware resources. Therefore, we specified the bit width of the fractional part as 12-bit which can maintain relatively high accuracy as well as relatively small resource utilization. Fig. 16 is here. 6. Discussions In this work, we demonstrated this concept that an FPGA-based digital neural activity could operate as a neural spike encoder using the compressing information on the firing rate and spikes amplitude. By extending the idea to the network level, we demonstrate that our approach is potentially suitable for the implementation of the nonlinear biological system. The proposed STN-GPe system can generate the patterns that are used to control the robotic locomotion with different kinds of gaits. While the proposed implementation is not the first to demonstrate a neuron circuit built with reconfigurable device (Cassidy, Georgiou & Andreou, 2013; Grossi & Pedersini, 2008;

Beyeler, Oros, Dutt & Krichmar, 2015; Basu et al., 2010; Indiveri, Chicca & Douglas, 2009; Indiveri et al., 2011; Nazari, Faez, Amiri & Karami, 2015; Yang et al., 2015; Yang et al., 2016; Joucla et al., 2016), our proposed work is the first digital design of the STN-GPe system to cut down the hardware cost and improve the computational efficiency relevant to applications in efficient neural network realization in fields of neuromorphic engineering, to co-integrate hardware neurons with digital synapses in a fully embedded system, and to analyze and demonstrate the Parkinsonian state over characteristic dispersions in device behavior. These advances significantly enhance the prospects for nonlinear neural systems built from reconfigurable digital circuit. The

existing

studies

for

the

implementation

of

the

high-dimension

conductance-based neuron model relied upon a LUT that realize the nonlinear functions, yet a cost-efficient implementation can be achieved by replacing the nonlinear parts with the piece-wise linear approximations. Additionally, the proposed digital circuit of the STN cell could be coupled with the digital GPe neuron to form a STN-GPe system as shown in Fig. 10. In this scheme, each neuron is also implemented on a reconfigurable FPGA that uses the "shift MUL" block to implement a multiplier-less circuit while the digital GPe neuron is also coupled by itself using a chemistry synapse. This circuitry suffices to reproduce the beta rhythm of the BG with a high precision. While an on-board STN-GPe system comes with only two coupled neurons, one such system can be extended to a large scale using the network-on-chip technique. Additionally, our system may be integrated with low cost biocompatible interface as a neuro-prosthetic device (Bonifazi et al., 2013; Moxon &

Foffani, 2015; Heer et al., 2007). Associating with high spatiotemporal-resolution techniques and an ideally seamless interaction with biological cells could allow a considerable cost and energy effective “smart prosthetic” feature, essential for many biomedical and ubiquitous neuro-machine interface applications. The proposed architecture allows for the applications in the closed-loop electrophysiological experiment to explore brain dynamics at the cellular and microcircuit levels. Electrophysiological experiment is one of the broadest and most successful applications which use intracellular experimental techniques including voltage clamp and dynamic clamp (Williams & Mitchell, 2008; Raikov, Preyer & Butera, 2004). As demonstrated, the proposed HOSTN model can regenerate the biological dynamics accurately under a stimulation, which makes them especially attractive to closed-loop dynamic-clamping experiment. Such systems may be able to solve the challenges that can limit the performance and application of dynamic clamp techniques. Hardware-based dynamic clamp systems are limited by their poor programmability, while software-based systems cannot guarantee the real-time performance. All software-based systems based on general personal computers (PCs) suffer from design limitations due to the performance restraints of the traditional PC. Digital implementation of the STN-GPe system can achieve promising simulation results. Another approach to use is to combine with the estimation algorithm to predict the hidden properties. The online estimation algorithm is always complex and needs to be conducted in real time. Previous study has proposed the method to estimate the hidden properties based on the digital neuron model using the UKF algorithm, which

can be used in the estimation of the subthalamic hidden dynamics (Yang et al., 2017). Since the proposed STN-GPe model is cost-efficient and convenient for digital implementation, it will reduce the computational burden of the estimation algorithm. Additional applications would be realized using a dedicated embedded system. The model-based STN-GPe modulator can be used as an alternative solution of the open-loop DBS that improve the control performance and decrease control energy, and exploited its time-dynamics as a minimum BG pacemaker prosthetic. This approach has been shown to perform well with the digital circuits (Yang et al., 2016; Piri, Amiri & Amiri, 2015; Amiri et al., 2016), and may offer a flexible solution space between static and dynamical solutions. Our CPG study shows that the STN-GPe system can be implemented on-chip to form a CPG network for the locomotion control by a number of pacemakers corresponding to the number of outputs (legs). It paves the way for the further investigation of the nonlinear dynamics of the CPG system, and in turn facilitates the exploration of the pathological state of the BG system.

7. Conclusions In this study, the STN-GPe system has been initially modeled for the high-performance hardware implementation and realized on a reconfigurable FPGA platform with high computational precision. The digital structure has been presented and implemented based on cost-efficient design of hardware and computation. A HOSTN model targeting low cost digital implementation have been proposed, which

has lowered computational and hardware cost in comparison with the original conductance-based STN neuron model. The proposed model is hardware-oriented and conveniently implemented on digital circuit, which can reproduce the relevant biological dynamics of STN cells. The results demonstrated that the presented STN-GPe system can reproduce relevant dynamical behaviors and are suitable for digital implementation. A novel CPG network is realized using the coupled STN-GPe systems, which can reproduce the primary bipedal traits. It can be applied in some fields including biomedical engineering and neuro-robotic control. In all these cases the system shows high performance. Our physically realized neural system is the primitive for a computationally universal system, for which studies are currently ongoing. Ultimately, our work aims towards the implementation of high-performance neuromorphic systems with low production cost, an electrophysiological tool for the modulation of the biological subthalamic-pallidal system, and a dedicated neuro-robotic controller to tackle unconventional problems in unconventional electronics environments. Appendix Table A1. The parameter values for the conductance-based STN neuron model. Parameters

Values

Parameters

Values

gL

2.25 nS/μm2

θh

-39.0

gK

45.0 nS/μm2

θn

-32.0

gNa

37.5 nS/μm2

θr

-67.0

gT

0.5 nS/μm2

θa

-63.0

gCa

0.5 nS/μm2

θb

0.4

gAHP

9.0 nS/μm2

θs

-39.0

EL

-60.0 mV

-57.0

EK

-80.0 mV

-80.0

ENa

55.0 mV

σm

15.0

ECa

140.0 mV

σh

-3.1

kCa

22.5

σn

8.0

k1

15.0

σr

-2.0

Ε

3.75×10-5 msec-1

σa

7.8

Φh

0.75

σb

-0.1

Φn

0.75

σs

8.0

Φr

0.2

-3.0

500.0 msec

-26.0

100.0 msec

-2.2

θm

1.0 msec

α

5.0 msec-1

1.0 msec

β

1.0 msec-1

17.5 msec

θg

30.0

40.0 msec

-39.0

-30.0

8.0

ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (Grant Nos.61471265, 61372010 and 61671320).

REFERENCE Ambroise, M., Levi, T., Joucla, S., Yvert, B., & Saïghi, S. (2013). Real-time biomimetic central pattern generators into FPGA for hybrid experiments. Frontiers in Neuroscience, 7, 215 Amiri, M., Amiri, M., Nazari, S., & Faez, K. (2016). A new bio-inspired stimulator to suppress hyper-synchronized neural firing in a cortical network. Journal of Theoretical Biology, 410, 107-118. Albin, R. L., Young, A. B., & Penney, J. B. (1989). The functional anatomy of basal ganglia disorders. Trends in neurosciences, 12(10), 366-375. Afsharpour, S. (1985). Topographical projections of the cerebral cortex to the subthalamic nucleus. Journal of Computational Neurology, 236(1), 14-28. Alexander, G. E, & Crutcher, M. D. (1990). Functional architecture of basal ganglia circuits: neural substrates of parallel processing. Trends in neurosciences, 13(7), 266-271. Bakken, T. E., Miller, JA., Ding, SL., Sunkin, SM., Smith, KA., & Ng, L. et al. (2016). A comprehensive transcriptional map of primate brain development. Nature, 535(7612), 367-375 Bergman, H., Wichmann, T. & DeLong, M. R. (1990). Reversal of experimental parkinsonism by lesions of the subthalamic nucleus. Science, 249(4975), 1436–1438. Bevan, M. D., Magill, P. J., Terman, D., Bolam, JP., & Wilson, CJ. (2002). Move to the rhythm: oscillations in the subthalamic nucleus–external globus pallidus network. Trends in neurosciences, 25(10), 525-531. Bogacz, R., & Gurney, K. (2007). The basal ganglia and cortex implement optimal decision making between alternative actions. Neural computation, 19(2), 442-477. Bonabi, S.Y., Asgharian, H., Safari, S., Ahmadabadi, M.N. (2014). FPGA implementation of a biological neural network based on the Hodgkin-Huxley neuron model. Frontiers in Neuroscience, 8, 379. Bonifazi, P., Difato, F., Massobrio, P., Breschi, G. L., Pasquale, V., & Levi, T., et al. (2013). In vitro large-scale experimental and theoretical studies for the realization of bi-directional brain-prostheses. Frontiers in Neural Circuits, 7(6), 40. Brown, P. (2003). Oscillatory nature of human basal ganglia activity: relationship to the pathophysiology of Parkinson's disease. Movement Disorders, 18(4), 357-363.

Beyeler, M., Oros, N., Dutt, N., & Krichmar, JL. (2015). A GPU-accelerated cortical neural network model for visually guided robot navigation. Neural Networks the Official Journal of the International Neural Network Society, 72, 75-87. Basu, A., Ramakrishnan, S., Petre, C., & Koziol, S. (2010). Neural dynamics in reconfigurable silicon. IEEE transactions on biomedical circuits and systems, 4(5), 311-319. Bose, A., Golowasch, J., Guan, Y., & Nadim, F. (2014). The role of linear and voltage-dependent ionic currents in the generation of slow wave oscillations. Journal of computational neuroscience, 37(2), 229-242. Cassidy, A. S., Georgiou, J., & Andreou, A. G. (2013). Design of silicon brains in the nano-CMOS era: Spiking neurons, learning synapses and neural architecture optimization. Neural Networks, 45(3), 4-26. Eliasmith, C., Stewart, T. C., Choo, X., Bekolay, T., DeWolf, T., Tang, Y., & Rasmussen, D. (2012). A large-scale model of the functioning brain, Science 338(6111), 1202-1205. Golubitsky, M., Stewart, I., Buono, P. L., & Collins, J. J. (1998). A modular network for legged locomotion. Physica D: Nonlinear Phenomena, 115(1-2), 56-72. Grossi, G., & Pedersini, F. (2008). FPGA implementation of a stochastic neural network for monotonic pseudo-Boolean optimization. Neural Networks, 21(6), 872-879. Ghosh-Dastidar, S., & Adeli, H. (2009). Spiking neural networks. International Journal of Neural Systems, 19(4), 295-308 Hamani, C,. Saint‐Cyr, J., A., Fraser, J., Kaplitt, M., & Lozano, AM. (2004). The subthalamic nucleus in the context of movement disorders. Brain, 127(1), 4-20. Holgado, A. J. N., Terry, J. R., & Bogacz, R. (2010). Conditions for the generation of beta oscillations in the subthalamic nucleus–globus pallidus network. The Journal of Neuroscience, 30(37), 12340-12352. Heer, F., Hafizovic, S., Ugniwenko, T., Frey, U., Franks, W., & Perriard, E., et al. (2007). Single-chip microelectronic system to interface with living cells. Biosensors and Bioelectronics, 22(11), 2546-2553. Indiveri, G., Linaresbarranco, B., Hamilton, T. J., Schaik, A. V., Etiennecummings, R., & Delbruck, T., et al. (2010). Neuromorphic silicon neuron circuits. Frontiers in neuroscience, 5(5), 73. Indiveri, G., Chicca, E., & Douglas, R. J. (2009). Artificial cognitive systems: from VLSI networks of spiking neurons to neuromorphic cognition. Cognitive Computation, 1(2), 119-127. Joucla, S., Ambroise, M., Levi, T., Lafon, T., Chauvet, P., & Saïghi, S., et al. (2016). Generation of Locomotor-Like Activity in the Isolated Rat Spinal Cord Using Intraspinal Electrical Microstimulation Driven by a Digital Neuromorphic CPG. Frontiers in Neurosciences, 10, 67. Levi, T., Lewis, N., Saïghi, S., Tomas, J., Bornat, Y., & Renaud, S. (2008). Neuromimetic integrated circuits. VLSI Circuits for Biomedical Applications, ed. I. Kris (Boston: Artech House), 241-264. Lu, J., Yang, J., Kim Y-B., & Ayers, J. (2012). Low power, high PVT variation tolerant central pattern generator design for a bio-hybrid micro robot. IEEE, International Midwest Symposium on Circuits and Systems IEEE, 59, 782–785. Krinskiĭ, V. I., & Kokoz, I. M. (1973). Analysis of the equations of excitable membranes. i. reduction of the hodgkins-huxley equations to a 2d order system. Biofizika, 18(3), 506.

Limousin, P., Krack, P., Pollak, P., Benazzouz, A., Ardouin, C., & Hoffmann, D., et al. (1998). Electrical stimulation of the subthalamic nucleus in advanced Parkinson's disease. New England Journal of Medicine, 339(16), 1105–1111. Magill, P J., Bolam, J. P., & Bevan, M. D. (2000). Relationship of activity in the subthalamic nucleus–globus pallidus network to cortical electroencephalogram. The Journal of neuroscience, 20(2), 820-833. Moxon, K. A., & Foffani, G. (2015). Brain-machine interfaces beyond neuroprosthetics. Neuron, 86(1), 55-67. Mallet, N., Pogosyan, A., Márton, L. F., Bolam, J. P., Brown, P., & Magill, P. J. (2008). Parkinsonian beta oscillations in the external globus pallidus and their relationship with subthalamic nucleus activity. The Journal of neuroscience, 28(52), 14245-14258. Matsumura, M., Kojima, J., Gardiner, T. W., & Hikosaka, O. (1992). Visual and oculomotor functions of monkey subthalamic nucleus. Journal of Neurophysiology, 67(6), 1615–1632. Nanami, T., & Kohno, T. (2016). Simple cortical and thalamic neuron models for digital arithmetic circuit implementation. Frontiers in Neuroscience, 10(181), 1--12. Nazari, S., Faez, K., Amiri, M., & Karami, E. (2015). A digital implementation of neuron-astrocyte interaction for neuromorphic applications. Neural Networks, 66(c), 79-90. Nazari, S., Amiri, M., Faez, K., & Amiri, M. (2015), Multiplier-less digital implementation of neuron–astrocyte signalling on FPGA. Neurocomputing, 164(c), 281-292. Nambu, A, Tokuno, H., Inase, M., & Takada, M. (1997). Corticosubthalamic input zones from forelimb representations of the dorsal and ventral divisions of the premotor cortex in the macaque monkey: comparison with the input zones from the primary motor cortex and the supplementary motor area. Neuroscience Letters, 239, 13-16. Parent, A., & Hazrati, L. N. (1995). Functional anatomy of the basal ganglia. I. The cortico-basal ganglia-thalamo-cortical loop. Brain Research Reviews, 20(1), 91-127. Pinto, C. M. (2007). Numerical simulations in two CPG models for bipedal locomotion. Journal of Vibration and Control, 13(9-10), 1487-1503. Piri, M., Amiri, M., & Amiri, M. (2015). A bio-inspired stimulator to desynchronize epileptic cortical population models: a digital implementation framework. Neural Networks, 67(c), 74-83. Plenz, D., & Kital, S. T. A basal ganglia pacemaker formed by the subthalamic nucleus and external globus pallidus. (1999). Nature, 400(6745), 677-682. Raikov, I., Preyer, A., & Butera, R. J. (2004). MRCI: a flexible real-time dynamic clamp system for electrophysiology experiments. Journal of neuroscience methods, 132(2), 109-123. Reato, D., Cammarota, M., Parra, L. C., & Carmignoto, G. (2012). Computational model of neuron-astrocyte interactions during focal seizure generation. Frontiers in Computational Neuroscience, 6(41), 81. Rubin, J. E., & Terman D. (2004). High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model. Journal of computational neuroscience, 16(3), 211-235. Rubchinsky, L. L., Park, C., & Worth, R. M. (2012). Intermittent neural synchronization in Parkinson’s disease. Nonlinear dynamics, 68(3), 329-346. Smith, Y., Bevan, M. D., Shink, E., & Bolam, J. P. (1998). Microcircuitry of the direct and indirect pathways of the basal ganglia. Neuroscience, 86(2), 353-388.

Thibeault, C. M., & Srinivasa, N. (2013). Using a hybrid neuron in physiologically inspired models of the basal ganglia. Frontiers in computational neuroscience, 7(7), 88. Terman, D., Rubin, J. E., Yew, A. C., & Wilson, C. J. (2002). Activity patterns in a model for the subthalamopallidal network of the basal ganglia. The Journal of neuroscience, 22(7), 2963-2976. Wichmann, T., Bergman, H. & DeLong, M. R. (1994). The primate subthalamic nucleus. I. Functional properties in intact animals. Journal of Neurophysiology, 72 (2), 494–506. Williams, S. R., & Mitchell, S. J. (2008). Direct measurement of somatic voltage clamp errors in central neurons. Nature neuroscience, 11(7), 790-798. Yu, Q., Tang, H., Tan, K. C., & Li, H. (2013). Rapid feedforward computation by temporal encoding and learning with spiking neurons. IEEE transactions on neural networks and learning systems, 24(10), 1539-1552. Yang, S., Wang, J., Li, S., Li, H., Wei, X., & Yu, H., et al. (2016). Digital implementations of thalamocortical neuron models and its application in thalamocortical control using FPGA for Parkinson‫ ׳‬s disease. Neurocomputing, 177, 274-289. Yang, S., Wang, J., Li, S., Deng, B., Wei, X., & Yu, H., et al. (2015). Cost-efficient FPGA implementation of basal ganglia and their Parkinsonian analysis. Neural Networks, 71(c), 62-75. Yang, S., Deng, B., Wang, J., Li, H., Liu, C., & Fietkiewicz, C., et al. (2017). Efficient implementation of a real-time estimation system for thalamocortical hidden Parkinsonian properties. Scientific Report,.9(7), 40152

Figure Captions: Fig. 1. Schematic representation for the connectivity of STN-GPe system. (a) A high level 3D plot of the STN-GPe system. (b) The main axis of the basal ganglia. The blue lines represent the GABAergic-inhibitory connections and the red line stands for the glutamatergic-excitatory connections. (b) The main axis of the basal ganglia (c) The equivalent circuit representation of a patch of the STN neuron membrane. Fig. 2. Modified functions in the neuron model. Blue lines represent the original functions. Red dots

are the knot points for the piecewise linear approximation and the segments connecting the red dots plotted by blue lines are the linear approximation segments for each linear function. Fig. 3. The phase portraits of the HOSTN model. (a) The V-n phase portrait. (b) The V-h phase portrait. (c) The V-r phase portrait. (d) The V-[Ca] phase portrait. (e) The phase plane of the original STN neuron model. (f) The phase plane of the HOSTN model. Fig. 4. The dynamics of the original and modified STN neuron models. The blue line and red line represent sodium ionic current and T-type calcium current respectively, and the high-threshold calcium current as well as the potassium ionic current are depicted by green and yellow lines. The black lines depict the leak current. (a) The steady-state current-voltage relation of the original model. (b) The steady-state current-voltage relation of the modified model. (c) The bifurcation analysis of the original model. (d) The bifurcation analysis of the HOSTN model. Fig. 5. The dynamical behaviors of the modified neuron models. (a) The theoretical analysis of the original model represented by the average spiking frequency and external applied current (f-Iapp). (b) The theoretical analysis of the HOSTN model represented by the average spiking frequency and external applied current (f-Iapp). (c) The dynamical response of the STN neuron model. (d) The dynamical response of the HOSTN model. The green lines represent the stimulation currents, and the black and blue lines represent the responses of the original and HOSTN models separately. Fig. 6. The dynamical behaviors and phase portraits of the coupled STN and GPe neurons with different values of synaptic conductance coefficient. The black and red lines represent the membrane potential of the GPe and STN nuclei respectively. (a) The synchronized biological activities of the coupled STN-GPe system. (b) The random phase relationship of the STN-GPe system. (c) The phase-locked biological activities of the coupled STN-GPe system.

Fig. 7. The phase portraits of the subthalamic dynamical behaviors. (a) The synchronization (a) The V-n phase portrait. (b) The V-h phase portrait. (c) The V-r phase portrait. (d) The V-[Ca] phase portrait. Fig. 8. Dependency of the firing frequency on the coupling strength gGPe→STN, gGPe→GPe, and gSTN→GPe. (a) Dependency of the firing frequency on the coupling strength gGPe→STN and gSTN→GPe for the GPe neuron. (b) Dependency of the firing frequency on the coupling strength gGPe→STN and gSTN→GPe for the STN neuron. (c) Dependency of the firing frequency on the coupling strength gGPe→GPe, and gSTN→GPe for the GPe neuron. (d) Dependency of the firing frequency on the coupling strength gGPe→GPe, and gSTN→GPe for the STN neuron. Fig. 9. The digital topology of the HOSTN neuron model. (a) General overview of pipelining structure of the HOSTN neuron model. (b) The V pipeline of the HOSTN model. (c) The n pipeline of the HOSTN model. (d) The h pipeline of the HOSTN model. (e) The r pipeline of the HOSTN model. (f) The [Ca] pipeline of the HOSTN model.

Fig. 10. The detailed digital topology of the HOSTN model. (a) The digital implementation of the "shift MUL" module. (b) The digital implementation of the "f1_block" module. Fig. 11. The proposed system architecture of the subthalamic nucleus–external globus pallidus pacemaker. (a) The pacemaker coupled by STN and GPe neurons. (b) The detailed digital structure of the GPe neuron. (c) The structured, sparsely connected architecture of the STN-GPe network. Fig. 12. The applications of the digital subthalamic nucleus–external globus pallidus pacemaker. (a) The digital implementation of the subthalamic modulator. (b) The central pattern generator composed by the proposed pacemaker. Different colors denote different coupling weights. Fig. 13. Biological behaviors of the proposed models and the subthalamic nucleus–external globus pallidus pacemakers. (a) The phase-locked mode. (b) The synchronization mode. (c) The

random-coupled mode. (d) The phase portrait of the phase-locked biological activities based on the STN-GPe system. (e) The phase portrait of the synchronized biological activities based on the STN-GPe system. The oscilloscope displays the patterns with x, 100 mV/div; y, 100 mV/div; time/div, 10ms. Fig. 14. Bipedal primary gaits based on the proposed CPG system. (a) Structure of the CPG network based on the STN-GPe pacemaker. Different colors denote different coupling weights. (b) The result of the proposed CPG system that is implemented in Altera Cyclone-Ⅳ EP4CE115 FPGA and shown on the oscilloscope via a digital-to-analog convertor. Fig. 15. The comparison between the outputs of the proposed CPG network with FPGA implementation and MATLAB simulation. The black lines and red lines stand for the CPG outputs based on software and FPGA respectively. The enlarged curves in blue dotted box are also depicted. Fig. 16. The error evaluation of the piecewise linear functions. (a) The CF RE of each piecewise linear function. (b) The CFNMAE of each piecewise linear function. (c) The color graph for the relative error of the multiplicand value with the bit width of the fractional part.

Table1

Functions

Values

Functions

Values





IL

g L (V  EL )

r V 

1 1  exp   V   r   r 

IK

g K n4 (V  EK )

m (V )

1 1  exp   V  m   m 

INa

g Na m3h(V  ENa )

a (V )

1 1  exp   V   a   a 

IT

gT a3b 2 (V  ECa )

b (r )

1 1  exp   r  b   b   1 1  exp  b  b  

ICa

gCa s 2 (V  ECa )

s (V )

1 1  exp   V   s   s 

 r (V )

 r0   r1 1  exp  V  r   r

Ca Ca  15

IAHP

g AHP (V  EK )

h V 

1 1  exp   V   h   h 

n V 

1 1  exp   V   n   n 









































 n (V )

 n0   n1 1  exp  V  n   n





 h (V )

 h0   h1 1  exp  V  h   h

Table2

Parameters

Values

Parameters

Values

α

2.0 msec-1

θg

20.0

β

0.08 msec-1

-57.0

E STN→GPe

0.0 mV

2.0

EGPe→GPe

-100.0 mV

g STN→GPe

0.3 nS/μm2

gGPe→GPe

0.1 nS/μm2

Table3

a1

a2

a3

a4

a5

a6

a7

a8

f1(V)

-3/256

-3/32

-11/32

-19/32

-1

-0.25

0.25

0.75

f2(V)

1.5

6.125

75/32

0

-31/32

--

--

--

f3(V)

3/256

23/16

39/8

13/8

-9/32

-31/32

--

--

g1(n)

0

3/64

3/16

51/128

--

--

--

--

g2(r)

0

1/16

7/16

20/16

--

--

--

--

g3(V)

3/4096

5/2048

3/512

9/1024

--

--

--

--

h1(V)

1/512

3/256

7/256

9/1024

1/4096

--

--

--

h2(V)

-1/4096

-1/64

-9/128

-1/64

0

--

--

--

h3(V)

-3/128

-13/128

-3/128

-6/4096 0

--

--

--

h4(V)

1/16384 7/512

9/128

1/64

1/8192

--

--

--

h5(V)

0

--

--

--

--

--

--

--

b1

b2

b3

b4

b5

b6

b7

b8

f1(V)

-0.8856

-5.709

-17.4

-27.3

-41.2

-37.35

-39.3

-46.18

f2(V)

112.3

411.5

245.95

161.8

138.56

--

--

--

f3(V)

0.9

81.1

238.2

164.3

138.785

138.56 --

--

g1(n)

0

-0.007

-0.05

-0.137

--

--

--

--

g2(r)

0

-0.001871 -0.0394 -0.165

--

--

--

--

g3(V)

0.076

0.148

0.19

0.153

--

--

--

--

h1(V)

0.147

0.690

1.37

1.02

0.991

--

--

--

h2(V)

0.985

0.18

-2.24

-0.4

0

--

--

--

h3(V)

-0.7757

-6.31

-1.374

-0.078

0

--

--

--

h4(V)

0.0064

0.709

3.2

1.4

0.999

--

--

--

h5(V)

21/512

--

--

--

--

--

--

--

Table4

Original model

HOSTN model

STN-GPe system

Resource

Avaliable

Total LEs

114480

16978(15%)

12590(11%)

18469(16%)

Total CFs

114480

12499(11%)

11216(10%)

16223(14%)

DLC

114480

7888(7%)

8433(7%)

12387(11%)

TMBs

3981312

289886(7%)

2632(<1%)

47507(<1%)

EM9Es

532

512(96%)

0(0%)

0(0%)

1(25%)

1(25%)

1(25%))

68.05MHz

125.38MHz

96.03MHz

Cost function

Original model

HOSTN model

Improvement

Neuron cost NC

1.1867

0.5983

49.58%

Delay D

2.1014

1.1405

45.73%

Resource cost H

27.2%

5.6%

79.41%

PLLs Frequency

4 143MHz

Resource Utilization

Table5

Gait

Left leg

Right leg

Two-legged hop

x1(t), x1(t)

x1(t), x1(t)

Walk

x1(t), x1(t+1/2)

x1(t+1/2), x1(t)

Two-legged jump

x1(t), x1(t+1/2)

x1(t), x1(t+1/2)

Run

x1(t), x1(t)

x1(t+1/2), x1(t+1/2)

Table6

Gait

Paramters

Two-legged hop

V1, V2, V3, V4 are not changed g=30

Walk

V1=V4=51 V2=V3=-20

Coupling weight matrix W 0.15 0.15 0.15  0  0.15 0 0.15 0.15   0.15 0.15 0 0.15    0.15  0.15  0.15 0   0.1 0.1 0.12  0  0.1 0 0.12 0.1   0.1 0.12 0 0.1   0.12  0.1  0.1 0  

g=30 Two-legged jump

V1=V3=51 V2=V4=-20

0.12 0.1 0.1  0 0.12 0 0.1 0.1   0.1 0.1 0 0.12   0   0.1 0.1 0.12

g=30 Run

V1=V4=51 V2=V3=-20 g=30

0.1 0.12 0.1  0  0.1 0 0.1 0.12  0.12 0.1 0 0.1    0.1  0.1 0.12 0  

Figure1

Figure2

Figure3

Figure4

Figure5

Figure6

Figure7

Figure8

Figure9

Figure10

Figure11

Figure12

Figure13

Figure14

Figure15

Figure16