2019 RISP International Workshop on Nonlinear Circuits, Communications and Signal Processing (NCSP2019)
Honolulu, Hawaii, USA, March 4-7, 2019

Time Series Analysis of Neurons and Visualization of Network Characteristics
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

MaxWell Biosystems
Mattenstrasse 26, 4058 Basel, Switzerland
Phone: +41-61-551-1070
E-mail: {marie.obien, urs.frey}@mxwbio.com

Abstract

respectively. By using computer simulations, we confirm thatthe characteristics of these neuronal networks change acrossIn this study, we use nonlinear time-series analysis applied to different DIV.
neuronal data for visualizing specific characteristics of neuronal activity. We used extracellular spike recordings of a
2. Neuron data set
neuronal cell culture, observed at three different days in vitro.
By applying three nonlinear time-series analysis, we confirmed that the neuronal culture has the highest activity during 2.1 Measurement platform
early network development and the neuronal activity settles
We used a CMOS-based high-density microelectrode array
down as the culture becomes mature.
(HD-MEA) system called MaxOne (MaxWell Biosystems,
1. Introduction
Understanding how brain circuits develop and operate is amajor goal for many neuroscience projects. Burst patterns inneuronal networks may have an important role in information
processing in the brain. Therefore, detecting and analyzing
burst patterns are investigated in various fields [1]. Although
it is important to study burst patterns in order to understandthe correlation and communication processes of neurons, unveiling a structure of the whole neuronal network is also required.
Nonlinear time-series analysis is a useful tool for characterizing the dynamics behind the observed time-series data [2].
The neuronal data obtained from living neurons should behigh-dimensional and of dynamic nature. In such a case,nonlinear time-series analysis can be used to characterize theneuronal data. Here, three techniques are applied: attractor
reconstruction, recurrence plot and Lyapunov exponent.
In this study, we focus on time-series data generated fromspike time data of neurons in culture. The spike times are detected using a high-density microelectrode array (HD-MEA)
system [3] called MaxOne (MaxWell Biosystems, Basel,Switzerland). The system features 26,400 electrodes in total and a sampling rate of 20 kHz. We chose a recording timeof 60 seconds. First, we detect the spikes per electrode andvisualize the spike times in a raster plot. Next, we convert
the spike time data to time-series data by calculating the spikerate of all electrodes per time bin. Finally, we apply nonlinear
time-series analysis to the three time-series data from a neuronal culture measured at days in vitro (DIV) 15, 20 and 30,Basel, Switzerland) composed of 26,400 platinum electrodes
in a regular grid format with 17.5 m electrode center-to-center
distance (3,265 electrodes/mm2). Up to 1,024 electrodes
could be simultaneously recorded by routing the electrodes tolow-noise read-out channels through a flexible switch-matrix
approach [4]. On-chip circuitry was used to amplify (0-78 dB
programmable gain), filter, and digitalize (10-bit, 20 kHz) therecorded signals. Online spike detection was performed by
the MaxLab Live software (threshold: x5 RMS noise).
2.2 Cortical cell cultures
Primary cell cultures were prepared as described inRef. [5], 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 electrode array surface was pre-coated with a thinlayer of poly (ethyleneimine) (Sigma, Missouri, USA),0.05 % by weight in borate buffer (Chemie Brunschwig,Basel, Switzerland) at 8.5 pH, followed by a 10 ml drop of0.02 mg/ml laminin (Sigma) in Neurobasal medium (Invitrogen, California, USA) for cell adhesion. 20000-30000 cells
in an 8-l drop were seeded over the array, and 1 ml of platingmedium was added after 20 minutes. After 24 hours, the plating medium was changed to growth medium. Plating medium
consisted of 850 ml of Neurobasal, supplemented with 10 %horse serum (HyClone, Utah, USA), 0.5 mM GlutaMAX (Invitrogen, California, USA) and 2 % B27 (Invitrogen, California, USA). Growth media consisted of 850 ml of DMEM Dulbecco’s Modified Eagle Medium (Invitrogen, California,450

2019 RISP International Workshop on Nonlinear Circuits, Communications and Signal Processing (NCSP2019)
Honolulu, Hawaii, USA, March 4-7, 2019

10000USA), supplemented with 10 % horse serum, 0.5 mM GlutaMAX and 1 mM sodium pyruvate (Invitrogen, California,USA). The cultures were maintained inside an incubator tocontrol environmental conditions (37◦ C, 65 % humidity, 9 %O2, 5 % CO2) in 1 ml of growth medium (partially replaced
twice per week).

"Day15.dat"

9000

spike rate (Hz)

8000
7000
6000
5000
4000
3000
2000
1000
00102030t (s)405060(a) DIV 15.

2.3 HD-MEA extracellular recordings
10000

The MaxOne HD-MEA recording setup was placed in arecording 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.
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.

"Day20.dat"

9000

spike rate (Hz)

8000
7000
6000
5000
4000
3000
2000
1000
00102030t (s)405060(b) DIV 20.
10000"Day30.dat"

9000

spike rate (Hz)

8000
7000
6000
5000
4000
3000
2000
1000
00102030t (s)405060(c) DIV 30.

Figure 2: Time series data of the neurons.
(a) DIV 15.

3. Simulation results using nonlinear time-series analysis
In this section, we explain three techniques of nonlinear
time-series analysis and show the results.
3.1 Attractor reconstruction(b) DIV 20.

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

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

(1)

where x(t) means a scalar time series and τ is the delay time.
Figure 3 shows the simulation results when neuronal timeseries data is embedded in 3-dimensional space with time delay τ = 10. From this figure, we confirm that all three cultures exhibit a clear structure, because the orbit draws in certain range and does not move about randomly. Furthermore,in the case of the neurons at DIV 15 (Fig. 3(a)), the size ofthe attractor 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.

4512019 RISP International Workshop on Nonlinear Circuits, Communications and Signal Processing (NCSP2019)
Honolulu, Hawaii, USA, March 4-7, 2019

3.2 Recurrence plot
In 1987, Eckmann et al. introduced the method of recurrence plots (RPs) to visualize the recurrences of dynamicalsystems [7]. RPs shows when the time series visits the same
region on the phase space. The recurrence of states in thephase space by a 2-dimensional plot is described by the following equation.

"result_d15_t10.dat"

10000
9000
8000
7000
6000
5000
4000
3000
2000
1000
0R(i, j) =0 1000
2000

3000 4000

5000 6000

7000 8000

9000100000

10000
9000
8000
7000
6000
5000
4000
3000
2000
1000

(1 (||~x(i) − ~x(j)|| < ε)0 (otherwise)(2)

A recurrence is then defined when the distance between
two states i and j (points on the trajectory) is smaller than athreshold ε. RP is a N × N matrix of black and white dotsin two time-axes. Using RPs for the analysis of time series,
allows not only to visualize, but also to quantify the structures
hidden in the data [8], [9].
The simulation results of RPs are shown in Fig. 4. The
number of recurrence points increases with days. First twocases (Fig. 4(a) and (b)) show that the pattern of RPs has no
regular form. In the last case (Fig. 4(c)), the pattern of RPs
becomes more regular.

(a) DIV 15.
"result_d20_t10.dat"

10000
9000
8000
7000
6000
5000
4000
3000
2000
1000
0

3.3 Lyapunov exponent

0 10002000

3000 4000

5000 6000

7000 8000

9000100000

10000
9000
8000
7000
6000
5000
4000
3000
2000
1000

Lyapunov exponents provide a qualitative and quantitative
characterization of the dynamical behavior. It is related to theexponentially fast divergence or convergence of nearby orbits
in phase space. The Lyapunov exponent is used to determine
whether a time series has a chaotic property in order to determine the sensitivity of the chaotic route to its initial conditions. The Lyapunov exponent is calculated by the following
equation [10]. It is well known that the orbit is unstable andchaotic when λ has a positive value.

(b) DIV 20.
"result_d30_t10.dat"

n−1

10000
9000
8000
7000
6000
5000
4000
3000
2000
1000
01X
∆x0
log
x→∞ n∆xni=1

λ = lim

0 10002000 3000

4000 5000

6000 7000

8000 9000
100000(3)

The calculated Lyapunov exponents are summarized inTable 1. The Lyapunov exponents of all three networks
have positive value. Namely, these neuronal time-series has
chaotic property. The youngest culture (DIV 15) has thelargest value of the Lyapunov exponent, and the value decreases as the culture matures.

10000
9000
8000
7000
6000
5000
4000
3000
2000
1000

(c) DIV 30.

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

452

Table 1: Lyapunov exponent.
neuron typeλDIV 153.82DIV 20
3.20
DIV 302.54

2019 RISP International Workshop on Nonlinear Circuits, Communications and Signal Processing (NCSP2019)
Honolulu, Hawaii, USA, March 4-7, 2019

3000

References

2500[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.

time

2000

1500

1000

[2] E. Bradley and H. Kantz, “Nonlinear Time-Series Analysis Revisited,” Chaos, vol. 25, doi: 10.1063/1.4917289,
2015.

500

005001000

1500
time

2000

2500

3000

[3] 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.

(a) DIV 15.
3000

2500

time

2000

[4] 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.

1500

100050000500

1000

1500

2000

2500

3000time[5] D. J. Bakkum U. Frey, M. Radivojevic, T. L.
Russell, J. Mller, M. Fiscella,A. Hierlemann,“Tracking Axonal Action Potential Propagation ona High-Density Microelectrode Array Across Hundreds of Sites,” Nature Communications, 4, 2181.
http://doi.org/10.1038/ncomms3181, 2013.

(b) DIV 20.
3000

2500

time

2000

1500

[6] 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. SpringerVerlag, Berlin.

100050000500

1000

1500

2000

2500

3000time(c) DIV 30.

[7] J. P. Eckmann, S. O. Kamphorst and D. Ruelle, “Recurrence Plots of Dynamical Systems,” Europhys. Lett., vol.
4, pp. 973, 1987.

Figure 4: Recurrence plot (ε = 1.0).

[8] M. Thiel, M. C. Romano and J. Kurths “How Much Information is Contained in a Recurrence Plot?,” Phys. Lett. A,
vol. 330, pp. 343–349, 2004.

4. Conclusions[9] N. Marwan, M. C. Romano, M. Thiel and J. Kurths, “Recurrence Plots for the Analysis of Complex Systems,”
Physics Reports, vol. 438, pp. 237–329, 2007.
In this study, we propose a visualization method for neuronal data using nonlinear time-series analysis. We used threerecordings from neurons observed at different ages. By ap- [10] A. Wolf, J. B. Swift, H. L. Swinney and J. A. Vastano,“Determining Lyapunov Exponents from a Time Series,”
plying three nonlinear time-series analysis, we confirmed thatPhysica D, vol. 16, pp. 285–317, 1985.
the youngest neuron has the strongest activity and the neuronal behavior settles down as the culture matures.
In future works, we would like to calculate the characteristics of RPs, and to use data from different types of cells andconditions of the neuronal cultures. We are also interested toclassify neuronal data using the techniques described here.

453