SpiFoG: an efficient supervised learning algorithm for the network of spiking neurons

There has been a lot of research on supervised learning in spiking neural network (SNN) for a couple of decades to improve computational efficiency. However, evolutionary algorithm based supervised learning for SNN has not been investigated thoroughly which is still in embryo stage. This paper introduce an efficient algorithm (SpiFoG) to train multilayer feed forward SNN in supervised manner that uses elitist floating point genetic algorithm with hybrid crossover. The evidence from neuroscience claims that the brain uses spike times with random synaptic delays for information processing. Therefore, leaky-integrate-and-fire spiking neuron is used in this research introducing random synaptic delays. The SpiFoG allows both excitatory and inhibitory neurons by allowing a mixture of positive and negative synaptic weights. In addition, random synaptic delays are also trained with synaptic weights in an efficient manner. Moreover, computational efficiency of SpiFoG was increased by reducing the total simulation time and increasing the time step since increasing time step within the total simulation time takes less iteration. The SpiFoG is benchmarked on Iris and WBC dataset drawn from the UCI machine learning repository and found better performance than state-of-the-art techniques.

www.nature.com/scientificreports/ the network in computationally efficient manner as well as to add biological plausibility for a couple of decades. Bohte et al. 8 introduced a supervised learning algorithm for multilayer feed forward SNN called SpikeProp akin to the back-propagation algorithm. In SpikeProp algorithm, time-to-first-spike input encoding scheme was used to overcome the problem of discontinuity i.e., input, hidden and output neurons can fire only a single spike and the spikes after the first spike are discarded. Later, it has been found that although Bohte et al. 8 used a small value of learning rate, a large value of learning rate can also lead to a successful training 15 . To speed up the convergence rate, performance and to allow multiple spikes, SpikeProp is further investigated in [15][16][17][18][19][20] . The demerits of gradient rule based learning is the surge or sudden jumps, and discontinuities in the error function as discussed by Shrestha et al. 16 . Sporea et al. 21 had used spike time dependent plasticity (STDP) and anti-STDP to device more similar biological supervised learning for multilayer SNN. This method allows neuron in all the layers to fire multiple times. However, authors did not consider hidden neurons spike time while training, which is a barrier to the biological plausibility. In 22 , a survey on supervised learning for SNN is carried out where various supervised learning methods was compared. Quiang et al. 23 proposed another supervised learning algorithm that uses temporal encoding. Ponulak et al. 24 proposed remote supervised method (ReSuMe) and argued that the method possesses interesting properties of the supervised Hebbian approach while avoiding its drawback. The main advantage of ReSuMe method is that the method is independent of the neuron models used and can be used in a variety of neuron models. In extended delay based ReSuMe 25 , and in DL-ReSuMe 12 , a delay learning based remote supervision method is investigated for the network of spiking neurons. Note that in 12,25 , constant synaptic delays was used. In 26 , multiple neurons are trained instead of a single neuron. Aboozar et al. 27 proposed BPSL algorithm that can fire multiple desired spikes from a single neuron. John et al. 28 , proposed another STDP based supervised learning algorithm SWAT. In 29 , a supervised training algorithm is proposed where a neuron learns spike timing-based decisions called tempotron which allows neurons to learn whether to fire or not in response to a particular input stimuli. However, in tempotron neuron did not learn precise desired output spike train, rather decides the firing ability of the neuron. Ammar Mohemmed et al. 30, 31 developed a supervised learning algorithm using extended delta learning rule where each spike sequence is convolved with a suitable kernel function. J. Wang et al. 32 proposed OSNN algorithm having adaptive structure based supervised learning ability. There exist membrane potential driven learning algorithms proposed for SNN such as MPD-AL 33 , EMPD 34 , and MemPo-Learn 35 . In the development of MPD-AL 33 algorithm, two approaches were used namely firing less spikes than desired and firing more spikes than desired. In 34 , two process namely at desired output times and at undesired output times were used. The error function proposed in 35 is formed by taking difference between membrane potential of output neuron and the firing threshold of its own. The ability of evolutionary strategies (ES) to work directly with real numbers without complex encoding technique drives researchers towards the investigation of ES to train SNN. In 36,37 , SNN was trained using evolutionary techniques. Evolutionary methods also had been used to optimise SpikeProp in 38 where particle swarm optimization (PSO) technique was used. In 37 , a self regulating evolving SNN named SRESN was proposed. Differential evolution is also used to develop a supervised learning algorithm for SNN in 39 called DEPT-ESNN. The limitation of ES based learning algorithm is that these algorithm demands the neuron in SNN to fire exactly once. Therefore, ES technique cannot be used efficiently where there is multiple spiking in SNN. In addition, ES based learning technique is very time consuming. There is an existing gap between neuroscience and machine learning to date. Although SNN is well capable to bridge this gap, but remained unexplored for decades as it was considered too complex to implement practically and difficulty in analysing the performance. Apart from that, adding biological plausibility is a much bigger challenge in terms of computational cost. In addition, due to the lack of efficient supervised learning algorithm and difficulty in software implementation, SNN became less popular although it is energy efficient, provides good computational power and biologically plausible. Moreover, SNN is very challenging due to its discrete spiking nature. The biological plausibility of unsupervised learning such as STDP, long term potential (LTP) and long term depression (LTD) in SNN is investigated and used by many researchers over the last couple of decades 5,40 . However, unsupervised learning paradigms are not suitable for classification problems and the classification of spatio-temporal pattern becomes very essential in the current scenario. Also from literature review it is found that the development of an efficient supervised learning algorithm for SNN is still in its infancy stage. Therefore, this paper presents an efficient algorithm to train multilayer Spiking neural network which uses elitist Floating pointGenetic algorithm (SpiFoG).
The main contributions of this paper is as follows: (1) An efficient supervised learning algorithm (SpiFoG) for multilayer feed-forward SNN is proposed.
(2) Random delays are introduced in the synapse model of LIF neuron which are also trained with synaptic weights. (3) Both positive and negative values are allowed to add more biological plausibility to the system. (4) The total simulation time is reduced, and the value of time step is increased to improve computational efficiency. (5) Hybrid crossover method is used for the faster convergence providing well exploration of the search space.

Results
To evaluate the performance, SpiFoG is benchmarked against the Iris (multiclass) and WBC (binary) dataset drawn from the UCI machine learning repository. For every dataset, the experiment is repeated 10 times randomly.
Dataset description. The Iris dataset 41  Performance measures. The performance of the SpiFoG is evaluated in terms of Accuracy (Training set and Test set), and computational efficiency (CE) defined in the Eq. (1). Table 1 shows the comparison in terms of mean accuracy along with standard deviation, and in terms of CE which is a better way to compare among the algorithms since running time in different machine with different configuration is different. The more the value of CE, better is an algorithm.
where G is the epoch at which the algorithm converges and it is independent of total simulation time T (it is the maximum time at which a neuron must issue atleast a single spike and to do that efficient setting of parameters are very necessary). There is no correlation between G and T, former is related to tuning parameters at training phase and later is related to issuing spikes after which the tuning is possible.

Discussion
The selection of hyper-parameters such as dt, T, τ m , τ s , R m in an efficient manner is very essential and necessary in case of SNN since each one is somehow related to each other whether it is the production of spikes or in training phase. A learning algorithm can only tune parameters to reduce overall network error after the production of spike times otherwise there is no way to compute the value of error. In our case, the selection of few parameter effects the performance of SpiFoG both in terms of accuracy (training set and test set) and CE. In this experiment, the total simulation time denoted by T is set to 20 ms which is less than state-of-the-art techniques and it is directly responsible for providing better CE. The problem with in the reduction of the value of T is that we have to ensure that within that time neurons fire atleast a spike and it is directly dependent on hyper-parameter setting. On the other hand, due to the efficient synaptic weight distribution, the neurons needs less iteration to fire a spike. Time step dt is taken as 1 ms i.e., at a step size of 1 ms each neuron updates the membrane potential that provides more speed and thereby CE increases. But, if the parameter settings are not good such as initialisation of random weights then whatever may be the value of time step there is a very less chance for a neuron to fire within the value of T. Therefore, before training, the setting of the range for synaptic weights is needed to be done providing a good trade off between the threshold value and the average width of the distribution. A very low value of synaptic weight results in no firing of the neurons and a very high value of synaptic weight results in early firing of the neurons. The computation of error is not possible if the neurons do not fire at all and early firing forces the algorithm to convergence prematurely. Therefore, to make the algorithm effective we have chosen the values of synaptic weights from uniform distribution within the range U(-0.25, 1) allowing both positive and negative values and a mixture of inhibitory and excitatory neurons to make the system more similar to the biological neuron. In this experiment, random synaptic delays are introduced with LIF model having values from random integer distribution within the range (1, 6) which is also trained along with synaptic weights. The www.nature.com/scientificreports/ selection of minimum and maximum values of delays effect the speed of the algorithm. A very high value of delay unnecessarily delays the firing of a post-synaptic neuron and thereby provides a poor value of CE. There are other parameters responsible for the efficiency of SpiFoG such as double decaying kernel function, hybrid crossover method (discussed in "Methods" section) which provides better convergence of the algorithm. The kernel function having double decay was used for the experiment and the value of time constants τ m and τ s effects the change in membrane potential. We have chosen the value of τ m as 10 ms and the value of τ s as 5 ms. A small value of τ m and τ s results in a narrow shaped PSP and in that case, a high weight value is required to reach the threshold value as well as faces a problem of non overlapping of PSPs from pre-synaptic neurons. On the other hand, higher value of τ m and τ s results in a wide PSP and in that case, internal state of a neuron increases very fast which takes very long time to come to the resting potential as well as the overlapping of PSPs occur. From the experimental point of view, we have seen a value of τ m equal to or slightly greater than T and a value of τ s as half of τ m adds more stability to the shape of PSP. The utilisation of hybrid crossover leads the algorithm towards faster convergence although the search space is well explored. Figure 1a (top), shows the mean accuracy (training set) for SpiFoG algorithm in case of Iris dataset having the value of dt as 1 ms where, it seems that the algorithm converges at 299th generation out of 1000 generations. Note that, up to 50 generations, accuracy (training set) increases rapidly and after that there is a slow change. This phenomenon demonstrates both exploitation and exploration of search space. Figure 1a (bottom), shows the search space exploration (i.e., big jumps) in case of Iris dataset and we refer this as unique or distinct change in accuracy (training set). When no distinct value occurs while training at a generation an algorithm is said to be stagnate at that value which justifies exploitation. All 30 generations shown in Fig. 1a (bottom) have a high variation in accuracy (training set) that depicts the fast convergence behaviour of the algorithm but still provides a better exploration with a population size 100. Figure 1b (top), shows the mean accuracy (training set) in case of WBC dataset having the value of dt as 1 ms where, it seems that the algorithm converges at 896th generation out of 1,000 generations. It is observed from Fig. 1b (top) that after 200 generation the change in accuracy (training set) is minute still it is allowed to evolve to settle down at a static value. Figure 1b (bottom) shows the search space exploration (i.e., big jumps) in case of WBC dataset where all 53 generations have a high variation in accuracy (training set). Figure 2c (top), shows the behaviour of fitness values without elitism for 80 chromosomes out of 100 (i.e., 20 best chromosomes are kept for the next generation and it is clearly visible in Fig. 2c (top)) in case of Iris dataset. The fitness values of non elitist chromosomes keeps increasing and decreasing at each generation. But, after elitism best chromosomes are combined with the chromosomes of the next generation results in a faster convergence, because best solutions never lost. Figure 2c (bottom) shows, unique or distinct values of fitness (exploration) after elitism in case of Iris dataset where the value of fitness only raises upwards which always guarantees the selection of parents having best fitness values for the next generation. Similarly, Fig. 2d (top) and (bottom) shows the behaviour of fitness values without elitist selection for 80 chromosomes in case of WBC dataset and distinct values of fitness (exploration) after elitism in case of WBC dataset. Figure 3a,b shows the confusion matrix for the Iris dataset and the WBC dataset where classified and misclassified samples are clearly visible for each dataset respectively. In addition, Fig. 3c,d shows the trade-off between true positive rate (sensitivity) and false positive rate (1-specificity) for every possible cut-off or threshold values in case of Iris and WBC dataset respectively. Note that, the ROC curve for Iris dataset shown in Fig. 3c, represents area of the individual classes where Class 1 shows the highest value in terms of area and both micro and macro average area is also shown as Iris is a multiclass classification problem. But, in case of WBC dataset (Fig. 3d), a single area under ROC curve is shown since it is a binary classification problem where the area is 0.95. The co-ordinates shown in Fig. 3c,d represents the vales of sensitivity vs. (1-specificity) for different threshold values. We have carried out the experiment using two values of dt one is 0.1 ms and other is 1 ms to compare with state-of-the-art algorithms specially in terms of CE. Table 1 shows the accuracy (training set and test set) for both the datasets where SpiFoG outperforms other algorithms in case of Iris dataset both in terms of accuracy

Methods
The primary objective of the proposed algorithm SpiFoG is to train a multilayer feed forward SNN efficiently. The encoding of real valued input features is carried out using the time-to-first-spike encoding scheme 8 . The LIF spiking neuron with an additional term random synaptic delay, and a double decay kernel function is given in Eqs. (4), (5), and (6). Both weights and delays are optimised using SpiFoG algorithm.

Encoding of inputs and outputs.
Using time-to-first-spike encoding scheme discussed in 8 , each real valued feature is converted to a number of overlapping Gaussian functions called encoding neurons having centred at µ as given in Eq.
(2) and a spread σ as given in Eq. where β is the adjustment factor. The value of β should be in the range [1,2], but 1.5 is found to be experimentally best suited.  Fig. 4a (bottom). The output spike times are restricted to 13 ms so that the total simulation time can be minimised. Three classes in case of Iris dataset are represented by 11 ms, 12 ms, and 13 ms and two classes in case of WBC dataset are represented by 11 ms, and 12 ms (trial and error basis).
Proposed neuron model. The dynamics of spiking neurons is described by simple LIF neuron model 13,14 .
In this paper, the LIF neuron model having random synaptic delays is introduced which uses a double decaying synaptic model. The dynamics of the sub-threshold regime of a post-synaptic LIF neuron j is characterised in the Eq. (4).
where dν j (t) is the small change in membrane potential of a post-synaptic neuron j at time t, dt is the unit time step, R m is the membrane resistance, τ m is the membrane time constant, ν j (t − 1) is the membrane potential of post-synaptic neuron j at time t − 1 (when t = 1 ms, it becomes the resting potential of that post-synaptic neuron j), and ξ(t) is the input stimuli at the post-synaptic neuron j from pre-synaptic neurons i given in Eq. (5). The value of R m is taken as 100 k .
where M is the total number of pre-synaptic neurons and N is the total number of post-synaptic neurons, w ij ∈ R is the synaptic weights which describe the connection strength between two neurons i and j, t i denotes the spike time of pre-synaptic neuron i, and κ ij ∈ I represents random synaptic delays between the synapses of neurons i and j. The definition of synaptic current kernel function for an input x is expressed as �(x) which is given in Eq. (6). Finally, the membrane potential ν j at time t of the post-synaptic neuron j is updated using Eq. (8).
When the value of ν j reaches the threshold potential ϑ , the post-synaptic neuron j fires a spike and the spike time t j for neuron j is recoded. After firing a spike membrane potential resets to a reset potential in case of multiple spike firing scheme. In this paper, the first time when ν j reaches ϑ is considered and each neuron is allowed to fire only once. The spike times of neuron j is described by S j given in Eq. (9).
In Eq. (5), the term κ ij in the synapse model of the LIF neuron is introduced to make the network more biologically plausible and trained along with w ij to tune these parameters for the production of a better set of output spike times those are more closer to the target spike times.
Network topology. Figure 5a shows the network topology of the proposed multilayer feed forward SNN that consists of 17 input neurons, 12 hidden neurons, and 1 output neuron. The number of hidden neurons are selected on the basis of trail-and-error and only one hidden layer is used. Figure 5b,c show the synapse model where synaptic delay between pre-synaptic neuron j and post-synaptic neuron k denoted as κ jk is associated with the spike times t j , and the change in membrane potential upon receiving the input stimuli which uses a double decay kernel function, respectively. The amplitude of �(x) increases with the increase in the value of synaptic weights. The more the value of synaptic weights more rapidly the amplitude increases. Therefore, the value of weights should be sufficient enough to reach the threshold value, otherwise the post-synaptic neurons end up www.nature.com/scientificreports/ with no spike times due to no firing. The negative weight values decreases the amplitude of kernel rather than increasing.
Proposed learning method. The purpose of the proposed algorithm SpiFoG is to train the SNN architecture shown in Fig. 5a by optimising the initially generated random parameters such as synaptic weights w and synaptic delays κ . The optimised weights and delays forces the neurons to fire in a more sophisticated manner and then the firing spike times from hidden neurons acts as pre-synaptic spike times for the output neuron. The objective is to minimise the difference between actual spike times t k and desired spike times t d . After the generation of actual spike times for all instances at gth generation t g k , where there are few spike times equal to the desired spike times t d , the value of w and κ are adjusted to produce more fine tuned spike times at (g + 1) th generation t g+1 k , using SpiFoG. The training continues until t k very closely converges to t d . Firstly, upon receiving the inputs, the PSP of the hidden neurons changes. When the value of PSP of a hidden neuron j crosses the threshold ϑ at time t then the neuron j fires a spike at time t. The value of threshold ϑ is taken as 1 mV and the value of resting membrane potential ν r is taken as 0 mV at time t = 0 ms. It is essential to optimise t k to match t d i.e., t k = f (w jk , κ jk ) = g(w ij , κ ij ) = h(t j ) . The values of w and κ are tuned by forming the objective function f i.e., the fitness value of chromosomes and it is inversely proportional to the sum squared error (SSE) denoted by e and later it is converted to mean squared error (MSE) by taking mean in case of all training instances. In Eq. (10), the evaluation of fitness function is shown as the objective function which is the accuracy of the training set.
The definition of e is given in Eq. (11) where, t k i is the output spike times for the only output neuron for all training samples and t d i is the desired spike times for all training samples.
The main goal of SpiFoG is to maximise f because maximising f will minimise e. It uses elitist floating point genetic algorithm. Elitism is done to retain the best chromosomes from the current generation in order to utilise in the next generation 43 . Elitist chromosomes do not participate in crossover operation. This process is done in the selection phase of GA where a portion of chromosomes is removed from the parent population selected for crossover. Thus it retains the best chromosomes and guarantees to produce better solution in the next generation (g + 1) even if the fitness of all chromosomes in the (g + 1) th generation are weaker than its previous generation g since best solutions never lost in this manner. In addition, since floating point gene values are used in this research, the optimisation method used by SpiFoG is termed as elitist floating point genetic algorithm. We have used 20% elitist chromosomes from the total population size 100. Maximising f tunes w and κ and thus the value of e decreases. Note that, decrease in e means SNN learns more precisely from the given training dataset. In this research, rank selection method is used to select fittest chromosomes. SpiFoG uses combination of two different crossover methods on total population to produce better chromosomes in the next generation in an efficient manner. Hybrid crossover includes single point crossover method 44 and another crossover method proposed by Haupt et al. 45 . The advantage of the later crossover method is that it always produce feasible solutions and therefore it is used to produce better solutions in terms of synaptic weights w. The definition of this crossover method is given in Eq. (12). The total genes (i + k) × j + (i + k) × j of each chromosome is divided into two parts; one part is for the genes of w and other part is for the genes of κ . Single point crossover is used upon the values of κ so that the value of delays remain integer even after crossover and to prevent the drastic change in delays. Single point crossover just alters the gene positions by selecting a single crossover point. After the generation of new chromosomes, uniform mutation is carried out to add some diversity to the search space. The definition of uniform mutation is given in Eq. (12).
where r 1 , r 2 , and r 3 are random numbers between [0, 1], P1 p pop , P2 p pop are first and second set of selected parents for crossover respectively, C1 p pop , C2 p pop are two new set of solutions after crossover from first and second parents, lb and ub are the lower bound and upper bound of the population, where mutation is to be carried out, and µ c pop represents the mutated solutions of the crossover solutions. The value of crossover rate and mutation rate was taken as 0.8 (i.e., 80% non-elitist selection of parents from totalpopulation) and 0.1, respectively. Training stops when SpiFoG reaches maximum generation. Algorithm 1 shows the mechanism of the adjustment of w ij and κ ij between input and hidden layer neurons. The same mechanism is followed for the neurons reside in the hidden layer and output layer where there is a change in dimension of weights and delays denoted by w jk and κ jk . Algorithm 2 shows the evaluation of fitness values. www.nature.com/scientificreports/

Conclusion
This paper proposed SpiFoG algorithm to train multilayer feed forward SNN. It uses elitist floating point GA with hybrid crossover method discussed in this paper for fast convergence. The dynamics of the spiking neurons were described by LIF neuron introducing random synaptic delay. The delays are also tuned by the SpiFoG. The dimension of synaptic delay and synaptic weight were kept same to allow only a single synaptic terminal between any two pre-synaptic and post-synaptic neurons. To provide a better dimension in terms of biological plausibility both inhibitory and excitatory neurons were allowed. Moreover, the total simulation time was reduced to 20 ms which is directly responsible for better computational cost. SpiFoG shows best accuracy (test set) in case of Iris dataset which is 97.24% and best accuracy (98.32% training set and 97.92% test set) in case of WBC dataset. Also, best CE in case of both the datasets were observed for SpiFoG. The algorithms where CE could not be computed due to less information about the total simulation time, time step and convergence generation, those were compared with SpiFoG in terms of accuracy only. To summarise it is found that SpiFoG is computationally efficient as well as biologically plausible. However, SpiFoG allows the emission of only a single spike and therefore it can not be applied to SNNs where there it is more important to allow multiple spike times for each individual neuron. It is possible to use other spiking neuron models to work with SpiFoG. For the future work we will use SpiFoG to classify time series data with different neuron model and evaluate the performance. www.nature.com/scientificreports/