An SNN consists of multiple pre-synaptic neurons x₀, …, x_N−1 connected to a post-synaptic neuron z (see Figure 1A). Let T be an experimental time period and δ be the time step. The discrete times are denoted by t = 0, δ, 2δ, …  (t < T). The binary variable x_i,s is set to 1 when a pre-synaptic neuron x_i generates a spike at time s and 0 otherwise. Owing to the discrete time, a spike has a width equal to the time step δ. Hence, each set x = {x_is} represents a spatio-temporal spike pattern (see Figure 1B). The pre-synaptic neuron x_i is connected to the post-synaptic neuron z via a synapse with a synaptic weight W_i ≥ 0 and a conduction delay τ_i ≥ 0. A pre-synaptic spike x_is = 1 arrives at the post-synaptic neuron z after a conduction delay τ_i. After τ_i, the spike evokes an excitatory post-synaptic potential (EPSP) proportional to the synaptic weight W_i. The time course of the EPSP is expressed by a temporal relation function g(Δt) that depends on the temporal difference s−t and the conduction delay τ_i (namely, Δt = s + τ_i − t). For multiple pre-synaptic spikes, the EPSP is the linear sum of the EPSPs evoked by the pre-synaptic spikes. Then, the post-synaptic membrane potential v_t at time t is then given by

The binary variable z_t is set to 1 when the post-synaptic neuron z generates a spike at time t and 0 otherwise. For the parameters τ = {τ_i} and W = {W_i} and a given spatio-temporal spike pattern x = {x_is}, the post-synaptic neuron z is assumed to generate a single spike within the experimental time period 0 ≤ t < T, i.e., ∑tzt=1. Such a constraint is important and is commonly applied in studies of temporal coding (e.g., Bohte et al., 2002) where the output of the SNN in temporal coding is defined as the timing of the post-synaptic spike. The present study initially follows the previous studies, but later proposed an alternative model without the constraint (see section 2.5). This study also assumes that for a spatio-temporal spike pattern x = {x_is}, the probability of generating a post-synaptic spike z_t = 1 at time t is exponentially proportional to the EPSP normalized by a constant Z=∑t′exp(vt′) as follows:

Note that p(z_t = 1|x) > 0 at any time t. In other words, the post-synaptic neuron z can elicit a spike even before the first pre-synaptic spike arrives. In biological neural networks, such spikes are induced by physiological noise. Note also that the normalizing constant Z includes a term vt′ denoting the membrane potential at a future time t′ > t. Therefore, sampling the post-synaptic spike z_t = 1 requires anti-causal information. The alternative model proposed in section 2.5 overcomes this limitation.

This subsection introduces the mixture probabilistic model. Let x be a set of visible binary random variables {x_is}, each following a Bernoulli distribution, and z be a set of latent binary random variables {z_t} following a Multinoulli distribution, i.e., ∑tzt=1. Also let ω = {ω_t} denote the mixture weights, and π = {π_ist} denote the posterior probability of x_is = 1 given z_t = 1. The generative Multinoulli-Bernoulli mixture model p(z, x; ω, π) is then expressed as

Here, the posterior probability π_ist is substituted by sigm(V_ist) using the sigmoid function sigm(u) = (1+exp(−u))⁻¹. In addition, V_ist is modeled by the function W_ig(s + τ_i − t)−v where v is a bias parameter and the mixture weight w_t is assumed to have a constant value e^b. The log-probability is then rewritten as

When the experimental time period T is sufficiently long and the time t is sufficiently distant from the temporal boundaries (0 and T), b^ is independent of the post-synaptic spike timing t; that is, b^=b-∑islog (1+eWig(s+τi-t)-v). The posterior probability of z_t = 1 given x = {x_is} is then equivalent to Equation (2);

The SNN model Equation (2) can be optimized through optimization of the Multinoulli-Bernoulli mixture model Equation (4) under the certain constraints. Let θ be the parameter set and X be a dataset of spatio-temporal spike patterns x. In general, training a generative model involves maximizing the model evidence 𝔼_x~X[log p(x; θ)] (Murphy, 2012). For an arbitrary distribution q(z), the evidence is expressed as

where the evidence lower bound L(θ) equals 𝔼_x~X𝔼_z~q(z)[log p(z, x; θ)−log q(z)] and D_KL(·) is the Kullback-Leibler divergence. The model evidence 𝔼_x~X[log p(x; θ)] can be maximized by the Expectation-Maximization (EM) algorithm. In the E-step, q(z) is substituted by the posterior probability p(z|x; θ^old) of the mixture model with the currently estimated parameter set θ^old. The M-step updates the estimated parameter set θ to maximize the evidence lower bound L(θ). The present study employs a stochastic version of the EM algorithm (Sato, 1999; Nessler et al., 2009). Given a spatio-temporal spike pattern x^, the stochastic EM algorithm first generates a post-synaptic spike z^ by Equation (2). This spike generation corresponds to the E-step of the EM algorithm. Next, the parameters θ = {W, τ} are updated to maximize the evidence lower bound L(θ) using a gradient ascent algorithm. This parameter update corresponds to the M-step of the EM algorithm. From an SNN perspective, the parameter update corresponds to synaptic plasticity with delay learning. Note that, unlike the original EM algorithm, the M-step of the stochastic EM algorithm is not guaranteed to maximize the model evidence 𝔼_x~X[log p(x; θ)] but at least maximizes the evidence lower bound L(θ).

Based on Equation (4), the gradients of the evidence lower bound L(θ) w.r.t. the conduction delay τ_i and synaptic weight W_i are respectively given by;

As Equation (4) and the first term in Equation (5) are independent of the time step δ, the second term in Equation (5) (which does depend on δ) is multiplied by δ: This normalization makes Equation (5) robust to the size of the time step δ. The parameters θ are then updated as follows:

where η is a learning rate. The conditions are similar in a supervised manner, except that the post-synaptic spike z^ is sourced externally.

A spatio-temporal spike pattern x is associated with one of the timings t = 0δ, 1δ, …  depending on the timing t of the generated post-synaptic spike z_t = 1. However, a given dataset X of spatio-temporal spike patterns x should be classifiable into smaller groups. In this study, the decision boundaries between N groups are defined as the n-th N-quantiles (n = 1, 2, …, N − 1) of the generated post-synaptic spike timings t for which z_t = 1. For example, suppose that a dataset X contains 100 spatio-temporal spike patterns x falling into estimated two groups. The spatio-temporal spike patterns x are then sorted by the timings t of the post-synaptic spikes z_t = 1. After sorting, patterns 1–50 are classified into group 1 and the remainder are classified into group 2.

In the SNN introduced above, the Multinoulli distribution p(z|x; τ, W) must be sampled over time t to generate a post-synaptic spike z_t = 1. In other words, whether a post-synaptic spike z_t = 1 is generated at time t anti-causally depends on the EPSP vt′ and post-synaptic spike zt′ at a future time t′ > t. To remove this unrealistic assumption, this subsection slightly modifies the SNN formulation. The new formulation replaces the Multinoulli distribution with multiple independent Bernoulli distributions and introduces the following Bernoulli–Bernoulli mixture model with posterior probability ρ = {ρ_ist} of x_is = 1 given z_t = 0;

In this case, given a spatio-temporal spike pattern x, the posterior probabilities p(z_t|x; ω, π, ρ) of the post-synaptic spikes z_t = 1 at times t are completely independent. As described above, the prior probability ω_t of z_t = 1 is replaced with a time-independent term e^b and the posterior probability π_ist of x_is = 1 given z_t = 1 is modeled as π_ist = sigm(V_ist) = sigm(W_ig(s + τ_i − t)−v). In addition, the posterior probability ρ_ist is substituted by a constant sigm(v). The log-probability function becomes

where b^=log ωt1-ωt+∑islog (1+e-v) and ĉ=log (1-ωt)-v∑isxis+∑islog (1+e-v). The prior probability ω_t is replaced by an independent variable b^ and a variable ĉ that depends on b^ and v.

The posterior probability of z_t = 1 given x is expressed as

where b~=b^-∑islog (1+eWig(s+τi-t)-v). The variable b^ is replaced by a variable b~, which is treated as an independent parameter hereafter. The independent parameter b in the Bernoulli–Bernoulli mixture model is additional to the parameters in the Multinoulli-Bernoulli mixture model.

The evidence lower bound L(θ) is given by 𝔼_x~X𝔼_zt~q(zt)[log p(z_t, x; θ)−log q(z_t)]]. The gradients of L(θ) w.r.t. the parameters τ_i and W_i are, respectively given by

Equations (7) and (8) are equivalent to Equations (4) and (5) respectively, but are summed over time t (they sum the post-synaptic spikes z_t = 1). Therefore, given a spatio-temporal spike pattern x and a single post-synaptic spike z_t = 1, the parameters are updated as in the Multinoulli-Bernoulli mixture model.

Unlike the SNN based on the Multinoulli-Bernoulli mixture model, the above SNN is unconstrained by the number of post-synaptic spikes (i.e., ∑tzt∈ℤ). To prevent an excessive number or the complete absence of post-synaptic spikes, this subsection introduces homeostatic plasticity (Turrigiano et al., 1999; Turrigiano and Nelson, 2000), a biological mechanism that presumably maintains the excitability of a neuron within a regular range (Van Rossum et al., 2000; Abraham, 2008; Watt and Desai, 2010; Matsubara and Uehara, 2016). The additional parameter b~ in Equation (6) embodies the intrinsic excitability of the Bernoulli neuron z. The intrinsic excitability b~ is adjusted after every receipt of the spatio-temporal spike pattern x. When post-synaptic neuron z generates at least one spike (i.e., ∑tzt>0), its intrinsic excitability b~ reduces by b~-; conversely if no spikes are generated (i.e., ∑tzt=0), its intrinsic excitability b~ increases by b~+. Algorithmically, this conditional expression is given by

The spatio-temporal spike patterns x are classified as described in section 2.4, with an important difference; each post-synaptic spike z_t = 1 is weighted by the inverse number of generated post-synaptic spikes, i.e., 1/∑tzt. In other words, a spatio-temporal spike pattern x is associated with a group by majority voting of multiple post-synaptic spikes z_t = 1.

In this study, the temporal relation function g(Δt) was expressed as the following Gaussian function:

Here, the parameters μ and σ denote the mean and standard deviation of the distribution, respectively. Figure 2A plots the temporal relation function g(Δt) when μ = 1.5 ms and σ = 1 ms. The EPSP is commonly modeled by the double-exponential function (e.g., Shouval et al., 2010), which has zero gradient for all Δt < 0. As the proposed SNN is trained according to the gradient of the temporal relation function g(Δt), the double-exponential function is unsuitable for the purpose. In this study, the causality was preserved by clamping the temporal relation function g(s + τ_i − t) to 0 for all s − t < 0.

The experimental time period T was 50 ms, the time step δ was set to 0.05 ms. The remaining parameters were set to μ = 1.5 ms, σ = 1.0 ms, v = 10, and η = 0.001, unless otherwise stated. The gradients of the evidence lower bound L(θ) w.r.t. the conduction delay τ_i and synaptic weight W_i were considered to be changes by the delay learning and synaptic plasticity, respectively. When plotted against the temporal difference Δt = s+τ_i−t, the gradients represent the temporal windows of the delay learning and synaptic plasticity. As shown in Figure 2B, the conduction delay τ_i increases (decreases) when the temporal difference Δt is smaller (larger) than 1.5 ms. The EPSP peaks at 1.5 ms, indicating that a post-synaptic spike x_t = 1 attracts the nearby pre-synaptic spikes x_is = 1. The temporal window of the synaptic weight W_i is similar to that of the STDP (Figure 2C). Specifically, the synapse is potentiated (depressed) when the temporal difference Δt is positive (negative), but is always depressed when the temporal difference Δt is large and positive, as noted in previous electrophysiological studies (Markram et al., 1997; Bi and Poo, 1998) and theoretical studies (Shouval et al., 2002; Wittenberg and Wang, 2006; Shouval et al., 2010). Figure 2D plots the amount of synaptic modification ΔW_i versus the current synaptic weight W_i. The decreasing trend is again consistent with previous electrophysiological (Bi and Poo, 1998) and theoretical (Van Rossum et al., 2000; Gütig et al., 2003; Shouval et al., 2010; Matsubara and Uehara, 2016) studies.

First, the SNN and the proposed learning algorithm based on the Multinoulli-Bernoulli mixture model were evaluated on toy spike patterns generated by three pre-synaptic neurons. Figure 3A shows typical spatio-temporal spike patterns x. The pre-synaptic neurons x₀, x₁, and x₂ elicited spikes at 1+ξ₀, 5+ξ₁, and 13+ξ₂ ms respectively in spike pattern A, and at 13+ξ₀, 9+ξ₁, and 1+ξ₂ ms, respectively in spike pattern B. Here, ξ₀, ξ₁, and ξ₂ are the noise terms following a uniform distribution U(-1,1). Fifty samples from each of spike patterns A and B were extracted as the training dataset; other fifty samples were reserved as the test dataset. All synaptic weights W_i were uniformly initialized to 1, and the conduction delays τ_i (ms) were sampled from the uniform distribution U(5,15). The post-synaptic neuron z was repeatedly fed with the spatio-temporal spike patterns x, and generated spikes z_t = 1. The gray shaded areas in Figure 3A depicts the posterior probability p(z_t = 1|x; W, τ) of spiking each millisecond. Figure 3B shows the probability distribution of post-synaptic spike timings z_t = 1 for the test dataset before learning. The unlearned model cannot discriminate the spike patterns A and B using the post-synaptic spike timings z_t = 1. Note that the proposed SNN is intrinsically probabilistic, meaning that a post-synaptic neuron can elicit a spike at an arbitrary time t. During the learning procedure, the synaptic weights W_i and the conduction delays τ_i were gradually changed by the proposed learning algorithm (Figure 3E). The conduction delays τ_i (ms) were clamped to the range 0–20 ms but never reached the limits. Figure 3F shows random samples of the post-synaptic spike timings z_t = 1, separated into spike patterns A and B by the decision boundary (the solid black line). After sufficiently many samples, the spike patterns converged into limited temporal ranges and the two groups were clearly demarcated. The post-synaptic spike timings z_t = 1 and their distributions after the learning procedure are exemplified in Figures 3C,D, respectively. In this trial, the classification accuracy of both the training and test datasets converged to 100 %. The average classification accuracy over 100 trials was 99.6 ± 2.6% for the training dataset and 99.6 ± 2.5% for the test dataset (Figure 3G). The results are summarized in Table 2.

Next, the proposed learning algorithm was evaluated on the iris flower dataset (Fisher, 1936), which consists of the data of three species (Setosa, Versicolor, and Virginica). The dataset of each species comprises 50 data points, and each data point consists of four features (the lengths and widths of the sepals and petals). Accordingly, the number N of pre-synaptic neurons was set to 4. The spike timing of each feature was normalized to the range 0–10 ms. As the sepal length k₀ ranges from 4.3 to 7.9 cm, the pre-synaptic neuron x₀ generated a single spike at s0=10×k0-4.37.9-4.3 ms. The generated post-synaptic spikes z_t = 1 were then classified into three groups corresponding to the three species. Fifteen randomly chosen data points were assigned as the test set; the remainder were used as the training dataset. The other conditions were those described in section 3.2. The results are summarized in Figure 4. The average classification accuracy over 100 trials was 89.5 ± 5.3% for the training dataset and 89.5 ± 7.5% for the test dataset.

For comparison, this subsection introduces the results of the multilayer ReSuMe algorithm (Sporea and Grüning, 2013), a supervised learning algorithm for multilayer SNNs in rate coding. The multilayer ReSuMe algorithm was also evaluated on the iris flower dataset. In this evaluation, a trial was deemed successful if the classification accuracy exceeded 95% on the training dataset. After weeding out the unsuccessful trials, the proposed learning algorithm and multilayer ReSuMe algorithm achieved a classification accuracy of 94.7 ± 6.5 and 94.0 ± 0.79% respectively, on the test dataset (see Table 3). Despite its single-layer architecture, the SNN trained by the proposed learning algorithm classified the iris flower dataset at least as accurately as the multilayer SNN trained by the multilayer ReSuMe algorithm.

In this subsection, the proposed learning algorithm was evaluated on the MNIST dataset (LeCun et al., 1998), which contains 28 × 28 grayscale images of 70,000 handwritten digits. In standard rate coding, the pixel intensities are represented by 28 × 28 = 784 pre-synaptic neurons (Nessler et al., 2009; Beyeler et al., 2013; O'Connor et al., 2013; Querlioz et al., 2013; Neftci et al., 2014; Diehl and Cook, 2015; Zambrano and Bohte, 2016). To emphasize the temporal coding, this study instead represented the rows and columns of the images by 28 pre-synaptic neurons and their 28 corresponding spike timings, respectively. When a pre-synaptic neuron x_i generated a spike at time s, the intensity of the pixel in row i and column s of the image was considered to exceed 0.5 (see Figure 5A). For simplicity, only the digits 0 and 8 were included in the analysis. After removing the other digits, 13,728 images were available. This experiment employed a 10-fold cross-validation: 10 % of images were randomly chosen for the test set. The results are summarized in Figure 5. This experiment yielded much more pre-synaptic spikes x_is = 1 than the previous two experiments. During transmission, the pre-synaptic spikes x_is = 1 that largely contributed to the post-synaptic spike z_t = 1 were emphasized. More specifically, if μ−1 < Δt = s+τ_i−t < μ+1 for the post-synaptic spike z_t = 1, the pre-synaptic spikes x_is = 1 was emphasized with a thick color line, whereas other pre-synaptic spikes x_is = 1 were attenuated (see Figures 5A,C). According to Figure 5C, the C-shaped edges are heightened, indicating that the SNN selectively detects and responds to these edges. In general, the C-shaped edge appears on the left side of digit 0 and in the center of digit 8, so the post-synaptic spikes z_t = 1 responded earlier to “0” than to “8”. In some trials, the SNN responded to reversed C-shaped edges (or to edges in digits such as “7”), and the post-synaptic spikes z_t = 1 responded earlier to “8” than to “0” (see Figure 6). The average classification accuracy over 100 trials was 88.7 ± 3.4% for the training dataset and 88.7 ± 3.7% for the test dataset.

The proposed learning algorithm was also evaluated without delay learning. More specifically, the synaptic weights W_i were updated by Equation (5), and the conduction delays τ_i were clamped to their initial values (i.e., were not updated by Equation 4). The other conditions were those described in previous subsections. All of the accuracies were drastically reduced (see Table 4).

The proposed learning algorithm is adaptable to supervised learning, in which synaptic plasticity and delay learning are not self-driven by generated post-synaptic spikes, but are driven by external spikes z~t=1 generated at the desired timings (Bohte et al., 2002; Ponulak, 2005; Gütig and Sompolinsky, 2006; Pfister et al., 2006; Paugam-Moisy et al., 2008; Taherkhani et al., 2015; Matsubara and Torikai, 2016). However, the desired timings z~t=1 are generally unknown. In this study, the spatio-temporal spike patterns x are classified by whether the post-synaptic spike timing z_t = 1 arrives earlier or later than the boundary decision. The group with the latest average post-synaptic spike timing should generate the most delayed post-synaptic spike. Therefore, the “late” group were given an external post-synaptic spike z~t+d=1 (where d > 0) after the generated post-synaptic spike z_t = 1. Meanwhile, the “early” group were given a post-synaptic spike z~t-d=1 before the generated post-synaptic spike z_t = 1. These post-synaptic spikes z_t+d = 1 and z_t−d = 1 tend to delay or advance the generated post-synaptic spike z_t = 1, respectively. In this study, the time parameter d was set to δ = 0.05 ms. The results of supervised learning are summarized in Tables 2, 3. In almost every case, the classification accuracy was improved, confirming that the proposed learning algorithm can be adapted to supervised learning without manually adjusting the spike timings.

This subsection evaluates the SNN based on the Bernoulli–Bernoulli mixture model described section 2.5. The homeostatic plasticity parameters were set to b₊ = 0.01 and b₋ = 0.0001, indicating that, on average, one in every hundred spatio-temporal spike patterns induces no post-synaptic spike z_t = 1. To determine the decision boundaries, the absence of the post-synaptic spike z_t = 1 was considered to indicate an exceptionally late post-synaptic spike (e.g., z_T = 1). As shown in Tables 3, 5, the classification accuracy was improved in the iris flower dataset but degraded in the MNIST dataset. As an example, Figure 7 shows the results of applying multiple post-synaptic spikes in the MNIST dataset.

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.