Computational modeling of spiking neural network with learning rules from STDP and intrinsic plasticity

Accepted Manuscript Computational modeling of spiking neural network with learning rules from STDP and intrinsic plasticity Xiumin Li, Wei Wang, Fangzheng Xue, Yongduan Song

PII: DOI: Reference:

S0378-4371(17)30789-6 http://dx.doi.org/10.1016/j.physa.2017.08.053 PHYSA 18499

To appear in:

Physica A

Received date : 21 April 2017 Revised date : 22 July 2017 Please cite this article as: X. Li, W. Wang, F. Xue, Y. Song, Computational modeling of spiking neural network with learning rules from STDP and intrinsic plasticity, Physica A (2017), http://dx.doi.org/10.1016/j.physa.2017.08.053 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.

Computational modeling of spiking neural network with learning rules from STDP and intrinsic plasticity 1, 2, ∗ Xiumin Li, 1, 2 Wei Wang,1, 2 Fangzheng Xue, 1, 2 Yongduan Song. 1 Key Laboratory of Dependable Service Computing in Cyber Physical Society of Ministry of Education, Chongqing University, Chongqing 400044, China. 2 College of Automation, Chongqing University,Chongqing 400044, China. ∗ [email protected]

Abstract Recently there has been continuously increasing interest in building up computational models of spiking neural networks (SNN), such as the Liquid State Machine (LSM). The biologically inspired self-organized neural networks with neural plasticity can enhance the capability of computational performance, with the characteristic features of dynamical memory and recurrent connection cycles which distinguish them from the more widely used feedforward neural networks. Despite a variety of computational models for brain-like learning and information processing have been proposed, the modelling of self-organized neural networks with multi- neural plasticity is still an important open challenge. The main difficulties lie in the interplay among different forms of neural plasticity rules and understanding how structures and dynamics of neural networks shape the computational performance. In this paper, we propose a novel approach to develop the models of LSM with a biologically inspired self-organizing network based on two neural plasticity learning rules. The connectivity among excitatory neurons is adapted by spike-timing-dependent plasticity (STDP) learning; meanwhile, the degrees of neuronal excitability are regulated to maintain a moderate average activity level by another learning rule: intrinsic plasticity (IP). Our study shows that LSM with STDP+IP performs better than LSM with a random SNN or SNN obtained by STDP alone. The noticeable improvement with the proposed method is due to the better reflected competition among different neurons in the developed SNN model, as well as the more effectively Preprint submitted to Nuclear Physics B

September 11, 2017

encoded and processed relevant dynamic information with its learning and self-organizing mechanism. This result gives insights to the optimization of computational models of spiking neural networks with neural plasticity. Keywords: STDP, Intrinsic Plasticity, Spiking Neural Network, reservoir computing 1. Introduction With the rapid development of theories and applications in brain science, spiking neural network (SNN) have very solid theoretical and experimental basis and show more powerful calculation and associative memory ability than traditional artificial neural network. Many recent advances have been made in modeling learning and computation with recurrently connected spiking neurons capable of performing a wide variety of tasks, such as working memory [1, 2], decision making [3] and predictive coding [4], et al (for a review, see [5]). However, the brain-inspired intelligence with computational application of SNN is still an important open challenge. The main difficulties lie in the configuration of model parameters and training of the massive synaptic weights. Liquid state machine (LSM) represents one of the most efficient computing models of SNNs proposed by Maass et al. [6]. LSM exploits the recurrent SNNs without training the entire network online. The liquid component (LC; or liquid) in LSM uses spiking neurons connected by synapses to project inputs into a high-dimensional feature space called ”liquid state”. For specific tasks, the liquid state can be mapped onto a target signal through a readout component (RC), which acts as a memory-less readout function [6]. Instead of updating all of the weights online, like in the above mentioned recurrent neural networks, synaptic connections in LC are usually selected randomly and fixed during the training process; only the RC is trained using a simple classification/regression technique according to specific tasks. This kind of network design, including the echo state networks (ESNs) [7], is often referred to as ”reservoir computing” [8, 9]. The key problem of reservoir computing is the design of reservoir or neural network, which can enhance computational performance and reduce computational complexities. Various methods for constructing topologies of neural networks have been developed [10]. Studies have shown that the optimization of SNN can enhance the computational accuracy of LSMs for a specific task [11, 12, 13, 14]. 2

Instead of a prior imposition of a specific topology for LSM, considering selforganizing neural networks based on neural plasticity is more natural. Neural plasticity is the ability of the brain to change synaptic strength in response to external stimuli and is an important phenomenon for neuroscience. One of the best-known forms of synaptic plasticity is the Hebbian learning or spike-timing-dependent plasticity (STDP) [15], which shows an asymmetric time window for cooperating and competing among developing retinotectal synapses. STDP has also been broadly applied in many neocortical layers and brain regions. Hebbian learning can be applied in SNN to improve the separation property when LSM is used to deal with real-world speech data [11]. In our previous work [12], a novel type of LSM was obtained via STDP learning. The heterogeneity in the excitability degrees of neurons caused the SNN to evolve into a sparse and active neuron-dominant structure, wherein strong connections were mainly distributed to the outward links of a few highly active neurons. This LSM with an active neuron-dominant structure has better computational capability for real-time computations of complex input streams than the traditional LSM model with a random SNN. However, topological changes induced by Hebbian or STDP learning will impact network dynamics and encourage altered firing patterns [16, 17, 18, 19]. These altered firing patterns can further induce changes in connectivity through the plasticity mechanisms. Understanding and modeling the interactions between network structure and dynamics are essential for establishing and applying SNNs. Besides synaptic plasticity, network dynamics will inevitably involve non-synaptic plasticity, which includes modification of neuronal excitability in the axon, dendrites, and soma of an individual neuron [20, 21]. Recent experimental results have shown that the intrinsic excitability of individual biological neurons can be adjusted to match the synaptic input by the activity of their voltage-gated channels [22, 23, 24]. This adaption of neuronal intrinsic excitability, called intrinsic plasticity (IP), has been observed in cortical areas and is important for cortical functions of neural circuits [25]. IP has been hypothesized to keep the mean firing activity of neuronal population in a homeostatic level [26], which is essential for avoiding highly intensive and synchronous firing caused by STDP learning. In fact, IP learning can be taken as an adaptive rule: a single neuron can strengthen the excitability when its input is weak and weaken the excitability when the input is boosted. This adaptive adjustment of the neuronal input-output response online is crucial for efficient information processing. However, IP learning is 3

not involved in most of current studies on computational modeling of neural network. Composite studies on multi-neural plasticity in combination with existing network learning algorithms are important to maximize information capacity [27, 28, 29]. In this work, We propose a novel approach to develop the LSM reservoir with a biologically inspired self-organizing SNN based on both STDP and IP learning. The main features and contributions of this work can be summarized as follows: 1) In the proposed method, the connectivity among excitatory neurons is self-organized by STDP learning, whereas IP learning regulates the degree of neuronal excitability to moderate the average activity level. 2) Temporal coding and post-synaptic potentials (PSPs) are incorporated into our IP model using spiking neurons called the Izhikevich model. This IP learning model has been verified to be consistent with experimental results; 3) The computational capability of LSM with different reservoir structures has been examined and evaluated using average normalized mean squared errors (NMSE), and the results show that the LSM with the proposed SNN model can significantly reduce the average NMSE. The refined reservoir structure by STDP+IP learning exhibits larger entropy of activity patterns than other networks, which may contribute to its higher efficiency in information processing. The rest of the paper is organized as follows. Section II describes the model and the learning rules used. Section III presents the simulation results. Conclusions are made in Section IV. 2. Methods 2.1. Neuron Model Izhikevich model, which has been shown to be both biologically plausible and computationally efficient[30], is used for all neurons in the SNN. Specific mathematical expressions are as follows: v˙i = 0.04vi2 + 5vi + 140 − ui + I + Iisyn u˙i = a(bvi − ui ) + Dξi { vi → c if vi > 30 mV, then spike occurs and ui → ui + d

4

(1) (2)

where vi represents the membrane potential and ui is a membrane recovery variable. The variable ξi is the independent Gaussian noise with zero mean and intensity D that represents the noisy background. I stands for the external applied current, and Iisyn is the total synaptic current through neuron i as explained below. The model can exhibit the firing patterns of all known types of cortical neurons with the choice of parameters a, b, c, and d given in [30]. In this model, b describes the sensitivity of the recovery variable to the subthreshold fluctuations of membrane potential. The critical value for Andronov−Hopf bifurcation is b0 ≈ 0.2. For b < b0 , the neuron is in the rest state and is excitable; for b > b0 , the system has a stable periodic solution generating action potentials. Hence, parameter b, which governs the degree of neuron′ s excitability, is a critical parameter that can significantly influence the dynamics of the system. Neurons with larger b are prone to exhibit larger excitability and fire with a higher frequency than others. Here we take into account the heterogeneity of neurons in the network by giving different value of the parameter b in the neuron model. That is, the threshold for generating action potential is different for each neuron. In this paper, to achieve heterogeneity, the initial values of b are randomly distributed in [bmin , bmax ]. By applying the IP learning, the values of b for each neuron will be adjusted in real time. Others parameters a, c, and d for each neuron are set to be typical values (Table IA) for regular spiking pyramidal neurons (a typical excitatory neuron). Besides, the application of noise can enhance the degree of network heterogeneity to avoid sensitive over-synchronous activity and it is also important for the generation of activeneuron-dominant structure in this model. Thus noise is considered in our neuron model and Poisson spike trains are used to mimic the realistic input synaptic currents (see Section III for details). 2.2. Synapse Model Two kinds of synapse model are widely used in computational neuroscience, namely, current- and conductance-based synapses. In real neural systems, synaptic current induced by an input spike is dependent on the voltage of the post-synaptic neuron. Current-based synapse models, which neglect the voltage-dependent properties of synaptic currents, are widely used for analyzing robust synchronization in neural networks because of their simplicity. However, the voltage-dependent properties of synaptic currents play a key role in robust network synchronization in most cases [31], implying that current-based synapses are oversimplified for analyzing robust network 5

synchronization. Thus, the conductance-based synapse model, in which postsynaptic currents are dependent on the membrane potential of post-synaptic neurons, is considered in this study: Iisyn = −(vi − Esyn )

N ∑

gji sj

(3)

j=1

where N is the number of neurons; Synaptic currents Iisyn transmit PSPs to the postsynaptic neuron when a spike is generated in the presynaptic neuron. Esyn is the reversal potential. The value of Esyn depends on the type of presynaptic cell. As only excitatory neurons are considered, Esyn can be set as 0. gji is the synaptic weight from neuron j to neuron i, with i ̸= j. During the off-line training of SNN with STDP rule, it will be updated and restricted to the interval [0, gmax ]. The synaptic variable sj in Eq. (3) is defined as: s˙j = α(vj ) · (1 − sj ) − sj /τf α(vj ) = α0 /(1 + e−vj /vshp )

(4)

where τf is the time delay constant for the fall duration of the response. The synaptic recovery function α(vj ) is similar to the Heaviside function. When the presynaptic cell is in the silent state vj < 0, s˙j can be reduced to s˙j = −sj /τf ; otherwise, sj rises quickly and acts on the postsynaptic cells. Parameters used in Eq. (3) and Eq. (4) are listed in Table I(B). 2.3. STDP learning and IP learning The SNN model or the reservoir in this study has a self-organized topology and varying neuronal excitabilities, essentially differing from most LSMs [6, 32, 12]. The reservoir in our model is a network of heterogeneous neurons with different behaviors or degrees of excitability. Its topology and intrinsic excitability are formed by two unsupervised learning rules: STDP and IP. The introductions on these learning rules are as follows. 2.3.1. A. STDP learning In brain networks, it has been observed connection between two neurons is reinforced if the post-synaptic neuron fires shortly after the pre-synaptic neuron, while it is weakened when the pre-synaptic neuron fires after the post-synaptic neuron [33]. It provides a function for long-term potentiation or depression of synapses based on the time difference ∆t = tj − ti between 6

a single pair of pre- and postsynaptic spikes, in neuron i and j, respectively (Fig. 1c). The synaptic conductance is updated by

F (∆t) =

{

△gji = gji F (△t)

(5)

A+ exp(−∆t/τ + ) if ∆t > 0 −A− exp(∆t/τ − ) if ∆t ≤ 0

(6)

with ∆t = tj − ti , where tj and ti are the spike time of the post-synaptic neuron j and pre-synaptic neuorn i respectively. τ + and τ − determine the temporal window for synaptic modifications. The parameters A+ and A− determine the maximum amount of synaptic modification, which occur when ∆t is close to zero. Parameters for equations introduced above are listed in Table I. The chosen parameter values are typical settings based on previous literatures [30]. After STDP learning synaptic weights in the liquid network will always converge to stable values. Simulation results are robust to these parameter values in a relatively large varying range. 2.3.2. B. IP learning The intensity of an average synaptic input in the brain may change dramatically. Neurons maintain responsiveness to both small and large synaptic inputs by regulating intrinsic excitability to promote stable firing. This way, neuronal activity can keep from falling silent or saturating when the average synaptic input falls extremely low or rises significantly high [20]. The evidence for persistent changes in intrinsic excitability called IP has been proven by training animal behavior and artificial patterns of activation in brain slices and neuronal cultures [34]. Specifically, IP may allow neurons to transmit the maximum information following a given stationary level of metabolic cost [35]. In [36], Stemmler first explored the idea that IP may contribute to an approximate exponential distribution of the firing rate level, where exponential distribution has the highest entropy among the distributions of a non-negative random variable with a stationary mean. These exponential distributions have been observed in visual cortical neurons [35]. In [26], an IP learning rule has been developed for analog neurons (neuron models that do not incorporate PSPs). The neuronal firing rates follow an exponential distribution by adjusting two parameters. The ideas of this adaptation method have been extended to a commonly used non-linear and Gaussian output distribution [37, 38]. The 7

a Depression of intrinsic excitability

Potentiation of intrinsic excitability

- depressing event

+ potentiating event

neuron

neuron

b % Control

ISI1

modification in intrinsic excitability

1

...

2

ISIN

Tmax neuron becomes PRUH sensitive

1 0 neuron becomes

-1 OHVV sensitive

N 10

250

ISI [ms]

c temporal code 4

 

STDP

3



intrinsic modification

2



synaptic connection

1 20ms

synaptic strength + g

ZOSK

Intrinsic Plasticity

Figure 1: STDP learning and Intrinsic Plasticity (IP). (a) Examples of plasticity in intrinsic neuronal excitability. Neuronal excitability would be potentated if neurons fire in low frequency, while it is depressed when neurons fire in high frequency. In particular, it encourages every unit to participate in the network dynamics on a regular basis. (b) The inter-spike interval (ISI) dependent IP learning rule proposed in this paper. During the learning process, the most recent ISI is used to adjust the neuronal excitability: If ISIi is larger than the threshold Tmax , the neuronal excitability is strengthened to make the neuron more sensitive to input stimuli; otherwise, the neuronal excitability is weakened to make the neuron less sensitive. In this way, the excitability of each unit can be adjusted in a homeostatic manner. (c) Combined learning in the reservoir (for clarity, only four neurons are drawn). In the proposed method, the connectivity among excitatory neurons is self-organized by STDP learning (red arrow lines), whereas IP learning (blue dashed lines) regulates the degree of neuronal excitability to moderate the average activity level.

classic ESN has shown to benefit from reservoirs, which are pretrained on the given input signal with the IP rule. In [39], a simple generalization of the IP model has been proposed in [26]; however, Weibull was used rather than exponential distribution, which is more suitable for describing the different firing-rate distributions of different neurons in various brain areas among species. 8

Although the biological mechanisms have not been precisely defined, IP will constrain the firing rate to obey maximum entropy probability distribution, which can help networks to maintain stable firing rate distribution in the existence of stable input changes. According to the principle of maximum entropy, i.e., the firing rate level with exponential distribution, the Gaussian or Weibull distribution is suitable for different conditions. However, the aforementioned IP models are for highly simplified artificial neuron models that do not incorporate PSPs. However, the IP rule used in [36] has been based on the Hodgkin-Huxley model with a high computational complexity that is unsuitable for network computation. Thus, in this paper we propose a new IP model based on the Izhikevich neuron [30]. The parameter b in Eq. (1) which governs the degree of neuron excitability is updated by bi → bi + bmax · ϕi , where ϕi is defined as (Fig. 1b):  T −ISIi ) if ISIi < Tmin  −ηIP · exp( min Tmin ISIi −Tmax ϕi = (7) η · exp( Tmax ) if ISIi > Tmax  IP 0, others where ηIP is learning rate. The neuronal inter-spike interval (ISI) is ISIik = tk+1 − tki , where tki is the time of the kth firing of neuron i; Tmin and Tmax i are thresholds, they determine the expected ranges of ISI. During the learning process, the most recent ISI is examined every tc time and used to adjust the neuronal excitability: If ISIi is larger than the threshold Tmax , the neuronal excitability is strengthened to make the neuron more sensitive to input stimuli; if ISIi is less than the threshold Tmin , the neuronal excitability is weakened to make the neuron less sensitive to input stimuli. Parameters for equations introduced above are listed in Table I. To get an insight into how the firing rate evolves during the IP learning process, a case study is performed and the histogram of firing rate during IP learning for a randomly connected network is shown in Fig. 2, from which an 2 /2σ 2 ] √ normal distribution of firing rate (p(x) = exp[−(x−µ) ) is observed. This σ 2π result is consistent with the maximum-entropy theory that distribution is Gaussian if the desired variance is fixed. It indicates that our IP model is reasonable.

9

Table 1: Network Parameters for all simulations described in Section II and III

(A) Neuronal parameters, used in (3) and (4) dt (time step) 0.05ms (B) Synaptic parameters, used in (5) and (6) Esyn (exc.) 0 (C) Parameters used in off-line learning of LC I 6 (D) Learning parameters, used in (5) τ+ 20ms (E) Learning parameters, used in (6) tck 50ms

a 0.02

bmin 0.12

bmax 0.2

c −65mV

α0 3

τf 2ms

vshp 5mV

gmax 0.03

τ− 20ms

A+ 0.005

A− 0.00525

Tmin 90ms

Tmax 110ms

ηIP 0.012

D 0.1

3. Results 3.1. Effects of STDP and IP on network structure and dynamics To better understand the underlying mechanisms of the effects of STDP and IP learning on network computational performance, we investigate the combined effect of STDP and IP on network structure and dynamics here. For simplicity, in this section all synaptic weights (gij ) are normalized in the interval of [0, 1]. In [40, 16], the authors proposed a novel neural network with an active neuron-dominant structure, which is self-organized via the STDP learning rule. In this model, strong connections are mainly distributed to the outward links of a few highly active neurons. Moreover, a recent experimental study [41] found that a small population of highly active neurons may dominate the firing in neocortical networks, suggesting the existence of active neuron-dominant connectivity in the neocortex. Such synaptic distribution has been shown to be beneficial for enhancing information transmission of neural circuits [40, 16]. In the initial state, each neuron (in the same group of LC) is considered bidirectionally connected with the same weight values and subject to the 10

d 8

a

b

c

Figure 2: (a) Distribution histogram of firing rate for: network without IP, network during the IP learning and the final network after IP learning. (b) Distribution histogram of firing rate for the final network after IP learning and its its Gaussian fitting with parameters: µ = 9.998, σ = 0.730. (c) The probability plot of firing rates after IP learning, which is well fitted to the Gaussian/normal distribution.

same external current. Thus at the learning beginning, the excitability of each neuron is totally dependent on the parameter b in the model, where a large value of b near bifurcation point leads to high excitability. That is, the threshold for generating action potential is different for each neuron. This feature is important for the emergence of the active-neuron-dominant structure (see Fig. 3). This process mediates the internal dynamical properties of different neurons and renders the whole network more synchronous and therefore more sensitive to weak input. Fig. 3 shows the updated network structure by STDP alone or STDP+IP after 10 s of learning (2x105 time step). Fig. 3a indicates the active neuron-dominant structure (the active neurons with high excitability have strong outward synaptic connections, and the less active neurons with low excitability have strong inward synaptic connections) obtained from STDP learning. The IP strengthens the competition among different neurons and makes the connectivity structure more complex. Fig. 3b shows that most of the synapses in the STDP condition 11

are rewired to be either 0 or 1; however, the distribution is not bimodal but rather skewed toward smaller values in the case of STDP+IP because of the heterogeneity of intrinsic neuronal properties. The degree distribution for different networks are also examined. Based on the definitions of degree distribution for a weighted but undirected network in [42], we extend it to be applicable to our weighted and directed (out) (in) and in-degree ki are defined, respectively, networks. The out-degree ki as the sum of the weights of all out-directed links and all in-directed links attached to node i , ∑ ∑ (out) (in) ki = gij , ki = gji . (8) ∏ j∈ (i)

j∈

∏

(i)

Fig. 3c demonstrates that the STDP+IP network or STDP alone has a wider degree distribution than the STDP+IP shuffled network. For the learned networks, neurons with larger out-degrees have smaller in-degrees (Fig. 3d). Moreover, the in- and out-degrees exhibit an approximately linear (STDP case) or exponential (STDP+IP case) relationship. When shuffling the weight distribution in STDP+IP case, these principles disappear, which partly explains the higher errors of STDP+IP shuffled network than that of STDP+IP network (Fig. 10). Examples of activity patterns during a 300 ms period for different networks (neuron number is N = 100) are shown in Fig. 4. Without STDP, neurons with high degrees of excitability tend to develop seizure-like activity bursts that may consume high energy during the signal transmission. STDP update can increase the synchronization degree of network activity. When IP is applied, the network becomes unsynchronized again but shows more complex dynamics of network activity. Therefore, the efficiency of the network obtained via STDP+IP in signal processing is investigated by comparing the information entropy of its network activity. The information entropy measures the complexity of activity patterns in a neural network. The entropy H is defined as H=−

n ∑

pi log2 pi

(9)

i=1

where n is the number of unique binary patterns and pi is the probability that pattern i occurs. As a first approximation, pi is given by the number 12

post neuron j with bj inceasing

a STDP alone

STDP+IP

g /g

STDP+IP shuffled

ij

max

1

100

100

100

80

80

80

0.8

60

60

60

0.6

40

40

40

0.4

20

20

20

0.2

20 40 60 80 100 pre neuron i with bi inceasing

20 40 60 80 100 pre neuron i with bi inceasing

20 40 60 80 100 pre neuron i with bi inceasing

0

b 6000

6000

6000 STDP+IP

number

STDP alone

STDP+IP shuffled

4000

4000

4000

2000

2000

2000

0

0

0.5 gij/gmax

0

1

0

0.5 gij/gmax

0

1

0

0.5 gij/gmax

1

c STDP alone

50

50

40 30 20

40 30 20

10

10

0

0

0

50 degree

100

out degree in degree

percent [%]

out degree in degree

percent [%]

percent [%]

STDP+IP shuffled

STDP+I P

50

out degree in degree

40 30 20 10

0

20

50 degree

80

0

100

0

20

40 degree

100

d 100

50

0 0

50 in degree

100

40 out degree

data fit

out-degree

out degree

STDP+IP shuffled

STDP+I P

STDP alone 100

50

0

0

50 in degree

100

35 30 25 20 20

30 in degree

40

Figure 3: Network structures obtained from learning rules. (a)Image of the normalized synaptic matrix GN = (gij /gmax ). Both the x-axis and the y-axis represent the index of neurons which is sorted with increasing value of the parameter b (i.e. increasing excitability degree). (b) Histogram of the normalized synaptic weights. (c) The degree distribution histogram. (d) Scatter plots of neurons’ in-degrees and out-degrees. Note that the activeneuron-dominant results for STDP alone are consistant with those obtained in [40].

of spikes occurring in a given bin of the interval histogram divided by the sum of all the spikes in the bin[43]. For calculation convenience neuronal activities are measured in pattern units consisting of a certain number of neurons. In each time bin, if any neuron of the unit is firing then the event of this unit is active; otherwise it is inactive. (e.g. if the unit size is set to 10, then n = [N/10].) Fig. 5 shows that the entropy of the STDP+IP network 13

b

Random

c

STDP alone

STDP+IP 100

80

80

80

60 40

60 40 20

20 0 450

neurons

100

100

neurons

neurons

a

500

550

600

650

700

750

0 450

60 40 20

500

time (ms)

550

600

650

700

0 450

750

time (ms)

500

550

600

650

700

750

time (ms)

Figure 4: Spiking trains of different networks with current injection (I = 6) onto all of the neurons. The red dot describes the top 10 most active neurons with higher b than the others. (a) Random network with uniformly distributed synaptic weights; (b) Network with only STDP learning; (c) Network obtained by STDP and IP learning.

is larger than those of the other networks, making it robust for the unit size changes. This finding indicates that the highly complex dynamics of neuronal activity of the STDP+IP-refined neural network benefit for improving the efficiency of information processing.

Entropy H (bit)

10 STDP+IP STDP alone Random

8 6 4 2 0 5

10

15

20

unit size Figure 5: Entropy of network activity for networks with different topologies, where time bin is set to be 4ms. It is noted that the STDP+IP network exhibits much larger entropy than the other networks, which is robust for the changes of unit size.

14

3.2. Computational Capability of LSM with STDP+IP The LSM model in this study comprises input component (IC), LC, and RC (See Fig. 6), similar to the traditional LSM [6]. In the model, inputs in the first layer (IC) are temporal sequences of spikes rather than firing rate, which are generated randomly by Poisson process. The second layer (LC) is a recurrent SNN. To facilitate the test processing of network computation capability to signals with different characteristics, LC (the reservoir) is equally divided into nine groups, and each group has been given an independent input stream. Four input streams are selected (See details in Section 3) and 400 neurons to obtain four independent groups with 100 neurons in each group. Input streams are applied to LC with synaptic weight wil. For the third layer (RC), only one linear readout neuron has connections from all of the reservoir neurons in LC. The output synaptic weights are trained by linear regression or different teaching signals (See details in Section 3). In the following section, the effect of connection styles between groups in LC on the computational performance are analyzed. IC

LC (reservoir)

RC

wout

wil

Input Streams

Readout Neuron

Figure 6: Network structure. In this model, neurons in LC are equally divided into four groups where each group receives one input stream independently, i.e. each group received input only from one of the four input streams. For each group, inputs in IC are fully connected to neurons in LC, and wil is the synaptic weight. In LC, synaptic weights and neuron parameters in each group are generated by the updating with different plasticity learning rules. Input streams are applied to LC with synaptic weight Wil . For RC, there is only one readout neuron with connections from all of the reservoir neurons in LC. The output synaptic weights are trained by linear regression for different targets.

15

Input Streams 120 r4(t)

25 24

0

17

0

120 r3(t)

16

120 r2(t)

9

0

8

firing rate (Hz)

spike train idx

32

120 r1(t)

0 0

0.2

0.4

0.6

0.8

1

0

time (s)

Figure 7: Four independent input streams. Each consists of eight spike trains generated by Poisson process with varying rates ri (t), i = 1, ...4.

In this subsection, comparisons of computational capability between LSMs with reservoir topologies are obtained randomly or generated via the proposed STDP+IP-based self-organizing and leaning rules. Simulations are conducted using CSIM package 1, which is a neural network simulator written in Matlab [44]. The computing procedure consists of three steps: (1) Input generation of IC circuits Fig. 7 shows that IC contains four independent input streams, each having eight independent spike trains generated by Poisson process with randomly varying rates ri (t), i = 1, 2, 3, 4. The time-varying firing rates ri (t) are chosen as follows. The baseline firing rates for streams 1 and 2 are set to 5 Hz, with random distributed bursts of 120 Hz for 50 ms. The rates for Poisson processes that generated the spike trains for input Streams 3 and 4 are periodically updated by randomly selecting two firing rates, namely, 30 Hz or 90 Hz. The time span of the sequence is 1 s. The lines in Fig. 7 are the actual firing rates (i.e., spikes are counted within a 30 ms window) resulting from the spike trains. (2) Offline learning of LC circuits Each neuron receives background input I, weak noisy input, and synaptic input Isyn as learning environment. Before training specific tasks, synaptic weights between neurons in each LC group are first updated offline using STDP learning, which has been discussed in the previous subsection (Fig. 1c). Meanwhile, in the Izhikevich model, parameter b governs the degree of neu16

ronal excitability; thus, the modifications of b are considered as a representative scheme describing IP mechanisms. This way, the sensitivity of neurons with low-firing activity to input stimuli will be improved, whereas neurons that fire in high frequency will have lower neuronal excitability. (3) Training readout neurons In order to investigate the computational capability of LSM, usually learning tasks could be any linear or nonlinear combinations of several independent inputs which are generated from eight Poisson spike trains with randomly varying rates r(t). In these tasks, target signals are selected as mathematical functions of the firing rate of input streams. The teaching NMSE signals are r1 , r2 , r3 , r4 , r2 + r4 , and r4 2 + r1 · r3 . The average normalized mean squared errors (NMSE) between target and observed values are calculated for a set of tasks that are similar to the original LSM in the study of Maass [44] to evaluate the general computational capability. The definition of NMSE is: 1 n

N M SE =

n ∑

(ˆ yi − yi )2

i=1 n ∑

1 n−1

(yi − y¯i

(10)

)2

i=1

where yi is the target signal, yˆi is the observed value, and y¯i is the average target value of n sets of testing data. Connections and parameters in LC remain unchanged and only readouts are trained through linear regression during the training of readout weights. The signals from LC to readout neurons consist of spikes rather than firing rates. The signal to each readout neuron should be transferred to synaptic currents from the presynaptic spike trains of all neurons in LC. Transformation occurs by applying a filter using an exponentially decaying kernel [44]. A single readout neuron is connected to all of the Izh neurons in the liquid network. Each neuron in the liquid network provides its final state value to the readout neuron, scaled by its synaptic weight. The final state value xi of the ith liquid neuron is calculated based on the latest spike time of the ith neuron: { x(i) + eδ/τ 0 ≤ δ/τ ≤ 30 x(i) = (11) x(i) + 0 δ/τ ≥ 30 where τ is a time constant set to be 30 ms. δ = tsample − tspike , tspike is the spike time and tsample is the current sample time. Then, the output of readout 17

neuron is y = xT ∗Wout , where Wout is the weight matrix trained by the linear regression algorithm (i.e., using the function of regress() in Matlab). Fig. 8 shows the output of the readout neuron and target for task training in LSM with STDP+IP. The results show that readouts can be well trained for simple (Fig. 8a) and complex nonlinear (Fig. 8b) combinations. The sensitivity of the model performance is also tested with respect to variations in the parameter wil (synaptic weights from IC to LC) to characterize and quantify the dynamics of networks systematically. Average NMSE for trained time series with different wil settings is shown in Fig. 9. LSM with STDP+IP performs better than those with random reservoirs for the varying ranges of synaptic weights. Moreover, LSM, whose reservoir is generated by STDP+IP learning, performs better than those generated by STDP alone. STDP learning can reduce a number of unnecessary synaptic connections. After STDP, the internal dynamics of neurons with different degrees of excitability are extracted into an active neuron-dominant topology (See details in Section 4), which contributes to the spike synchronization of the neural network. However, high-degree global intrinsic synchronization impairs computational performance. Adding IP can increase the information entropy, which is beneficial for reducing the likelihood of pathological synchronization, thereby contributing to the significant enhancement of transmitting time-varying signals (Fig. 9). a

r2(t)+r4(t)

300

output

250

target

200 150 100 50 0

0

0.1

0.2

0.3

0.4

0.5 time (s)

0.6

0.7

0.8

0.9

1

0.9

1

r4(t)2+r1(t)*r3(t)

b

output target

15000 10000 5000 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

time (s)

Figure 8: Firing activity of the readouts (dashed) and the targets (solid) for 2 sample tasks training in LSM with STDP+IP. (a) linear combination r2 + r4 ; (b) nonlinear combination r42 + r1 · r3 . Here we set wil = 0.2.

18

0.1

Average NMSE

0.08

0.06

0.04 random STDP STDP+IP

0.02

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

wil

Figure 9: Comparisons of computational capability of LSMs with different reservoirs (STDP+IP: LC constructed by STDP and IP learning; STDP: LC constructed by STDP alone; Random: LC with uniformly distributed synaptic weights in the range of [0,gmax]). wil is the synaptic weight from IC to LC and varies from 0.1 to 0.8 with a step of 0.1.

Five different connection types among groups have been examined to investigate the influence of dynamic communication among groups in LC on the computational performance of LSM (Fig. 7). These cases are 1) no connections between groups; 2) 10% randomly selected neurons in each group are interconnected; 3) 10% most excitable neurons (with the top 10% largest values of the parameter b) in each group are interconnected; 4) 20% randomly selected neurons in each group are interconnected; and 5) 20% most excitable neurons in each group are interconnected. Fig. 10 shows that for all these connection types, LSM with STDP+IP has the minimum errors. Moreover, NMSE is enlarged as connected percentage among groups in LC increases. Compared with the other cases, LSM with STDP+IP is more robust for changing the connections among groups. These results indicate that the good performance of STDP+IP LSM with independent groups is not only caused by the distribution of G matrix. Overall, connections among groups can weaken the computational capability because these connections interrupt the dynamical states of each subnetwork. This interruption causes the confusion in the information stored in each separate group, thereby impairing LSM computational capacity. 4. Conclusion In this study, a new SNN model based on STDP+IP learning rules is developed for LSM. The internal dynamics of different neurons are encoded 19

0.14 0.12

random STDP alone STDP+IP STDP+IP-shuffled

Average NMSE

0.1 0.08 0.06 0.04 0.02 0

no rand10% big10% big20% rand20% different connection styles between groups in LC

Figure 10: Computational capability of LSM with different connection styles between groups in LC: no connections between groups; connections happen in 10% or 20% of the randomly selected neurons in each group; connections happen in 10% or 20% of the most excitable neurons in each group. Here the shuffled STDP+IP condition means that the weight values are the same with the STDP+IP condition but their orders are randomly shuffled. (wil = 0.05 after normalization.)

in the topology of the emergent network after learning. LSM has smaller NMSEs for real-time computations than the traditional LSM model with a random SNN because of the existence of the self-organization structure in the SNN. Moreover, the LSM with a SNN generated by STDP+IP learning performs better than the that generated by STDP alone. Furthermore, the activity entropy of the SNN is examined, and results show that the network obtained using STDP+IP learning has more entropy than the other networks, indicating its high capability of signal propagation. As the STDP training of liquid network is independent from specific learning tasks, once the liquid updating is finished, the liquid connections always keep unchanged and only readouts weights are trained through linear regression, which can be achieved in one time step. Thus, our proposed SNN model optimized by STDP+IP is efficient in performing real-time computational tasks. The noticeable improvement of computational performance of the proposed method is due to the better reflected competition among different neurons in the developed SNN model, as well as the effectively encoded relevant dynamic information with its learning and self-organizing mechanism.

20

Another issue for investigation is the contribution of the biological structure of the SNN to the optimization of LSM computation when inhibitory neurons in the neural networks are considered, and commonly observed learning rules are included. At present, the number of neurons in our study is still small, however, the existing parallel CPU implementation simulator or GPU-based simulator may support for future design and analysis of very large spiking neural models. Further studies on the computational model of SNN with larger network size and much richer spatio-temporal structure are necessary and would be our future study. Acknowledgment This work is supported by the National Natural Science Foundation of China (Nos. 61473051 and 61304165), Natural Science Foundation of Chongqing (No. cstc2016jcyjA0015) and Fundamental Research Funds for the Central Universities (No. 106112017CDJXY170004). [1] M. Boerlin, S. Denve, Spike-based population coding and working memory, Plos Computational Biology 7 (2) (2011) e1001080. [2] A. Renart, P. Song, X. J. Wang, Robust spatial working memory through homeostatic synaptic scaling in heterogeneous cortical networks 38 (3) (2003) 473–485. [3] X. J. Wang, Probabilistic decision making by slow reverberation in cortical circuits., Neuron 36 (5) (2002) 955. [4] M. Boerlin, C. K. Machens, S. Denve, Predictive coding of dynamical variables in balanced spiking networks, Plos Computational Biology 9 (11) (2013) e1003258. [5] L. F. Abbott, B. Depasquale, R. M. Memmesheimer, Building functional networks of spiking model neurons., Nature Neuroscience 19 (3) (2016) 350. [6] W. Maass, T. Natschl¨ager, H. Markram, Real-time computing without stable states: A new framework for neural computation based on perturbations, Neural computation 14 (11) (2002) 2531–2560. 21

[7] H. Jaeger, The echo state approach to analysing and training recurrent neural networks-with an erratum note, Bonn, Germany: German National Research Center for Information Technology GMD Technical Report 148 (2001) 34. [8] J. J. Steil, Backpropagation-decorrelation: online recurrent learning with O(N) complexity, in: Neural Networks, 2004. Proceedings. 2004 IEEE International Joint Conference on, Vol. 2, IEEE, 2004, pp. 843– 848. [9] D. Verstraeten, B. Schrauwen, D. Stroobandt, Reservoir computing with stochastic bitstream neurons, in: In Proceedings of the 16th Annual ProRISC Workshop, 2005, pp. 454–459. [10] D. Verstraeten, B. Schrauwen, M. dHaene, D. Stroobandt, An experimental unification of reservoir computing methods, Neural networks 20 (3) (2007) 391–403. [11] D. Norton, D. Ventura, Preparing more effective liquid state machines using hebbian learning., in: IJCNN, 2006, pp. 4243–4248. [12] F. Xue, Z. Hou, X. Li, Computational capability of liquid state machines with spike-timing-dependent plasticity, Neurocomputing 122 (2013) 324–329. [13] F. Triefenbach, K. Demuynck, J.-P. Martens, Large vocabulary continuous speech recognition with reservoir-based acoustic models, Signal Processing Letters, IEEE 21 (3) (2014) 311–315. [14] E. A. Antonelo, B. Schrauwen, On learning navigation behaviors for small mobile robots with reservoir computing architectures, IEEE transactions on neural networks and learning systems 26 (4) (2015) 763–780. [15] H. Markram, J. L¨ ubke, M. Frotscher, B. Sakmann, Regulation of synaptic efficacy by coincidence of postsynaptic aps and epsps, Science 275 (5297) (1997) 213–215. [16] X. Li, M. Small, Enhancement of signal sensitivity in a heterogeneous neural network refined from synaptic plasticity, New Journal of Physics 12 (2010) 083045. 22

[17] S. Klampfl, W. Maass, Emergence of dynamic memory traces in cortical microcircuit models through stdp, The Journal of Neuroscience 33 (28) (2013) 11515–11529. [18] J. Zhang, K. M. Kendrick, G. Lu, J. Feng, The fault lies on the other side: Altered brain functional connectivity in psychiatric disorders is mainly caused by counterpart regions in the opposite hemisphere, Cerebral Cortex 25 (10) (2014) 3475–3486. [19] J. Zhang, W. Cheng, Z. Liu, K. Zhang, X. Lei, Y. Yao, B. Becker, Y. Liu, K. M. Kendrick, G. Lu, Neural, electrophysiological and anatomical basis of brain-network variability and its characteristic changes in mental disorders, Brain 139 (8) (2016) 2307–2321. [20] N. S. Desai, L. C. Rutherford, G. G. Turrigiano, Plasticity in the intrinsic excitability of cortical pyramidal neurons, Nature neuroscience 2 (6) (1999) 515–520. [21] A. Destexhe, E. Marder, Plasticity in single neuron and circuit computations, Nature 431 (7010) (2004) 789–795. [22] G. Daoudal, D. Debanne, Long-term plasticity of intrinsic excitability: learning rules and mechanisms, Learning & Memory 10 (6) (2003) 456– 465. [23] R. H. Cudmore, G. G. Turrigiano, Long-term potentiation of intrinsic excitability in lv visual cortical neurons, Journal of neurophysiology 92 (1) (2004) 341–348. [24] J. Naud´e, J. T. Paz, H. Berry, B. Delord, A theory of rate coding control by intrinsic plasticity effects, PLoS computational biology 8 (1) (2012) e1002349. [25] E. Marder, L. Abbott, G. G. Turrigiano, Z. Liu, J. Golowasch, Memory from the dynamics of intrinsic membrane currents, Proceedings of the National Academy of Sciences 93 (24) (1996) 13481–13486. [26] J. Triesch, Synergies between intrinsic and synaptic plasticity in individual model neurons, Advances in neural information processing systems 17 (2004) 1417–1424. 23

[27] A. Lazar, G. Pipa, J. Triesch, Fading memory and time series prediction in recurrent networks with different forms of plasticity, Neural Networks 20 (3) (2007) 312–322. [28] A. Lazar, G. Pipa, J. Triesch, Sorn: a self-organizing recurrent neural network, Frontiers in computational neuroscience 3 (2009) 23. [29] P. Zheng, C. Dimitrakakis, J. Triesch, Network self-organization explains the statistics and dynamics of synaptic connection strengths in cortex, PLoS computational biology 9 (1) (2013) e1002848. [30] E. Izhikevich, Simple model of spiking neurons, Neural Networks, IEEE Transactions on 14 (6) (2003) 1569–1572. [31] Z. Wang, W. Wong, Key role of voltage-dependent properties of synaptic currents in robust network synchronization, Neural Networks 43 (2013) 55–62. [32] A. L. H. Burgsteiner, M. Krll, G. Steinbauer, Movement prediction from real-world images using a liquid state machine, Applied Intelligence 26 (2) (2007) 99–109. [33] S. Song, K. Miller, L. Abbott, et al., Competitive hebbian learning through spike-timing-dependent synaptic plasticity, nature neuroscience 3 (2000) 919–926. [34] W. Zhang, D. J. Linden, The other side of the engram: experience-driven changes in neuronal intrinsic excitability, Nature Reviews Neuroscience 4 (11) (2003) 885–900. [35] R. Baddeley, L. F. Abbott, M. C. Booth, F. Sengpiel, T. Freeman, E. A. Wakeman, E. T. Rolls, Responses of neurons in primary and inferior temporal visual cortices to natural scenes, Proceedings of the Royal Society of London. Series B: Biological Sciences 264 (1389) (1997) 1775– 1783. [36] M. Stemmler, C. Koch, How voltage-dependent conductances can adapt to maximize the information encoded by neuronal firing rate, Nature neuroscience 2 (6) (1999) 521–527.

24

[37] J. J. Steil, et al., Online reservoir adaptation by intrinsic plasticity for backpropagation-decorrelation and echo state learning, Neural Networks 20 (3) (2007) 353–364. [38] B. Schrauwen, M. Wardermann, D. Verstraeten, J. J. Steil, D. Stroobandt, Improving reservoirs using intrinsic plasticity, Neurocomputing 71 (7) (2008) 1159–1171. [39] C. Li, A model of neuronal intrinsic plasticity, Autonomous Mental Development, IEEE Transactions on 3 (4) (2011) 277–284. [40] X. Li, J. Zhang, M. Small, et al., Self-organization of a neural network with heterogeneous neurons enhances coherence and stochastic resonance, Chaos 19 (1) (2009) 3126. [41] L. Yassin, B. L. Benedetti, J.-S. Jouhanneau, J. A. Wen, J. F. Poulet, A. L. Barth, An embedded subnetwork of highly active neurons in the neocortex, Neuron 68 (6) (2010) 1043–1050. [42] I. Antoniou, E. Tsompa, Statistical analysis of weighted networks, Discrete Dynamics in Nature and Society 2008. [43] P. Matthews, Relationship of firing intervals of human motor units to the trajectory of post-spike after-hyperpolarization and synaptic noise., The Journal of physiology 492 (Pt 2) (1996) 597–628. [44] W. Maass, P. Joshi, E. Sontag, Computational aspects of feedback in neural circuits, PLOS Computational Biology 3 (1) (2007) e165.

25