Journal of Robotics, Networking and Artificial Life, Vol. 5, No. 1 (June 2018) 32–36___________________________________________________________________________________________________________A Metaheuristic Approach for Parameter Fitting in Digital Spiking Silicon Neuron ModelTakuya Nanami

The University of Tokyo, Institute of industrial ScienceTokyo, Japan
Filippo Grassia

LTI Lab., University of Picardie Jules VerneSaint-Quentin, FranceTakashi Kohno

The University of Tokyo, Institute of industrial ScienceTokyo, Japannanami@sat.t.u-tokyo.ac.jp, filippo.grassia@u-picardie.fr, kohno@sat.t.u-tokyo.ac.jp

Abstract
DSSN model is a qualitative neuronal model designed for efficient implementation in digital arithmetic circuit. In
our previous studies, we developed automatic parameter fitting method using the differential evolution algorithm
for regular and fast spiking neuron classes. In this work, we extended the method to cover low-threshold spiking
and intrinsically bursting. We optimized parameters of the DSSN model in order to reproduce the reference ionicconductance model.
Keywords: Spiking neuron model, Low-threshold spiking, Intrinsically bursting, Differential evolution, FPGA

fire (I&F) based model which approximates a spiking
process by the resetting of the state variables. The I&Fbased models can be implemented resource-efficiently,
but are capable of reduced reproducibility of complex
neuronal activities. For example, these models assume
that the amplitude of spikes is fixed, whereas the spike
intensity is nonuniform in the nervous system. 4
We have been studying qualitative neuronal models
for digital as well as analog circuit implementation thatsatisfy both the reproducibility of neuronal activities
and low computational cost. 5, 6, 7, 8 The core idea of ourqualitative-modeling-based approach is to reproduce thecore mathematical structures that a wide variety ofneuronal activities. In our previous studies, 9, 10 weextended the DSSN models7 to support various neuronalclasses; regular spiking (RS), fast spiking (FS),intrinsically bursting (IB), low-threshold spiking (LTS),

1. Introduction
The silicon neuronal network is composed of analog and(or) digital circuits that solve the differential equationsof neuronal and synaptic models. Different kinds ofneuronal models have been used in silicon neuronal
networks because of the trade-off between
reproducibility of neuronal activities and computational
efficiency. For example, Merolla et al.1 constructed adigital silicon neuronal network comprising 1 million
spiking neurons and 256 million synapses. Theyadopted Leaky integrate-and-fire (LIF) model that is
one of the most commonly used and simple neuronmodels. Thomas et al.2 implemented 1024-spiking
neurons of Izhikevich (IZH) model3 on a fieldprogrammable gate array (FPGA). Those LIF, expanded
LIF, and IZH models are categorized as integrate-and-Copyright © 2018, the Authors. Published by Atlantis Press.
This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).
32

Journal of Robotics, Networking and Artificial Life, Vol. 5, No. 1 (June 2018) 32–36___________________________________________________________________________________________________________elliptic bursting (EB), and parabolic bursting (PB). We
also studied an automatic parameter fitting method
using the differential evolution (DE) algorithm11 for theRS and FS classes. 12
In this work, we developed the method for the lowthreshold spiking and intrinsically bursting classes.
Parameters were optimized to reproduce waveforms of
the ionic-conductance models in Ref. 13 and reduce thecircuit resource requirements for implementation.
The remainder of this paper is organized as follows.
Section 2 introduces our neuron model, the DE
algorithm, and details of our parameter optimization
procedure. The result is shown in Section 3. Section 4summarizes this work and suggests ideas for the future.

constant andrepresents the input stimulus.
Parameters ∅ , ∅ , ∅ , ε , , and τ control the time
constants of the variables. Parameters, , , ,,, and , where is , ,,, , or , areconstants that adjust the shape of the -, -, - and nullclines. Parameters
andare used to select ∅depending on . All of the variables and constants inthis qualitative model are purely with no physical unit.
Cubic shaped -nullcline is a key to replicate thespiking dynamics. Because multipliers are resourceconsuming in a digital arithmetic circuit, the DSSN
model adopts a piecewise quadratic function for the nullcline so that its numerical integration step requires
only one multiplication between variables.

2. Method

2.2. Differential evolution algorithm
The DE algorithm11 is a popular heuristic method tosolve optimization problems using a real numberfunction. It is one of the population-based optimization
algorithms in which a number of parameter sets areconverged to quasi-optimal solutions with which theDSSN model reproduces the target waveform. The
initial random parameter sets have to be selected
carefully so that the behavior of the DSSN model does
not diverge. Detailed explanation is written in theprevious study. 12

2.1. Digital Spiking Silicon Neuron model
The 4-variable DSSN model is a qualitative neuronmodel that can simulate various neuronal classes (RS,FS, IB, LTS, EB, and PB) by Euler's method with fixedpoint operation. Its equations are given by∅(1),1,(2),(3)

∅ u∅
∅∅

≡

The LTS class neurons are found in the thalamus. It has
the spike-frequency adaptation as with the RS class, andhas the rebound bursting; neurons in the LTS classgenerate a burst firing just after the termination of a
sufficiently strong and long hyperpolarizing stimulus.
The IB class is a well-known bursting neuron class inthe cortex. Neurons in the IB class generate a burst
firing immediately at the onset of a sufficiently strong
stimulus, then continue spiking until its termination.
The DSSN model has 28 parameters, and exploring
their spanning space is not realistic for the computerswith average power. To avoid this problem, we divided
the state variables of the DSSN model into two groups
based on their time scale, and tuned each group insequence. In our previous study, 12 the RS and FS
classes were characterized as the degree of spikefrequency adaptation, and we just optimized parametersof fast state variables in order to reproduce the saddle-(5),00 ,≡≡2.3. Parameter optimization procedure(4),(6)

c,(7)

c(8)
,
wherecorresponds to the membrane potential.
Variables and are the fast and slow variables which
abstractly describe the activity of the fast and slow ionchannels, respectively. Variable
is the slowest
variable that modifies the structure of the fast subsystemcomprising v and n, and used in only the IB mode. In
the RS, FS, and LTS modes, is fixed. Parameter is a33

Journal of Robotics, Networking and Artificial Life, Vol. 5, No. 1 (June 2018) 32–36___________________________________________________________________________________________________________Fig. 1. Waveforms of the ionic-conductance model (red) andthe DSSN model (blue) of the LTS class in response to aninhibitory (top) and weak excitatory (middle) and strong
excitatory (bottom) step stimulus input.

Fig. 2. Waveforms of the ionic-conductance model (red) andthe DSSN model (blue) of the IB class in response to a weak(top) and middle (middle) and strong (bottom) excitatory step
stimulus input.

node on invariant circle bifurcation. However, LTS andIB classes have more complex structure in their fasts
subsystem. For the LTS class, we firstly explored theparameter group responsible for the dynamics of the faststate variables and so that the fast subsystem of ourmodel is consistent with that of the ionic conductancemodel of the LTS class. Then, the parameters related tothe slow state variable were determined so that the
dynamical behavior of the whole ionic conductancemodel of the LTS class is reproduced. For the IB class,we firstly fitted the parameters responsible for thedynamics of the state variables , , and so that the
dynamical behavior of the corresponding subsystem in
the ionic conductance model of the IB class isreproduced. Then, the parameters related to theremaining variable
were determined so that thedynamical behavior of the whole ionic conductancemodel of the IB class is reproduced.

In the DSSN model’s circuit implementation,
multiplication between a coefficient and a state variable
was realized by shifters and adders, and we needed alarger number of adders and shifters as the number ofdigits with value 1 in the fixed point representation of
the coefficient increases. Parameters were optimized for
not only reproducing waveforms of the ionicconductance model but also reducing the circuit size.
We defined the cost functions and as follows,
(9)
,
(10)
,whereis the function that calculates cost from thetarget waveforms of the ionic conductance model andthe DSSN model with parameter vector , based on thedifference of spike timing as well as the mean squareerror of the state variable corresponding the membranepotential. In this function, minimum and maximum
values of waveforms are normalized to 0 and 1 for34

Journal of Robotics, Networking and Artificial Life, Vol. 5, No. 1 (June 2018) 32–36___________________________________________________________________________________________________________comparison. Function
is a number of digits withvalue 1 in the binary fixed-point expression of thecoefficients in . Constantis to balance these twocost functions. Parameters responsible for the dynamics
of the fast variables were determined to minimize thecost function , using the DE algorithm.
In the LTS class, parameters related to state variables
and were firstly determined to minimize the costfunction . Then the slow state variables are fixed in
the ionic conductance model, which is used as the target
waveforms in the. Secondly, the remaining
parameters related to the slow variable
weredetermined to minimize the cost function. The
original ionic conductance model is used in the. Inthe IB class, parameters related to state variables , ,and were firstly determined to minimize the costfunction . Secondly, the remaining parameters related
to the slowest variable were determined to minimize
the cost function .

Fig. 3. Statistical properties of the ionic-conductance model(red) and the DSSN model (blue). In the LTS mode (left), thex−axis and y-axis are
and . In the IB mode (right), xaxis corresponds to the number of spikes in the bursting and yaxis corresponds firing rate after the bursting.

3. Result
We found two parameter sets with which the DSSNmodel reproduces the activities of the ionic-conductancemodels for LTS and IB classes (Fig. 1). We evaluatedthe similarity of the spiking patterns between the ionicconductance model and the DSSN model by calculatingsome statistics while changing the stimulus intensity(Fig. 3). For the LTS class, we usedand 14 andthese are explained in our previous study. 12 Our fitted
model’scharacteristics are not completely thesame as that of the reference model, but share the trend
thatanddecreases as the stimulus intensity
increases. For the IB class, we measured the number ofspikes in the bursting duration and average firing rate
after the bursting. With both models, the number of
spikes in the bursting duration increases linearly as theaverage firing rate increases, although the exact value
does not coincide completely. The similarity of thesecharacteristics will be further estimated based on therequirements specified by the circuits' applications.
The results in Figs. 1, 2, and 3 were obtained bynumerical integration with fixed-point operation and theEuler's method (2s) in C language. The same
results are expected to be observed in the following
FPGA implementation, because this condition is thesame as our VHDL codes. We compiled the DSSN
models for Xilinx XC7VX690T using Xilinx Vivado

Table. 1 : Device Utilization.

Design Suite. Device utilization is listed in Table 1. In
the column labeled LTS (hand-fitted) and IB (handfitted), the resource usage for LTS and IB settings inRef. 9 is shown. The requirement for LUTs was reduced
without any loss of the reproducibility of the neuronalactivities.
4. Conclusion
In this work, we extended the parameter optimization
method for the DSSN model in the LTS and IB classmodes. These classes have complex mathematicalstructures of the fast subsystem, thus we carefully
designed our fitting procedure in order to reproduce
important structures in these classes. We could findparameter sets with which the DSSN model reproduces
the characteristics activities of LTS and IB classes using35

Journal of Robotics, Networking and Artificial Life, Vol. 5, No. 1 (June 2018) 32–36___________________________________________________________________________________________________________average-power PC ( a Xeon E5-2667 CPU with a
NVIDIA P100 GPU ). Cost functions
wereintroduced to reduce the necessary circuit size. We
confirmed that the DSSN model with determinedparameter set can precisely reproduce the ionic
conductance model in C simulation, and the necessary
circuit size was simulated by Xilinx Vivado Design
Suite. In our future work, we will apply our siliconneuronal networks to the spike pattern recognition. 15

10.

11.

12.

Acknowledgement
13.

This work was partially supported by JSPS SAKURA
Program and JST PRESTO and CREST

References
14.

1. P. A. Merolla, J. V Arthur, R. Alvarez-Icaza, A. S.
Cassidy, J. Sawada, F. Akopyan, B. L. Jackson, N. Imam,C. Guo, Y. Nakamura, B. Brezzo, I. Vo, S. K. Esser, R.
Appuswamy, B. Taba, A. Amir, M. D. Flickner, W. P.
Risk, R. Manohar, and D. S. Modha, “A million spikingneuron integrated circuit with a scalable communication
network and interface”, Science, Vol. 345, pp. 668–673,
2014
2. D. B. Thomas and W. Luk, “FPGA accelerated
simulation of biologically plausible spiking neuralnetworks,” Proc. - IEEE Symp. F. Program. Cust.
Comput. Mach. FCCM 2009, pp. 45–52, 2009.
3. E. M. Izhikevich, “Which model to use for corticalspiking neurons ?”, IEEE Trans. Neural Netw. vol. 15, pp.
1063–1070, 2004
4. H. Alle and J. R. P. Geiger, “Combined analog and action
potential coding in hippocampal mossy fibers”, Science,
vol. 311, pp. 1290–1293, 2006
5. T. Kohno and K. Aihara “A three-variable ultralowpower analog silicon neuron circuit” 2016 International
Symposium on Nonlinear Theory and Its Applications,Japan, November 27th–30th, pp. 190–193, 2016
6. T. Kohno, M. Sekikawa, J. Li, T. Nanami, and K. Aihara,"Qualitative-Modeling-Based Silicon Neurons and Their
Networks," Frontiers in Neuroscience, Vol. 10, No.273,
pp. 1–16, Jun., 2016.
7. T. Kohno and K. Aihara. Digital spiking silicon neuron:
Concept and behaviors in gj-coupled network.
Proceedings of International Symposium on Artificial
Life and Robotics 2007, 2007.
8. J. Li, Y. Katori, and T. Kohno, "An FPGA-based siliconneuronal network with selectable excitability siliconneurons," Frontiers in Neuroscience, Vol. 6, pp. 1–13,Dec., 2012.
9. T. Nanami and T. Kohno. “Simple cortical and thalamic
neuron models for digital arithmetic circuit

15.

36implementation”, Frontiers in Neuroscience, Vol. 10,No.181, pp. 1–12, 2016
T. Nanami, K. Aihara, and T. Kohno, “Elliptic andparabolic bursting in a digital silicon neuron model”,
2016 International Symposium on Nonlinear Theory andIts Applications, Japan, November 27th–30th, 2016
R. Storn, K. Price, “A Simple and Efficient Heuristic for
global Optimization over Continuous Spaces”, Journal ofGlobal Optimization, Vol. 11, pp. 341–359, 1997.
T. Nanami, F. Grassia, and T. Kohno, “A parameter
optimization method for Digital Spiking Silicon Neuronmodel”, Journal of Robotics, Networking and Artificial
Life, Vol. 4, pp. 97–101, 2017.
M. Pospischil, M. Toledo-Rodriguez, C. Monier, Z.
Piwkowska, T. Bal, Y. Fregnac, H. Markram, and A.
Destexhe. “Minimal Hodgkin-Huxley type models fordifferent classes of cortical and thalamic neurons”.
Biological Cybernetics, Vol. 99, pp. 427–441, 2008.
S. Shinomoto, K. Shima, and J. Tanji, “Differences inspiking patterns among cortical neurons.,” NeuralComput., vol. 15, no. 12, pp. 2823–2842, 2003.
T. Masquelier, R. Guyonneau, and S. J. Thorpe,“Competitive STDP-based spike pattern learning.,”
Neural Comput., vol. 21, pp. 1259–1276, 2009.