2020 RISP International Workshop on Nonlinear Circuits,Communications and Signal Processing (NCSP 2020)
Honolulu, Hawaii, USA, February 28 - March 2, 2020

Modeling of Bursting Neurons and Its Characteristic using Nonlinear Time Series Analysis
Yoko Uwate† , Marie Engelene J. Obien‡ , Urs Frey‡ and Yoshifumi Nishio††‡

Dept. of E. E. Eng., Tokushima University
2-1 Minami-Josanjima Tokushima 770–8506 Japan
Phone: +81-88-656-7662
E-mail: {uwate, nishio}@ee.tokushima-u.ac.jp

AbstractMaxWell Biosystems
Albisriederstrasse 253, 8047 Zürich, Switzerland
Phone: +41-61-551-1070
E-mail: {marie.obien, urs.frey}@mxwbio.com

visualization method of network characteristics of burst patterns obtained from Izhikevich neuron model using nonlinear
time-series analysis. From the simulation results, we confirm
that the Izhikevich neuron model shows similar development
stages as biological neurons, by controlling the parameters of
the coupling probability and strength.

It is well known that burst patterns of neuronal networks may
play an important role in information processing in the brain.
We consider that it is advantageous to construct a model using mathematical neuronal models producing burst patterns,
because it is such models are easier to study and more accessible as compared to real biological neuronal data. In this
2. Attractor reconstruction of biological neuron data
study, we use the Izhikevich neuron model to produce burst
patterns and apply an attractor reconstruction to the Izhike2.1 Neuron data setvich neuron data.
1. Introduction
Burst patterns within neurons may play an important rolein information processing in the brain. Therefore, burst detection and burst pattern analyzes methods are developed,
which are used in various fields [1]. Although it is important
to study burst patterns, as these can relate to network synchronicity and synaptic connectivity, unveiling the structureof the whole neuronal network is also required. Nonlinear
time-series analysis is a useful tool for characterizing the dynamics behind observed time-series data [2]. Typically, neuronal data obtained from living neurons is high-dimensionaland dynamic data, therefore nonlinear time-series analysis issuitable to characterize neuronal data. Previously, we proposed a visualization method for bursting patterns of whole
neural networks using nonlinear time series analysis [3]. We
applied nonlinear time-series analysis to three sets of timeseries data from a neuronal culture measured at days in vitro(DIV) 15, 20 and 30, respectively. By using computer simulations, we confirmed that the characteristics of these neuronal
networks change across different DIV and that the proposed
visualization methods are suitable to capture the essence of
these changes.
However, it is difficult to extract the ground-truth of the underlying network morphology for data obtained from real biological neuronal networks. To relate these time-series analysis results to the underlying network morphology, it is required to construct a model using mathematical neuron models producing such burst patterns. In this study, we propose a

We used a CMOS-based high-density microelectrode array (HD-MEA) system MaxOne (MaxWell Biosystems AG,Zurich, Switzerland) composed of 26,400 platinum electrodes in a regular grid format with 17.5 m electrode centerto-center distance (3,265 electrodes/mm2)[4]. Up to 1,024
electrodes could be simultaneously recorded by routing theelectrodes to low-noise read-out channels through a flexible
switch-matrix approach [5]. On-chip circuitry was used toamplify (0-78 dB programmable gain), filter, and digitalize(10-bit, 20 kHz) the recorded signals. Online spike detection
was performed by the MaxLab Live software (threshold: x5RMS noise).
Primary cell cultures were prepared as described inRef. [6], in accordance to Swiss Federal Laws on animal welfare. Briefly, cells from embryonic day 18 Wistar rat cortices
were dissociated in 2 ml of trypsin with 0.25 % EDTA (Invitrogen, California, USA) with trituration. The MaxOne HDMEA recording setup was placed in a recording incubator(65 % humidity) to ensure controlled environmental conditions (36◦ C and 5 % CO2 ). During experiments, the MaxOne
chips were transferred to the recording incubator and covered
with sterilized lids to minimize media evaporation. Cultures
were recorded at DIV 15, 20, and 30.
2.2 Attractor reconstruction
Figure 1 shows the raster plot obtained from the same neuronal culture across multiple days; DIV 15, 20 and 30, using
1,024 electrodes. Next, we calculate the spike rate at timebins, then time-series data is obtained as shown in Fig. 2.

222

The attractor of dynamical systems can be reconstructed
topologically in the embedding space from Takens’ theorem [7]. The state vectors in the reconstructed m-dimensional
embedding space are defined byy(t) = {x(t), x(t + τ ), ..., x(t + (m − 1)τ )}(a) DIV 15.

(1)

where x(t) means a scalar time series and τ is the delay time.
Figure 3 shows the analysis results when neuronal timeseries data is embedded in 3-dimensional space with timedelay τ = 10. From this figure, we confirm that all threecultures exhibit a clear network structure, because the orbit
draws in a certain range and is not random. Furthermore, inthe case of the neurons at DIV 15 (Fig. 3(a)), the size of theattractor is larger and the behavior of the orbit is complex. As
the culture matures, the attractor size becomes smaller andthe complexity of the orbit becomes weaker.

(b) DIV 20.

(c) DIV 30.
"result_d15_t10.dat"

Figure 1: Raster plot of cultured neurons from 1,024 electrodes.

10000
90008000
7000
6000
5000
4000
3000
2000
1000
0

10000

0 1000
2000"Day15.dat"

9000

3000 4000

5000 6000

spike rate (Hz)

8000

7000 8000

9000100000

10000
9000
8000
7000
6000
5000
4000
3000
2000
1000

7000
6000

(a) DIV 15.

5000
4000
3000

"result_d20_t10.dat"

2000
1000
00102030t (s)

40

50

60
10000
9000
8000
7000
6000
5000
4000
3000
2000
1000
0

(a) DIV 15.
10000"Day20.dat"

9000

spike rate (Hz)

8000
7000

0 10002000 3000

6000

4000 5000

6000 7000

5000

8000 9000
100000

10000
9000
8000
7000
6000
5000
4000
3000
2000
1000

4000
3000

(b) DIV 20.

2000
1000
00102030t (s)

40

50

60

"result_d30_t10.dat"(b) DIV 20.
10000
9000
80007000
6000
5000
4000
3000
2000
1000
010000"Day30.dat"

9000

spike rate (Hz)

8000
7000
6000
5000
4000

0 10002000 3000

3000
2000

4000 5000

6000 7000

8000 9000
100000

10000
9000
8000
7000
6000
5000
4000
3000
2000
1000

1000
00102030t (s)405060(c) DIV 30.

Figure 2: Time series data of cultured neurons.

(c) DIV 30.

Figure 3: Attractor reconstruction (τ = 10). x-axis: x(t),y-axis: x(t + τ ), z-axis: x(t + 2τ ).

223

90000Instantaneous firing rate (sp/s)

3. Attractor reconstruction for Izhikevich neuron model
3.1 Izhikevich neuron modelIzhikevich neuron model efficiently produces a wide variety of neuron spiking and bursting dynamics. The Izhikevich
neuron model is described by the following equations.

80000
70000
60000
50000
40000
30000
20000
10000
001000

2000

3000

4000

5000

6000

5000

6000

t(ms)

(a) p=0.20, γ=1.6[mV].
Instantaneous firing rate (sp/s)

90000

dv
= 0.04v 2 + 5v + 140 − u + Iex
dtdu= a(bv − u)
dt

(2)

80000
70000
60000
50000
40000
30000
20000
10000
001000

2000

3000

4000

t(ms)

(b) p=0.08, γ=2.2[mV].

if v > 30 mV, then v ← c and u ← u + d.
Instantaneous firing rate (sp/s)

90000
80000
70000

where v represents the membrane potential of the neuron, urepresents a slow membrane recovery variable, accounting forthe activation of K+ ion currents and inactivation of Na+ ion
currents. Iex denotes the excitatory input current.
In the computer simulations, 1,000 neurons were coupled(c) p=0.03, γ=2.5[mV].
randomly. The ratio of excitability and inhibitory neuron is0.8 and 0.2, respectively. The parameters of the Izhikevich
neuron model are fixed with a=0.02, b=0.2, c =-50 and d=8. Figure 5: Time series data of the simulated neuronal network.
By using these parameters, the Izhikevich neuron model produces tonic spiking as shown in Fig. 4.
3.2 Attractor reconstruction
60000
50000
40000
30000
20000
10000
00

1000

2000

3000

4000

5000

6000t(ms)Figure 4: Tonic spiking of the neuron network model.

The coupling probability (p) of the network and the coupling strength (γ) of excitability neurons are important parameters for modeling the development of neurons. Figure 5shows the network activity of 1,000 neurons, resp. the firing
rate during a fixed bin, for different values of the coupling
probability and strength. Table 1 summarizes the three different parameter sets. Figure 5(a) uses a high coupling probability and weak coupling strength. While, Fig. 5(c) uses low
coupling probability and strong coupling strength.
Table 1: Parameterspγ [mV]
Figure 5 (a) 0.20
1.6Figure 5 (b) 0.08
2.2Figure 5 (c) 0.03
2.5

Finally, we apply the attractor reconstruction to neuron
data obtained from the simulated Izhikevich neuronal network model. The simulation results of the attractor reconstruction in a 3-dimensional space are shown in Fig. 6, where
the time delay is set to τ =3. Even though the mathematicalneuron model is used, a similar development process is observed as with the biological neuronal data.
Namely, we confirm that all three cultures exhibit a clear
structure, because the orbit draws in certain range and is not
random. Furthermore, in the case of the neurons of Fig. 5(a),the size of the attractor is larger and the behavior of the orbit is complex. The attractor size becomes smaller and thecomplexity of the orbit becomes weaker with a lower coupling probability and a stronger coupling strength, similar asobserved with the cultured neurons.
4. Conclusions
In this study, we applied the proposed technique, which isa visualization method of whole network characteristics forneuron data obtained from Izikevich neuron model. By using an attractor reconstruction technique, we obtained the attractors of actual measured neuronal data and also simulated
neuronal data in a 3-dimensional space. From the observed
results, we confirm that the Izhikevich neuron model, with224

appropriate parameters, may exhibit similar network properties as during the development of brain function in cultured
neuronal cells.
In future works, we would like to use neuron data with different parameters of the Izikevich neuron model. We are alsointerested to classify neuronal data using the techniques described here and apply such techniques also for applications
in drug development, such as observing and classifying effects of pharmacological compounds on neuronal networks.

90000
80000
70000
60000
50000
40000
30000
20000
10000
0References

010000
20000
30000
40000
50000
60000
70000
80000
900000

90000
80000
70000
60000
50000
40000
30000
20000
10000

(a) p=0.20, γ=1.6[mV].

[1] D. J. Bakkum, M. Radivojevic, U. Frey, F. Franke, A. Hierlemann and H. Takahashi, “Parameters for Burst Detection,”
Frontiers in Computational Neuroscience, doi: 10.3389/fncom.2013.00193, Jan. 2014.
[2] E. Bradley and H. Kantz, “Nonlinear Time-Series Analysis Revisited,” Chaos, vol. 25, doi: 10.1063/1.4917289, 2015.
[3] Y. Uwate, M. E. J. Obien, U. Frey and Y. Nishio, “Time Series
Analysis of Neurons and Visualization of Network Characteristics,” Proc. of NCSP19, 2019.

90000
80000
70000
60000
50000
40000
30000
20000
10000
0

[4] M. E. J. Obien, K. Deligkaris, T. Bullmann, D. J. Bakkum,
and U. Frey, “Revealing Neuronal Function through Microelectrode Array Recordings,” Frontiers in Neuroscience, 8(423).
http://doi.org/10.3389/fnins.2014.00423, 2015.
010000
20000
30000
40000
50000
60000
70000
80000
900000

90000
80000
70000
60000
50000
40000
30000
20000
10000

(b) p=0.08, γ=2.2[mV].

[5] U. Frey, J. Sedivy, F. Heer, R. Pedron, M. Ballini, J.
Mueller,
A. Hierlemann, “Switch-Matrix-Based HighDensity Microelectrode Array in CMOS Technology,”
Solid-State Circuits, IEEE Journal Of, 45(2), 467482.
http://doi.org/10.1109/JSSC.2009.2035196, 2010.
[6] D. J. Bakkum U. Frey, M. Radivojevic, T. L. Russell, J. Mller,
M. Fiscella, A. Hierlemann, “Tracking Axonal Action Potential Propagation on a High-Density Microelectrode Array
Across Hundreds of Sites,” Nature Communications, 4, 2181.
http://doi.org/10.1038/ncomms3181, 2013.

90000
80000
70000
60000
5000040000
30000
20000
10000
0010000
20000
30000
40000
50000
60000
70000
80000
900000

90000
80000
70000
60000
50000
40000
30000
20000
10000

(c) p=0.03, γ=2.5[mV].

[7] F. Takens, “Detecting Strange Attractors in Turbulence,” In D.A.
Rand and B.-S. Young, editors, Dynamical Systems and Turbulence, vol. 898 of Lecture Notes in Mathematics, pp. 366381,
Warwick 1980, 1981. Springer-Verlag, Berlin.
[8] J. P. Eckmann, S. O. Kamphorst and D. Ruelle, “Recurrence
Plots of Dynamical Systems,” Europhys. Lett., vol. 4, pp. 973,
1987.
[9] M. Thiel, M. C. Romano and J. Kurths “How Much Informationis Contained in a Recurrence Plot?,” Phys. Lett. A, vol. 330, pp.
343–349, 2004.

Figure 6: Attractor reconstruction of the simulated neuronalnetwork(τ = 3). x-axis: x(t), y-axis: x(t + τ ), z-axis: x(t +[10] N. Marwan, M. C. Romano, M. Thiel and J. Kurths, “Recur2τ ).

rence Plots for the Analysis of Complex Systems,” Physics Reports, vol. 438, pp. 237–329, 2007.

[11] A. Wolf, J. B. Swift, H. L. Swinney and J. A. Vastano, “Determining Lyapunov Exponents from a Time Series,” Physica D,
vol. 16, pp. 285–317, 1985.

225