A Radial Basis Function Spike Model for Indirect Learning via Integrate-and-Fire Sampling and Reconstruction Techniques

Abstract

This paper presents a deterministic and adaptive spike model derived from radial basis functions and a leaky integrate-and-fire sampler developed for training spiking neural networks without direct weight manipulation. Several algorithms have been proposed for training spiking neural networks through biologically-plausible learning mechanisms, such as spike-timing-dependent synaptic plasticity and Hebbian plasticity. These algorithms typically rely on the ability to update the synaptic strengths, or weights, directly, through a weight update rule in which the weight increment can be decided and implemented based on the training equations. However, in several potential applications of adaptive spiking neural networks, including neuroprosthetic devices and CMOS/memristor nanoscale neuromorphic chips, the weights cannot be manipulated directly and, instead, tend to change over time by virtue of the pre- and postsynaptic neural activity. This paper presents an indirect learning method that induces changes in the synaptic weights by modulating spike-timing-dependent plasticity by means of controlled input spike trains. In place of the weights, the algorithm manipulates the input spike trains used to stimulate the input neurons by determining a sequence of spike timings that minimize a desired objective function and, indirectly, induce the desired synaptic plasticity in the network.

1. Introduction

This paper presents a deterministic and adaptive spike model obtained from radial basis functions (RBFs) and a leaky integrate-and-fire (LIF) sampler for the purpose of training spiking neural networks (SNNs), without directly manipulating the synaptic weights. Spiking neural networks are computational models of biological neurons comprised of systems of differential equations that can reproduce some of the spike patterns and dynamics observed in real neuronal networks [1, 2]. Recently, SNNs have also been shown capable of simulating sigmoidal artificial neural networks (ANNs) and of solving small-dimensional nonlinear function approximation problems through reinforcement learning [3–5]. Like all ANN learning techniques, existing SNN training algorithms rely on the direct manipulation of the synaptic weights [4–9]. In other words, the learning algorithms typically include a weight-update rule by which the synaptic weights are updated over several iterations, based on the reinforcement signal or network performance.

In many potential SNN applications, including neuroprosthetic devices, light-sensitive neuronal networks grown in vitro, and CMOS/memristor nanoscale neuromorphic chips [10], the synaptic weights cannot be updated directly by the learning algorithm. In several of these applications, the objective is to stimulate a network of biological or artificial spiking neurons to perform a complex function, such as processing an auditory signal or restoring a cognitive function. In neuroprosthetic medical implants, for example, the artificial device may consist of a microelectrode array or integrated circuit that stimulates biological neurons via spike trains. Therefore, the device is not capable of directly modifying the synaptic efficacies of the biological neurons, as do existing SNN training algorithms, but it is capable of stimulating a subset of neurons through controlled pulses of electrical current.

As another example, light-sensitive neuronal networks grown in vitro can be similarly stimulated through controlled light patterns that cause selected neurons to fire at precise moments in time, in an attempt to induce plasticity, while their output is being recorded in real time using a multielectrode array (MEA) [11]. In this case, a digital computer can be used to determine the desired stimulation patterns for an in-vitro neuronal network with random connectivity, produced by culturing dissociated cortical neurons derived from embryonic day E18 rat brain [11, 12]. As a result, the cultures may be used to verify biophysical models of the mechanisms by which biological neuronal networks execute the control and storage of information via temporal coding and learn to solve complex tasks over time. In these networks, the actual connectivity and synaptic plasticities are typically unknown and cannot be manipulated directly as required by existing SNN learning algorithms.

This paper presents a novel indirect learning approach and algorithm that assume synaptic weights cannot be updated or manipulated at any time. The learning algorithm induces changes in the SNN weights by modulating spike-timing dependent plasticity (STDP) through controlled input spike trains. The algorithm adapts the input spike trains that are used to stimulate the input neurons of the SNN and, thus, are realizable through controlled pulses of electric voltage or controlled pulses of blue light. The main difficulty to be overcome in indirect learning is that the algorithm aims at adapting pulse signals, such as square waves, in place of continuous-valued weights. While available in closed analytic form, these signals typically are represented by piece-wise continuous, multi-valued (or many-to-one), and nondifferentiable functions that are difficult to adapt or update using optimization or reinforcement-learning algorithms. Furthermore, stimulation patterns typically are generated by spike models that are stochastic, such as the Poisson spike model [5, 13–15]. Thus, even when the spike model is optimized, it does not allow for precise timing of pre- and post-synaptic firings, and as a result, may induce undesirable changes in the synaptic weights.

In this paper, a deterministic spike model that allows for precise timing of neuron firings is obtained using adaptive RBFs to model the characteristics of the spike pattern through a continuous and infinitely differentiable function that also is one to one. The RBF model is combined with an LIF sampling technique originally developed in [16, 17] for the approximate reconstruction of bandlimited functions. It is shown that, by this approach, the spike trains generated by the LIF sampler display the precise characteristics specified by the RBF model. Furthermore, this deterministic spike model can be optimized to modulate STDP through controlled input spike trains that bring about the desired SNN weight change without direct manipulation. The indirect learning approach presented in this paper is applicable both in supervised and unsupervised settings. In fact, when the desired SNN output is unknown, it can be replaced by a reinforcement signal produced by a critic SNN, as shown by the adaptive critic method reviewed in Section 3.

The paper is organized as follows. The model of spiking neural network used to derive and demonstrate the training equations is presented in Section 2. In Section 3, an adaptive critic approach is described to illustrate how the proposed indirect training methodology can be applied using reinforcement learning, for example, to model or control a dynamical plant. The novel spike model and indirect training methodology are presented in Section 4 and demonstrated on a benchmark problem involving a two-node spiking neural network. The generalized form of gradient equations for indirect training is presented in Section 5 and demonstrated through a three-node spiking neural network. These gradient equations show how the methodology can be generalized to any spiking neural network with the characteristics described in the next section.

2. Spiking Neural Network Model

2.1. Models of Neuron and Synapse

Various models of SNNs have been proposed in recent years, motivated by biological studies that have shown complex spike patterns and dynamics to be an essential component of information processing and learning in the brain. The two crucial considerations involved in determining a suitable SNN model are the range of neurocomputational behaviors the model can reproduce, and its computational efficiency [18]. As can be expected, the implementation efficiency typically increases with the number of features and behaviors that can be accurately reproduced [18], such that each model offers a tradeoff between these competing objectives. One of the computational neuron models that is most biophysically accurate is the well-known Hodgkin-Huxley (HH) model [2]. Due to its extremely low computational efficiency, however, using the HH model to simulate large networks of neurons can be computationally prohibitive [18]. Recently, bifurcation studies have been used to reduce the HH dynamics from four to two differential equations, referred to as the Izhikevich model, which are capable of reproducing a wide range of spiking patterns and behaviors with much higher efficiency than the HH model [19].

In [15], the authors proposed an indirect training method based on a Poison spike model and demonstrated it on a network of Izhikevich neurons. In this work it was found that the adaptive critic architecture described in Section 3 could be implemented without modeling the neuron response in closed form. However, the effectiveness of the approach in [15] was limited in that, due to the use of a stochastic Poison model, the Izhikevich SNN could not converge to the optimal control law. Therefore, in this paper, a new deterministic spike model is proposed and implemented by deriving the training equations in closed form, using a leaky integrate-and-fire (LIF) SNN. The LIF is the simplest model of spiking neuron. It has the advantages that it displays the highest computational efficiency and is amenable to mathematical analysis [13, 14].

The LIF membrane potential, 𝑣(𝑡), is governed by 𝐶_𝑚𝑑𝑣(𝑡)𝑑𝑡=𝐼_leak(𝑡)+𝐼_𝑠(𝑡)+𝐼_inj(𝑡), (1) where 𝐶_𝑚 is the membrane capacitance, 𝐼_leak(𝑡) is the current due to the leak of the membrane, 𝐼_𝑠(𝑡) is the synaptic input to the neuron, and 𝐼_inj(𝑡) is the current injected to the neuron [20]. The leak current is defined as follow: 𝐼_leak(𝑡)𝐶=−^𝑚𝜏_𝑚𝑣(𝑡)−𝑉₀, (2) where 𝑉₀ is the resting potential and 𝜏_𝑚 is the passive-membrane time constant [20]. 𝜏_𝑚 is related to the capacitance, 𝐶_𝑚, and the membrane resistance, 𝑅_𝑚, of the membrane potential by 𝜏_𝑚=𝑅_𝑚𝐶_𝑚.

The response of the membrane potential is obtained by solving the differential equation (1) analytically for 𝑣(𝑡), such that 𝑣(𝑡)=𝑉₀+𝑒^{−𝑡/𝜏𝑚𝑡}_𝑡0𝐼_inj𝜌𝐶_𝑚𝑒^{𝜌/𝜏𝑚}𝑑𝜌, (3) where 𝜌 is a dummy variable used for integration and 𝑡₀ is the time at which the membrane potential equals 𝑉₀ [20]. Whenever the membrane potential, 𝑣(𝑡), reaches a prescribed threshold value, 𝑉_th, the neuron will fire. At this point, we have considered an isolated neuron that is stimulated by an external current, 𝐼_inj(𝑡). When the LIF neuron (1) is part of a larger network, the input current, referred to as the synaptic current, is generated by the activity of presynaptic neurons.

Let the synaptic current, denoted by 𝐼_𝑠(𝑡), be modeled as the sum of the currents of excitatory presynaptic neurons and of inhibitory presynaptic neurons, 𝐼_𝑠(𝑡)=𝐶_𝑚^𝑁𝐸_𝑘=1𝑎_𝐸,𝑘𝑆_𝐸,𝑘(𝑡)+𝐶_𝑚^𝑁𝐼_𝑘=1𝑎_𝐼,𝑘𝑆_𝐼,𝑘(𝑡), (4) where the subscript 𝐸 denotes inputs from excitatory neurons, while the subscript 𝐼 denotes inputs from inhibitory neurons. The amplitudes, 𝑎_𝐸,𝑘>0 and 𝑎_𝐼,𝑘<0, represent the change in potential due to a single synaptic event and depend on the weight of the synapse. 𝑁_𝐸 and 𝑁_𝐼 are the numbers of excitatory and inhibitory current synapses, respectively. 𝑆_𝐸,𝑘 and 𝑆_𝐼,𝑘 describe the excitatory and inhibitory synaptic inputs as a series of input spikes to each synapse. The synaptic inputs are modeled as a sum of instantaneous impulse functions, 𝑆_𝐸,𝑘(𝑡)=_{𝑡𝐸,𝑘}𝛿𝑡−𝑡_𝐸,𝑘,𝑆_𝐼,𝑘(𝑡)=_{𝑡𝐼,𝑘}𝛿𝑡−𝑡_𝐼,𝑘, (5)(6) where 𝑡_𝐸,𝑘 and 𝑡_𝐼,𝑘 are the firing times of the presynaptic neurons and 𝛿 represents the Dirac delta function [21].

In addition to the above deterministic properties, there are prevalent stochastic qualities in biological neural networks caused by effects such as thermal noise and random variations in neurotransmitters. For simplicity, these effects are assumed negligible in this paper. However, the reader is referred to [15] for a technique that can be used to incorporate these effects in the above SNN model.

2.2. Model of Spike-Timing-Dependent Plasticity (STDP)

A persistent learning mechanism known as spike-timing-dependent plasticity, recently observed in biological neuronal networks, is used in this paper to model synaptic plasticity in the LIF SNN. Synaptic plasticity refers to the mechanism by which the synaptic efficacies or strengths between neurons are modified over time, typically as a result of the neuronal activity. These changes are known to be driven in part by the correlated activity of adjacently connected neurons. The directions and magnitudes of the changes are dependent on the relative timings of the presynaptic spike arrivals and postsynaptic firings. For simplicity, in this paper, all changes in synaptic strengths are assumed to occur solely as a result of the spike-timing-dependent plasticity mechanism. The approach presented in [5] can be used to also incorporate a model of the Hebbian plasticity.

Spike-timing-dependent plasticity (STDP) is known to modify the synaptic strengths according to the relative timing of the output and input action potentials, or spikes, of a particular neuron. If the presynaptic neuron fires shortly before the postsynaptic neuron, the strength of the connection will be increased. In contrast, if the presynaptic neuron fires after the postsynaptic neuron, the strength of the connection will be decreased, as illustrated in Figure 1. Two constants 𝜏₊ and 𝜏₋ determine the ranges of the presynaptic to postsynaptic interspike intervals over which synaptic strengthening and weakening occur. Let the synaptic efficacy or strength be referred to as weight and denoted by 𝑤. 𝐴 is the maximum amplitude of the change of weight due to a pair of spikes. 𝑡_pre and 𝑡_post denote firing times of the presynaptic neuron and the postsynaptic neuron, respectively. Then, for each set of neighboring spikes, the weight adjustment is given by, Δ𝑤=𝐴𝑒^{[(𝑡pre−𝑡post)𝜏]}, (7) such that the connection weight increases when the postsynaptic spike follows the presynaptic spike and the weight decreases when the opposite occurs. The amplitude of the adjustment Δ𝑤 lessens as the time between the spikes becomes larger, as is illustrated in Figure 1.

Figure 1 

STDP term as a function of the time delay between the last spike of postsynaptic neuron and presynaptic neuron.

Different methods can be used to identify the spikes that give rise to the STDP mechanism. In this paper, the nearest-spike STDP model [22] is adopted, by which, for every spike of the presynaptic neuron, its nearest postsynaptic spike is used to calculate the timing difference in (7), regardless of whether it takes place before or after the presynaptic firing. Then, from (7), the weight change due to one of the spikes of the 𝑖th neuron is modeled by the ruleΔ𝑤_𝑖Δ𝑡_𝑖=𝐴₊𝑒^{Δ𝑡𝑖/𝜏+}ifΔ𝑡_𝑖𝐴≤0₋𝑒^{−Δ𝑡𝑖/𝜏−}ifΔ𝑡_𝑖>0, (8) where Δ𝑡_𝑖=𝑡_1,𝑖−argmin_𝑡2,𝑘𝑡_1,𝑖−𝑡_2,𝑘. (9)Δ𝑡_𝑖 is the firing time difference between the two neurons, 𝑡_1,𝑖 is the firing time of the presynaptic neuron, and 𝑡_2,𝑖 is the firing time of the postsynaptic neuron. The constants used in this paper are 𝜏₊=𝜏₋=20 ms, 𝐴₊=0.006, and 𝐴₋=1.5𝐴₊, where 𝐴₊ and 𝐴₋ determine the maximum amounts of synaptic modification that occur when Δ𝑡 is approximately zero [21].

From (8), the firing time of the postsynaptic neuron is chosen such that the absolute firing time difference between the presynaptic neuron and the postsynaptic neuron, |𝑡_1,𝑖−𝑡_2,𝑘|, is minimized. Therefore, the postsynaptic spike 𝑡_2,𝑘 that is nearest to the presynaptic firing time 𝑡_1,𝑖 is used to calculate the firing time difference Δ𝑡_𝑖. The value of 𝑡_2,𝑘 that minimizes |𝑡_1,𝑖−𝑡_2,𝑘| can be found by comparing the firing time difference of the two neurons. In order to obtain an indirect learning method that does not rely on the direct manipulation of the synaptic weights, in this paper, it is assumed that the weights can only be modified according to the STDP rule (8). Also, the training algorithm cannot specify any of the terms in (8) directly but can only induce the firing times in (8) by stimulating the input neurons, for example, through controlled pulses of electric voltage or blue light.

3. Adaptive Critic Architecture for Indirect Reinforcement Learning

Many approaches have been proposed to train SNNs by modifying the synaptic weights by means of reward signals or by update rules inspired by reward-driven Hebbian learning and modulation of STDP [5, 8]. However, these methods are not applicable to approaches involving biological neuronal networks (e.g., neuronal cultures grown in vitro), because biological synaptic efficacies cannot be manipulated directly by the training algorithm. Similarly, in nanoscale neuromorphic chips, the CMOS/memristor synaptic weights cannot be manipulated directly but may be modified via controlled programming voltages that induce a mechanism analogous to STDP, as demonstrated in [10]. This paper presents an approach for training an SNN through controlled input pulse signals (e.g., voltages or blue light) that are generated by a new deterministic and adaptive spike model presented in Section 4.1.

The derivation of the training equations used to adapt the proposed spike model and subsequent pulse signals are demonstrated in Section 4.2 using a two-node LIF SNN. It is shown that the synaptic weight can be updated by the proposed indirect training method and driven precisely to a desired value, by minimizing a function of the synaptic weight error. It follows that, using a simple chain rule, the same training equations can be used to minimize a function of the decoded SNN output error, as is typically required by supervised training. In this section, we show how the indirect training method can be used for reinforcement learning, using a critic SNN to modulate the STDP mechanisms in an action SNN with synaptic weights that cannot be manipulated directly (e.g., a biological neuronal network).

Let 𝑁𝑁_𝑎 represent an action SNN of LIF neurons, exhibiting STDP and modeled as described in Section 2. It is assumed that the weights of the action SNN cannot be manipulated directly, and the only controls available to the learning algorithm are the training signals (Figure 2), comprised of programming voltages that can be delivered using a square wave or the Rademacher function (Section 4.1). Now, let 𝑁𝑁_𝑐 denote a critic SNN of feedforward fully connected HH neurons. In this architecture, schematized in Figure 2, 𝑁𝑁_𝑎 is treated as a biological network and 𝑁𝑁_𝑐 is treated as an artificial network implemented on a computer or integrated circuit. Thus, the synaptic weights of 𝑁𝑁_𝑐, defined as 𝑤_𝑖𝑗, are directly adjustable, while the synaptic efficacies of 𝑁𝑁_𝑎 can only be modified through a simulated STDP mechanism which is modulated by the input spikes from 𝑁𝑁_𝑐. The algorithm presented in this section trains the network by changing the values of 𝑤_𝑖𝑗 inside 𝑁𝑁_𝑐, which are assumed to be bounded by a positive constant 𝑤_max such that −𝑤_max≤𝑤_𝑖𝑗≤𝑤_max, for all 𝑖,𝑗.

Figure 2 

Adaptive critic architecture comprised of a critic network, 𝑁𝑁𝑐, and, an action network, 𝑁𝑁𝑎, to control a dynamical system.

In this paper, spike frequencies are used to code continuous signals, and the leaky integrator 𝑢=𝛼_{𝑡𝑘∈𝑆𝑖(𝑇)}𝑒^{𝛽(𝑡𝑘−𝑡)}𝐻𝑡−𝑡_𝑘−𝛾 (10) is used to decode spike trains and convert them into continuous signals as required according to the architecture in Figure 2 Where, 𝛼, 𝛽, and 𝛾 are user-specified constants, 𝐻(⋅) is the Heaviside function, and 𝑆_𝑖(𝑇) is the set of spiking times of the output neuron. One advantage of the leaky integrator decoder is its effectiveness of filtering inevitable noise during the flight control without influencing the speed of matching the continuous value with the target function.

3.1. Adaptive Critic Recurrence Relations

Typically, adaptive critic algorithms are used to update the actor and critic weights by computing a synaptic weight increment through an optimization-based algorithm, such as backpropagation ([23], page 359). Therefore a new approach is required in order to apply adaptive critics to SNNs in which changes in the synaptic weights can occur only as a result of pre- and postsynaptic neuronal activity. The training approach presented in Section 4 can be combined with the policy and value-iteration procedures described in this section to supervise stimulation patterns in the action SNN such that the synaptic strengths, and subsequently their dynamic mappings, are adapted according to the STDP rule in (8).

Suppose the action SNN is being trained to control or model a dynamical system which obeys a nonlinear differential equation of the form𝑥̇(𝑡)=𝑓[𝑥(𝑡),𝑢(𝑡),𝑡], (11) where 𝑥∈𝑋⊂ℝ^𝑛 and 𝑢∈𝑈⊂ℝ^𝑚 denote the dynamical system state and control inputs, respectively, and 𝑥̇ denotes the derivative of 𝑥 with respect to time. The model of the dynamical system (11) may be unknown or imperfect and may be improved upon over time by a system identification (ID) algorithm. The macroscopic behavioral goals of the plant can be expressed by the value function or cost-to-go𝑉^𝜋𝑥𝑡_𝑘=^𝑡𝑓−1_{𝑡𝑘=𝑡}ℒ𝑥𝑡_𝑘,𝑢𝑡_𝑘,𝑥𝑡_𝑘+1, (12) where 𝑢(𝑡_𝑘)=𝜋[𝑥(𝑡_𝑘)] is the unknown control law or dynamic mapping to be learned by the action SNN, 𝑁𝑁_𝑎. Then, the cost-to-go in (12) can be used to represent the future performance of the action network and dynamical system, as is accrued from the present, 𝑡_𝑘, up to the final time, 𝑡_𝑓, if subject to the present control law 𝜋[⋅]. The Lagrangian ℒ[⋅] represents instantaneous behavioral goals as a function of 𝑥 and 𝑢.

Since the dynamic mapping 𝜋[⋅] must be adapted over time through the learning algorithm and an accurate plant model may not be available, the cost-to-go typically is unknown and is learned by the critic, 𝑁𝑁_𝑐. Then, by Bellman’s principle of optimality [24], at any time 𝑡_𝑘 the cost-to-go can be minimized online with respect to the present control law, based on the known value of 𝑥(𝑡_𝑘) and predictions of 𝑥(𝑡_𝑘+1). The value of 𝑥(𝑡_𝑘) is assumed to be fully observable from the dynamical system at any present time 𝑡_𝑘, and 𝑥(𝑡_𝑘+1) is predicted or estimated from the approximation of the system’s dynamics (11), based on the present values of the state and the control. In [25], Howard showed that if iterative approximations of the control law and optimal cost-to-go, denoted by 𝜋_ℓ and 𝑉_ℓ, respectively, are modified by a policy-improvement routine and a value-determination operation, respectively, they eventually converge to their optimal counterparts 𝜋^∗ and 𝑉^∗.

The policy-improvement routine states that, given a cost-to-go function 𝑉(⋅,𝜋_ℓ) corresponding to a control law 𝜋_ℓ, an improved control law 𝜋_ℓ+1 can be obtained as follows: 𝜋_ℓ+1𝑥𝑡_𝑘=argmin_𝑢∈𝑈ℒ𝑥𝑡_𝑘,𝑢𝑡_𝑘,𝑥𝑡_𝑘𝑥𝑡+1+𝑉_𝑘+1,𝜋_ℓ, (13) such that 𝑉[𝑥(𝑡_𝑘),𝜋_ℓ+1]≤𝑉[𝑥(𝑡_𝑘),𝜋_ℓ], for any 𝑥(𝑡_𝑘). Furthermore, the sequence of functions 𝐶={𝜋_ℓ∣ℓ=0,1,2,…} converges to the optimal control law 𝜋^∗. The value-determination operation states that given a control law 𝜋(⋅), the cost-to-go can be updated according to the following rule: 𝑉_ℓ+1𝑥𝑡_𝑘+1,𝜋=ℒ𝑥𝑡_𝑘,𝑢𝑡_𝑘,𝑥𝑡_𝑘+1+𝑉_ell𝑥𝑡_𝑘+1,𝜋, (14) such that the sequence of functions 𝑉={𝑉_ℓ∣ℓ=0,1,2,…} converges to 𝑉^∗.

Then, at every iteration cycle ℓ of the adaptive critic algorithm, the above policy and value-iteration procedures can be used to in tandem to supervise training of the actor and critic SNN; that is, 𝑢𝑡_𝑘⟵𝑁𝑁_𝑎ℓ,𝑥𝑡_𝑘≈𝜋_ℓ𝑥𝑡_𝑘,𝑉^𝜋𝑥𝑡_𝑘⟵𝑁𝑁_𝑐ℓ,𝑥𝑡_𝑘,𝜋≈𝑉_ℓ𝑥𝑡_𝑘,𝜋, (15)(16) respectively where ℓ=𝑀𝑘 and 𝑀 is a positive constant chosen based on the desired frequency of the updates. It can be seen that, at ℓ+1, the desired mappings for 𝑁𝑁_𝑎 and 𝑁𝑁_𝑐, provided by (13) and (14), respectively, are based on the actor and critic outputs at ℓ, and on 𝑥(𝑡_𝑘+1), which can be estimated from an available dynamical system model (11). To accomplish (15) and (16), two error metrics defined as a function of the actual (decoded) actor and critic outputs and of the desired actor and critic outputs (13) and (14) must be minimized by the chosen training algorithm, as explained in the following subsection.

3.2. Action and Critic Network Training Algorithm

The action and critic networks implemented in Figure 2 differ in that while the critic network can be trained by a conventional algorithm that manipulates the synaptic weights directly (e.g., see [5, 15]), the action network must be trained by inducing weight changes indirectly through STDP (8), such that its (decoded) output 𝑢 will closely match 𝑢 in (15), for all 𝑥∈𝑋. Since the weights within 𝑁𝑁_𝑎 cannot be adjusted directly, they are manipulated with training signals from 𝑁𝑁_𝑐. Connections are made between 𝑞 pairs of action/critic neurons which serve as outputs in 𝑁𝑁_𝑐 and inputs for 𝑁𝑁_𝑎 (Figure 2). To provide feedback signals to the critic, connections are also created between 𝑟 pairs of neurons from 𝑁𝑁_𝑎 to 𝑁𝑁_𝑐. The information flow through the networks is illustrated in Figure 2. Since it is not known what training signals will produce the desired results in 𝑁𝑁_𝑎, the critic must also be trained to match 𝑉^𝜋 in (16) for all 𝑥∈𝑋 and 𝑢∈𝑈.

Both updates (15) and (16) can be formulated as follows. Let SNN[⋅] denote an action or critic network comprised of multiple spiking neurons and STDP synapses. Then, as shown in (15) and (16), SNN[⋅] must represent a desired mapping between two vectors 𝜉∈ℝ^𝑝 and 𝑦∈ℝ^𝑟: 𝑦𝑡_𝑘=𝑔_ℓ𝜉𝑡_𝑘⟵SNN𝜉𝑡_𝑘,𝑊_ℓ𝑡_𝑘≜𝑦𝑡_𝑘. (17) The desired mapping 𝑔_ℓ∶ℝ^𝑝→ℝ^𝑟 can be assumed known and stationary during every iteration cycle ℓ and is updated by the policy/value-iteration routine. Let 𝑊_ℓ(𝑡_𝑘)={𝑤_𝑖𝑗(𝑡_𝑘)} denote a matrix containing the SNN synaptic weights at time 𝑡_𝑘. Then, for proper values of 𝑊_ℓ, if the SNN is provided with an observation of the input, 𝜉^𝑙, encoded as 𝑛 spike train sequences 𝑋^𝑙_𝑖={𝑡̂_𝑖,𝜅∶𝜅=1,…,𝑁^𝑙_𝑖}, 𝑖=1,…,𝑛, over a time interval [𝑡_𝑘,𝑡_𝑘+𝜏] with the firing times at 𝑡̂_𝑖,𝜅, the decoded output of the SNN must match the output 𝑦^𝑙=𝑔_ℓ[𝜉^𝑙], where since 𝜏≪(𝑡_𝑘+1−𝑡_𝑘), 𝑋^𝑙_𝑖 can be used to encode instantaneous values of 𝜉^𝑙 at any 𝑡_𝑘. Then, the chosen SNN training algorithm must modify the synaptic weights such that a figure of merit representing distance is optimized at the output space of the mapping in (17).

For the critic network, 𝑁𝑁_𝑐, the synaptic weights are adjusted manually, using the reward-modulated Hebbian approach presented in [5, 15]. Because the target function 𝑦 is known for all values of 𝜉, the SNN output error (𝑦−𝑦) is also known and can be used as a feedback to the critic in the form of an imitated chemical reward, 𝑟(𝑡), that decays over time with time constant 𝜏_𝑐. The critic reward is modeled as 𝑟(𝑡)=𝑏𝑦,𝑦+𝑟(𝑡−Δ𝑡)𝑒^−(𝑡−𝑡̂_𝑖)𝜏_𝑐, (18) where for the critic 𝑦∈ℝ. Therefore, the error function can be defined as 𝑏(𝑦,𝑦)=sgn(𝑦−𝑦), where sgn(⋅) is the signum function. Thus, the value of 𝑏(⋅) is positive when the critic’s output is too low, and it is negative when the critic’s output is too high.

For the action network, 𝑁𝑁_𝑎, the synaptic weights cannot be manipulated directly; therefore the distance between 𝑦 and the (decoded) SNN’s output 𝑦 is minimized with respect to the parameters of the RBF spike model described in the next section. From (17), we identify (tag) two sets of neurons referred to as input and output neurons (Figure 2), where each neuron in the set provides the response for one element of 𝜉 and 𝑦, respectively, in the form of a spike train. It is assumed that the set of input neurons can be induced to fire on command with very high precision over a training time interval [𝑡_𝑘,𝑡_𝑘+𝜏], with 𝜏<𝑇≪(𝑡_𝑘+1−𝑡_𝑘). The input neurons’ firings could be implemented in practice by using local programming voltages with controlled pulse width and height that are easily realizable both in CMOS/memristor chips [10] and in neuronal networks grown in vitro [26]. In order to induce STDP in a manner that will improve the SNN representation of the mapping in (17), the programming voltages are delivered based on an optimized spiking sequence S^𝜄_𝑖={𝑡̂_𝑖,𝜅∶𝜅=1,…,𝑁^𝜄_𝑖}, 𝑖=1,…,𝑞, during the 𝑖th time interval of the training algorithm, [𝑡_𝜄,𝑡_𝜄+1].

Existing spike models [13, 14] cannot be used to generate 𝑆^𝜄_𝑖 because they are stochastic and, as such, they do not allow for precise timing of pre- and postsynaptic firings. As a result, they may induce undesirable changes in the synaptic weights by virtue of the STDP rule (7), illustrated in Figure 1, which shows that a small difference in the arrival time of a pre- and postsynaptic spike can obtain the opposite effect in terms of the weight change Δ𝑤. Beside allowing for precise timing of neuron firings, the new deterministic spike model presented in the next section can be updated by the training algorithm, by optimizing a corresponding continuous signal modeled by a superposition of Gaussian radial basis functions (RBFs), which is characterized by a set of adjustable parameters 𝑃={𝑤_𝑘,𝑐_𝑘,𝛽_𝑘∣𝑘=1,…,𝑁}, where 𝑤_𝑘 is the height the RBF pulse, 𝑐_𝑘 is the center of RBF pulse, 𝛽_𝑘 is the width of the RBF pulse, and 𝑁 is the number of RBF pulses. As shown in the next section, the continuous RBF signal can be integrated against a suitable averaging function in a leaky integrator and then compared to a positive threshold, by means of a leaky integrate-and-fire (LIF) sampler. Subsequently, a precise pulse function with desired widths and intensity can be generated at a sequence of time instants, 𝑡̂_𝑖,𝜅, during the interval [𝑡_𝜄,𝑡_𝜄+1], that correspond to the centers of the RBF signal specified by 𝑃. Thus, a set of optimal RBF parameters 𝑃^∗ used to generate 𝑆^𝜄_𝑖 can be determined by minimizing a measure of the distance between the (decoded) SNN’s output 𝑦 and the desired output 𝑦, computed by the policy-improvement routine (13).

In the next section, the derivation of gradient equations that can be used to minimize the SNN’s output error (𝑦−𝑦) with respect to 𝑃 is illustrated by means of a two-node LIF neural network, with STDP modeled as shown in Section 2. Let 𝑃={𝑝_𝑙𝑘} denote a matrix of RBF parameters to be adjusted by the training algorithm. It can be easily shown that the minimization of an error function 𝐸(𝑦−𝑦) with respect to the elements of 𝑃 when 𝑦 is a known constant can be accomplished using the gradients 𝜕𝑦/𝜕𝑝_𝑙𝑘, or some function thereof, based on the chosen unconstrained minimization algorithm [27]. As shown in (17), for a given input 𝜉^𝑙 the only SNN variables to be optimized are the synaptic weights 𝑤_𝑖𝑗. From the STDP rule (7), the weights are a function of the parameters 𝑝_𝑙𝑘, which determine the input spike sequence (or training signal) to the SNN (Figure 3). By virtue of the chain rule of differentiation, it follows that 𝜕𝑦𝜕𝑝_𝑙𝑘=𝜕𝑦𝜕𝑤_𝑖𝑗𝜕𝑤_𝑖𝑗𝜕𝑝_𝑙𝑘. (19) Since 𝜕𝑦/𝜕𝑤_𝑖𝑗 can be computed from the SNN model (Section 2), it follows that if a training algorithm can successfully modify the synaptic weights 𝑤_𝑖𝑗 using the proposed spike model, it also can successfully modify 𝑦, simply by redefining the error function. In other words, if a training algorithm can successfully train the synaptic weights of an SNN using the gradient 𝜕𝑤_𝑖𝑗/𝜕𝑝_𝑙𝑘, it follows that it can also train the SNN using the gradient 𝜕𝑦/𝜕𝑝_𝑙𝑘, since the component 𝜕𝑦/𝜕𝑤_𝑖𝑗 is known and given by the SNN equations (1)–(5). Based on this property, the novel spike model and indirect training algorithm are illustrated by training the synaptic weight of a two-node LIF SNN to meet a desired value 𝑤^∗, without direct manipulation.

Figure 3 

Model of two-node LIF spiking neural network.

4. Indirect Training Methodology

Consider the two-node LIF SNN schematized in Figure 3, modeled using the approach described in Section 2. 𝑠(𝑡) represents the input given to the LIF sampler, and 𝑆(𝑡) denotes the corresponding LIF sampler's output. Based on the approach presented in Section 3, the SNN synaptic strength 𝑤₂₁, representing the synaptic efficacy for a pre-synaptic neuron (labeled by 𝑖=1) and a postsynaptic neuron (labeled by 𝑖=2) cannot be modified directly by the training algorithm but can only change as a result of the STDP mechanism described in Section 2. In place of controlling 𝑤₂₁, the goal of the indirect training algorithm is to determine a spike train that can be used to stimulate the input neuron (𝑖=1) using 𝐼_inj, thereby inducing it to spike, such that the synaptic weight 𝑤₂₁ changes from an initial (random) value to a desired value 𝑤^∗.

For this purpose, we introduce a continuous spike model comprised of a superposition of Gaussian radial basis functions (RBFs): 𝑠(𝑡)=^𝑁_𝑘=1𝑤_𝑘exp−𝛽_𝑘𝑡−𝑐_𝑘², (20) where 𝑤_𝑘 determines the height of the 𝑘th RBF, which will decide the constant current input 𝐼_inj. 𝛽_𝑘 determines the width of the 𝑘th RBF and 𝑐_𝑘 determines the center of the 𝑘th RBF, where 𝑘=1,…,𝑁 and 𝑁 is the number of RBFs. Thus, the set of RBF adjustable parameters is 𝑃={𝑤_𝑘,𝑐_𝑘,𝛽_𝑘∣𝑘=1,…,𝑁}. A spike train can be obtained from the continuous spike model (20) by processing 𝑠 using a suitable LIF sampler that outputs a square pulse function 𝑆(𝑡) with the same heights, widths, and centers, as the RBF spike model (20).

In this two-node SNN model, there is no synaptic current sent to the first neuron because there are no presynaptic neurons connected to it. Rather, the input for the first neuron is the square pulse function provided by the output of the LIF sampler. Therefore, during one square pulse, the input, 𝐼_inj(𝑡), to neuron 𝑖=1 can be viewed as a constant, which for our case is ℎ, the height of the RBF. For a fixed threshold and a constant input, the time it takes for the first neuron to fire (or a spike to be generated) is 𝑇=−𝜏_𝑚𝜃ln1−ℎ𝑅_𝑚ℎ𝑅_𝑚>𝜃, (21) where 𝜃 is the potential difference between the spike generating threshold, 𝑉_th, and the resting potential, 𝑉₀; that is, 𝜃=𝑉_th−𝑉₀. ℎ is the height of input square pulse, 𝑅_𝑚 is the resistance of the membrane, and 𝜏_𝑚 is the passive-membrane time constant [20]. Thus, from (21), the value of 𝑇 can be controlled by adjusting ℎ. We allow the first neuron to fire near the end of a square pulse by inputting a spike sequence generated by an RBF model with a suitable width and height. Then, the firing times of the first neuron can be easily adjusted by altering the centers of the RBF.

In this simple example, it is assumed that the synapse is of the excitatory type. Therefore, (4) can be written as, 𝐼_𝑠(𝑡)=𝐶_𝑚𝑎_𝐸^𝑘𝑡_𝑘=𝑘0𝛿𝑡−𝑡_1,𝑘−𝜏_𝑑, (22) where 𝜏_𝑑 is a known constant that represents the conduction delay of the synaptic current from neuron 𝑖=1 and 𝑡_1,𝑘 are the firing times of the first neuron. 𝑘₀ is the index of the spike of the first neuron, which occurs after the previous spike of the second neuron, and 𝑘_𝑡 is the index of the spike of the first neuron, which provokes the firing time of the second neuron. The only input to the second neuron is the synaptic current defined in (22). Therefore, substituting (22) and (2) into (1) results in the following equation for the membrane potential, 𝐶_𝑚𝑑𝑣(𝑡)𝐶𝑑𝑡=−^𝑚𝜏_𝑚𝑣(𝑡)−𝑉₀+𝐶_𝑚𝑎_𝐸^𝑘𝑡_𝑘=𝑘0𝛿𝑡−𝑡_1,𝑘−𝜏_𝑑. (23) Solving (23) for the membrane potential of the second neuron provides the response of neuron 𝑖=2, 𝑣(𝑡)=𝑉₀+^𝑘𝑡_𝑘=𝑘0𝑎_𝐸exp−𝑡−𝑡_1,𝑘−𝜏_𝑑𝜏_𝑚𝐻𝑡−𝑡_1,𝑘−𝜏_𝑑, (24) where 𝐻(⋅) is the Heaviside function. Whenever the membrane potential in (24) exceeds the threshold 𝑉_th, the second neuron fires, and the membrane potential 𝑣(𝑡) is then set instantly equal to 𝑉₀. It follows that the firing times of the second neuron 𝑡_2,𝑗 can be written as a function of the firing times of the first neuron, 𝑡_2,𝑗=𝑡_1,𝑘𝑡+𝜏_𝑑𝐻^𝑘𝑡_𝑘=𝑘0𝑎_𝐸−𝑡exp^2,𝑗−𝑡_1,𝑘−𝜏_𝑑𝜏_𝑚×𝐻𝑡_2,𝑗−𝑡_1,𝑘−𝜏_𝑑−𝑉_th−𝑉₀, (25) where 𝑗 is the index of the firing times. In addition, the membrane potentials of the first neuron and the second neuron can be calculated using (24) and (21).

An example of the membrane potential time history for the two neurons is plotted in Figure 4, where the square pulse function 𝑆 used to stimulate the first neuron is plotted along with the neurons’ membrane potentials 𝑣₁ and 𝑣₂. The red dots denote the firing times and potentials for the two neurons. It can be seen from Figure 4 that, with this input, the first neuron fires five times and the second neuron fires one time. In this SNN, the input to the second neuron is given only by the synaptic current caused by the firing of the first neuron. It can be seen that only when the second and third firing times of the first neuron are very close, the membrane potential of the second neuron increases over the threshold, causing it to fire, whereas the fourth and fifth firing times of the first neuron are too far apart to cause firing of the second neuron.

(a)

(b)

(c)

Figure 4 

Membrane potential of the two neurons in the SNN in Figure 3.

4.1. Deterministic Spike Model

The deterministic spike model consists of a continuous RBF model in the form (20) combined with an LIF sampler that converts the RBF into a square wave function. The LIF sampler developed in [28] for the approximate reconstruction of bandlimited functions is adopted, which converts any continuous signal 𝑓(𝑡) into a square wave function by integrating it against an averaging function 𝑢_𝑘,𝑧(𝑡). The integrated result is compared to a positive threshold and a negative threshold, such that when either one of the two thresholds is reached, a pulse is generated at time 𝑡_𝑘=𝑧. The value of the integrator is then reset and the process repeats. In this paper, the averaging function 𝑢_𝑘,𝑧(𝑡) is chosen to be the exponential 𝑒^{𝛼(𝑡−𝑧)}𝑋_{[𝑡𝑘,𝑧]}, where 𝑋_𝐼 is the characteristic function of 𝐼 and 𝛼>0 is a constant that models the leakage of the integrator, as due to practical implementations [28]. Then, the LIF sampler firing condition that generates the square wave (or sequence of pulses) is ±𝜃=^𝑡𝑘+1_𝑡𝑘𝑓(𝑡)𝑒^{𝛼(𝑡−𝑡𝑘+1)}𝑑𝑡≜𝑓,𝑢_𝑘, (26) where 𝑡_𝑘 is the time instant of the sample, 𝑡_𝑘+1 is the next time instant of the sample, 𝑢_𝑘 is the averaging function, 𝜃 is the threshold value, and 𝛼 is the leakage of the integrator. The output of the LIF sampler can be expressed in terms of the time instants at which the integral reaches the threshold, {𝑡₀,…,𝑡_𝑛}, and by the samples 𝑞₁,…,𝑞_𝑛, defined as 𝑞_𝑗≜^𝑡𝑗−1_𝑡𝑗𝑓(𝑥)𝑒^{𝛼(𝑥−𝑡𝑗)}𝑑𝑥,for1≤𝑗≤𝑛. (27)

From (26), it follows that |𝑞_𝑗|=𝜃. By adjusting the parameters of the LIF sampler, it is possible to convert the RBF signal in (20) to a square wave comprised of pulses with the same width, height, and centers as the RBFs in (20). Thus, by using the RBF spike model in (26) with suitable widths and heights, it is possible to induce the first neuron to fire shortly before the end of each pulse of the square-pulse function (also considering the refractory period of the neuron). For simplicity, in this example, it is assumed that the heights 𝑤_𝑘 and widths 𝛽_𝑘 are known positive constants of equal magnitudes for all 𝑘. Then, the centers of the RBF comprise the set of adjustable parameters, 𝑃={𝑐_𝑘∣𝑘=1,…,𝑁}, to be optimized. The same approach can be easily extended to a case where all of the RBF parameters are adapted.

It can be shown using the LIF SNN model in Section 2 that, under the aforementioned assumptions, the firing times of the first neuron satisfy the constraints, 𝑡_1,𝑘≤𝑐_𝑘+𝛽2,𝑡_1,𝑘+Δ^abs≥𝑐_𝑘+𝛽2, (28) where Δ^abs is an absolute refractory time of the neuron. Then, the firing times of the first neuron can be written as a function of the RBF centers, 𝑐_𝑘. By design, the first neuron fires after the same time interval, relative to each square pulse (Figure 4), due to the chosen RBF heights and widths. Then, the firing times of the first neuron can be written as, 𝑡_1,𝑘=𝑐_𝑘−𝛽2+𝑇, (29) where 𝑡_1,𝑘 denotes the firing times of the first neuron, 𝑐_𝑘 are the centers of the RBF, 𝛽 is the constant width of the RBF, and 𝑇 is given by (21).

The RBF spike model is demonstrated in Figure 5, where an example of RBF output obtained from (20) is plotted and, after being fed to the LIF sampler, produces a corresponding square wave which can be used as controlled pulse. It can be seen that the centers of the RBFs precisely determine the times at which a pulse occurs in the square wave and that the square wave maintains the desired constant width and magnitude for all 𝑘. The next subsection illustrates how, by adapting the continuous RBF spike model in (20), it is possible to minimize a desired error function and train the synaptic weight of the SNN without directly manipulating it.

(a)

(b)

Figure 5 

Deterministic spike model signals.

4.2. Indirect Training Equations

As explained in Section 2.2, it is assumed that the synaptic weight 𝑤₂₁ can only be modified by controlling the activity of the SNN input neuron(s), with 𝑖=1, and that it obeys the nearest-spike STDP model in (8). Over time, the synaptic weight can change repeatedly, and therefore its final value can be written as the sum of all incremental changes that have occurred over the time interval [𝑡_𝜄,𝑡_𝜄+1], 𝑤₂₁(𝑡)=^𝑀_𝑖=1𝑤₀1+Δ𝑤_𝑖, (30) where, for the two-neuron SNN in Figure 3, Δ𝑤₁,…,Δ𝑤_𝑀 are due to 𝑀 pairs of pre- and postsynaptic spikes that occur at any time 𝑡∈[𝑡_𝜄,𝑡_𝜄+1]. Every weight increment Δ𝑤_𝑖 is induced via STDP and, thus, obeys equation (8). Since the firing times in (8) depend on the centers of the RBF input through (21)–(29), it follows that every weight increment Δ𝑤_𝑖 is a function of the RBF centers. In this example, the constraints (28) can be written as, 𝛽2<𝑐₁,𝑐_𝑁+𝛽2<𝑡_𝑓,𝑐_𝑘+𝛽<𝑐_𝑘+1,for𝑘=1,…,𝑁, (31) where, from Figure 5, 𝑁=5.

A training-error function is defined in terms of the desired synaptic weight 𝑤^∗ and the actual value of the synaptic weight 𝑤₂₁(𝑡). While different forms of error function may be used, the chosen form determines the complexity of the derivation of the training gradients. It was found that the most convenient form of training-error function can be derived from the ratio of the actual weight over the desired weight, 𝑒(𝑡)=𝑤₂₁(𝑡)𝑤^∗, (32) where 𝑤₂₁(𝑡) is the weight value at time 𝑡∈[𝑡_𝜄,𝑡_𝜄+1], obtained from all previous spikes. As explained in Section 3, 𝑤^∗ can be viewed as the weight that leads to desired output 𝑦 for a given input 𝜉. It can be seen that when 𝑤₂₁ is equal to 𝑤^∗, 𝑒 is equal to one. Plugging (30) into (36), the weight ratio can be rewritten as, 𝑒(𝑡)=1𝑤^∗𝑀_𝑖=1𝑤₀1+Δ𝑤_𝑖, (33) and the error between 𝑤₂₁ and 𝑤^∗ can be minimized by minimizing the natural logarithm of the ratio (33): 𝐸(𝑡)≜ln[𝑒(𝑡)]=ln𝑤₀+ln1+Δ𝑤₁+⋯+ln1+Δ𝑤_𝑀−ln𝑤^∗. (34)

As is typical of all optimization problems, minimizing a quadratic form presents several advantages [29]. Therefore, at any time 𝑡∈[𝑡_𝜄,𝑡_𝜄+1], the indirect training algorithm seeks to minimize the quadratic training-error function: 𝐽(𝑡)≜𝐸(𝑡)²={ln[𝑒(𝑡)]}²=ln𝑤₀+ln1+Δ𝑤₁+⋯+ln1+Δ𝑤_𝑀−ln𝑤^∗2. (35) Then, indirect training can be formulated as an unconstrained optimization problem in which 𝐽 is to be minimized with respect to the RBF centers, or 𝑃={𝑐_𝑘∶𝑐_𝑘∈[𝑡_𝜄,𝑡_𝜄+1]}. Any gradient-based numerical optimization algorithm can be utilized for this purpose. The analytical form of the gradient 𝜕𝐽/𝜕𝑐_𝑘 depends on the form of the spike patterns. The following subsection derives this gradient analytically for one example of spike patterns and demonstrates that, using this gradient, the synaptic weight can be trained indirectly to meet the desired value 𝑤^∗ exactly. The same results were derived and demonstrated numerically for all other possible spike patterns, but they are omitted here for brevity. Another representation of error 𝑒′(𝑡) is defined below to make the convergence of error more understandable in figures, 𝑒^′(𝑡)=𝑤₂₁−𝑤^∗. (36) where, |⋅| denotes the absolute value.

4.3. Derivation of Gradient Equations for Indirect Training

Consider the case in which 𝑀=5, 𝑤^∗=1.72, and 2𝑎_𝐸>𝑉_th−𝑉₀>𝑎_𝐸, which results in a two-node SNN (Figure 3) in which neuron 𝑖=1 must fire at least two times in order for neuron 𝑖=2 to fire. An example of such a spike pattern is shown in Figure 6, where 𝑁₁ and 𝑁₂ denote spike trains of neuron 𝑖=1 and 𝑖=2, respectively. In this case, the first neuron fires five times, and the second neuron fires two times. The firing time of the first neuron is controlled by the RBF spike model described in Section 4.1. In Figure 6, 𝑡_𝑖,𝑘 denotes the 𝑘th firing time of the 𝑖th neuron with 𝑘∈ℐ_𝑖, where ℐ_𝑖={𝑘=1,…,5} is an index set for the firing times and 𝑖 is the index labeling the neurons. In this example, the second neuron fires after 𝑡_1,3 and 𝑡_1,5, because 𝑡_1,2, 𝑡_1,3 and 𝑡_1,4, 𝑡_1,5 are close enough to cause the membrane potential of the second neuron, 𝑣₂, to exceed the threshold.

(a)

(b)

Figure 6 

Optimized input spike train and indirect weight changes brought about by STDP.

Initially, the synaptic weight is equal to 1.7. The weight change is discontinuous due to the discrete property of spikes. As shown in Figure 6, the weight increases three times at 𝑡_1,3+𝜏_𝑑, decreases at 𝑡_1,4, and increases again at 𝑡_1,5+𝜏_𝑑. For this example, the synaptic weight changes in five increments: Δ𝑤₁=𝐴₊𝑡exp^1,1𝜏₊exp−𝑡_1,3𝜏₊exp−𝜏_𝑑𝜏₊Δ𝑤₂=𝐴₊𝑡exp^1,2𝜏₊exp−𝑡_1,3𝜏₊exp−𝜏_𝑑𝜏₊Δ𝑤₃=𝐴₊exp−𝜏_𝑑𝜏₊Δ𝑤₄=−𝐴₋𝑡exp^1,3𝜏₋exp−𝑡_1,4𝜏₋𝜏exp^𝑑𝜏₋Δ𝑤₅=𝐴₊exp−𝜏_𝑑𝜏₊. (37) When the equations above are substituted in (35), the training-error function can be written as 𝐽=ln𝑤₀+ln1+𝐴₊𝑡exp^1,1𝜏₊exp−𝑡_1,3𝜏₊exp−𝜏_𝑑𝜏₊+ln1+𝐴₊𝑡exp^1,2𝜏₊exp−𝑡_1,3𝜏₊exp−𝜏_𝑑𝜏₊+ln1+𝐴₊exp−𝜏_𝑑𝜏₊+ln1−𝐴₋𝑡exp^1,3𝜏₋exp−𝑡_1,4𝜏₋𝜏exp^𝑑𝜏₋+ln1+𝐴₊exp−𝜏_𝑑𝜏₊−ln𝑤^∗2. (38) Then, the gradients of 𝐽 with respect to the RBF centers are given by 𝜕𝐽𝜕𝑐₁=𝜕𝐸²𝜕𝑐₁=𝜕𝐸²𝜕𝑡_1,1=2𝐸Δ𝑤₁𝜏₊1+Δ𝑤₁𝜕𝐽𝜕𝑐₂=𝜕𝐸²𝜕𝑐₂=𝜕𝐸²𝜕𝑡_1,2=2𝐸Δ𝑤₂𝜏₊1+Δ𝑤₂𝜕𝐽𝜕𝑐₃=𝜕𝐸²𝜕𝑐₃=𝜕𝐸²𝜕𝑡_1,3=2𝐸−Δ𝑤₁𝜏₊1+Δ𝑤₁+−Δ𝑤₂𝜏₊1+Δ𝑤₂+Δ𝑤₄𝜏₋1+Δ𝑤₄𝜕𝐽𝜕𝑐₄=𝜕𝐸²𝜕𝑐₄=𝜕𝐸²𝜕𝑡_1,4=2𝐸−Δ𝑤₄𝜏₋1+Δ𝑤₄𝜕𝐽𝜕𝑐₅=𝜕𝐸²𝜕𝑐₅=𝜕𝐸²𝜕𝑡_1,5=0. (39) Using the above gradients, the optimal values of 𝑐₁,…,𝑐₅ can be obtained by minimizing 𝐽 using a gradient-based numerical algorithm such as Newton’s method.

An example of indirect learning algorithm implementation is shown in Figure 6, where the RBF-adjustable parameters (which, in this case, coincide with the firing times of neuron 𝑖=1) are optimized to induce a change in the synaptic weight 𝑤₂₁ from an initial value of 1.704, to a desired value 𝑤^∗=1.72. Another example, for which the gradient equations are omitted for brevity, is shown in Figures 7–9. In this case, the desired weight value 𝑤^∗=2.8 is far from the initial weight 𝑤₂₁(𝑡_𝜄)=3.5, and thus 𝑁=43 spikes are required to train the SNN. The optimal RBF input and corresponding controlled pulse used to stimulate neuron 𝑖=1 are shown in Figures 7 and 8, respectively. The time history of the weight 𝑤₂₁ is plotted in Figure 9(a), along with the corresponding deviation from 𝑤^∗ plotted in Figure 9(b).

Figure 7 

Optimized RBF spike model for the SNN in Figure 3.

Figure 8 

Square wave obtained by the LIF sampler, for the RBF spike model in Figure 7.

(a)

(b)

Figure 9 

Indirect weight changes brought about by STDP and corresponding deviation from 𝑤∗=2.8, for the SNN in Figure 3 stimulated using the input spike train in Figure 7.

5. Generalized Form of Gradient Equations for Indirect Training

A more general form of the gradient equations can be found by rewriting the response of the membrane potential for neuron 𝑖=2 in (24), as follows, 𝑣(𝑡)=𝑉₀+^𝑏_𝑘=𝑎𝑎_𝐸exp−𝑡−𝑡_1,𝑘−𝜏_𝑑𝜏_𝑚𝐻𝑡−𝑡_1,𝑘−𝜏_𝑑, (40) where 𝑡_1,𝑘 denotes the 𝑘th firing time of neuron 𝑖=1. Then, the gradients of the training objective function (35) with respect to the centers of the RBF spike model for 𝑀=5 are given by the equations in the Appendix, which are obtained in terms of the two functions, 𝑔₊𝑡_1,𝑖,𝑡_1,𝑗=𝐴₊𝑡exp^1,𝑖𝜏₊exp−𝑡_1,𝑗𝜏₊exp−𝜏_𝑑𝜏₊,𝑔₋𝑡_1,𝑖,𝑡_1,𝑗=−𝐴₋𝑡exp^1,𝑖𝜏₋exp−𝑡_1,𝑗𝜏₋𝜏exp^𝑑𝜏₋, (41) where 𝑡_1,𝑖, 𝑡_1,𝑗 are two distinct firing times of neuron 𝑖=1.

The methodology presented in Section 4 can also be extended to larger SNNs, although in this case it may be more convenient to compute the gradients of the objective function numerically. As an example, consider the three-neuron SNN in Figure 10, modeled by the approach in Section 2 and with two synaptic weights 𝑤₂₁ and 𝑤₂₃ that each obey the STDP mechanism in (8). Suppose that the desired values of the synaptic weights are 𝑤^∗₂₁=4.1 and 𝑤^∗₂₃=2.8. Using the indirect learning method presented in this paper, and the gradient provided in the Appendix, a training-error function formulated in terms of the deviations of 𝑤₂₁(𝑡) and 𝑤₂₃(𝑡) from 𝑤^∗₂₁ and 𝑤^∗₂₃, respectively, can be minimized with respect to the parameters (centers) 𝑃 of the RBF spike model (20).

Figure 10 

Model of three-node LIF spiking neural network.

Once the optimal RBF spike model, plotted in Figure 11, is fed to the LIF sampler, the controlled pulse plotted in Figure 12 is obtained and implemented via 𝐼_inj. The controlled pulse is thus used to stimulate neuron 𝑖=1 at precise instants in time that corresponds to centers of the optimal RBF spike model. By this approach, 𝑤₂₁(𝑡) can be made to converge to the desired value 𝑤^∗₂₁. The same method is implemented to stimulate neuron 𝑖=2 to make the value of 𝑤₂₃(𝑡) converge to 𝑤^∗₂₃. The time histories of the weight values 𝑤₂₁(𝑡) and 𝑤₂₃(𝑡), obtained by the indirect training algorithm are plotted in Figure 13(a). As is also shown by the corresponding training errors, plotted in Figure 13(b), the SNN weights over time converge to the desired values, that is 𝑤^∗₂₁=4.1 and 𝑤^∗₂₃=2.8. These results demonstrate that, even for larger SNNs, the indirect training method presented in this paper is capable of modifying synaptic weights until they meet their desired values, without direct manipulation. Since this indirect training algorithm only relies on modulating the activity of the input neurons via controlled input spike trains, it also has the potential of being realizable in vitro and in silico to train biological neuronal networks and CMOS/memristor nanoscale chips, respectively.

Figure 11 

Optimized RBF spike model for the SNN in Figure 10.

Figure 12 

Square wave obtained from the LIF sampler, using the optimized RBF spike model in Figure 11.

(a)

(b)

Figure 13 

Indirect weight changes brought about by STDP and corresponding error functions, for the SNN in Figure 10 stimulated using the input spike train in Figure 12.

6. Summary and Conclusion

Recently, several algorithms have been proposed for training spiking neural networks through biologically plausible learning mechanisms, such as spike-timing-dependent synaptic plasticity. These algorithms, however, rely on being able to modify the synaptic weights directly or STDP. In other words, they minimize a desired objective function with respect to weight increments that are assumed to be controllable and, thus, are decided by a weight update rule, and then are implemented directly by the training algorithm. In several potential applications of spiking neural networks, synaptic weights cannot be manipulated directly but change over time by virtue of pre- and postsynaptic neural activity. In these applications, the activity of selected input neurons can be controlled via programming voltages or pulses of blue light that induce precise spiking of the input neurons, at precise moments in time.

This paper presents an indirect learning method that induces changes in the synaptic weights by modulating spike-timing-dependent plasticity using controlled input spike trains, in lieu of the weight increments. The key difficulty to be overcome is that indirect learning seeks to adapt a pulse signal, such as a square wave or Rademacher function, in place of continuous-valued weights. Pulse signals that can be delivered to stimulate input neurons and cause them to spike, such as blue light patterns or programming voltages, are represented by piece-wise continuous, multi-valued (or many-to-one), and nondifferentiable functions that are not well suited to numerical optimization. Furthermore, stimulation patterns typically are generated by spike models that are stochastic, such as the Poisson spike model. Therefore, even when the spike model is optimized, it does not allow for precise timing of pre- and post-synaptic firings, and as a result, may induce undesirable changes in the synaptic weights. This paper presents a deterministic and adaptive spike model derived from radial basis functions and a leaky integrate-and-fire sampler. This spike model can be easily optimized to determine the sequence of spike timings that minimizes a desired objective function and, then, used to stimulate input neurons. The results demonstrate that this methodology is capable of inducing the desired synaptic plasticity in the network and modify synaptic weights to meet their desired values only by virtue of controlled input spike trains that are realizable both in biological neuronal networks and in CMOS/memristor nanoscale chips.

Appendix

The gradients of the training objective function can be derived as follows. The membrane potential of neuron 2 is denoted by 𝑣_𝑖.𝑗(𝑡), where the membrane potential can only be affected by presynaptic spikes that occur within the time interval [𝑡_1.𝑖,𝑡_1,𝑗].𝜕𝐸²𝜕𝑐₁=𝑔2𝐸⁺𝑡_1,1,𝑡_1,3𝜏₊1+𝑔₊𝑡_1,1,𝑡_1,3if𝑣𝑡_1,3+𝜏_𝑑1,3>𝑉_th𝑔2𝐸⁺𝑡_1,1,𝑡_1,2𝜏₊1+𝑔₊𝑡_1,1,𝑡_1,2if𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th𝑔2𝐸⁺𝑡_1,1,𝑡_1,4𝜏₊1+𝑔₊𝑡_1,1,𝑡_1,4if𝑣𝑡_1,4+𝜏_𝑑1,4>𝑉_th𝑔2𝐸⁺𝑡_1,1,𝑡_1,5𝜏₊1+𝑔₊𝑡_1,1,𝑡_1,5if𝑣𝑡_1,5+𝜏_𝑑1,5>𝑉_th.𝜕𝐸²𝜕𝑐₂=𝑔2𝐸⁺𝑡_1,2,𝑡_1,3𝜏₊1+𝑔₊𝑡_1,2,𝑡_1,3if𝑣𝑡_1,3+𝜏_𝑑1,3>𝑉_th2𝐸−𝑔₊𝑡_1,1,𝑡_1,2𝜏₊1+𝑔₊𝑡_1,1,𝑡_1,2if𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑3,5>𝑉_th𝑐₃>𝑐₅+𝑐₂+2𝜏_𝑑22𝐸−𝑔₊𝑡_1,1,𝑡_1,2𝜏₊1+𝑔₊𝑡_1,1,𝑡_1,2+𝑔₋𝑡_1,2,𝑡_1,3𝜏₋1+𝑔₋𝑡_1,2,𝑡_1,3if𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,4+𝜏_𝑑3,4>𝑉_th,𝑐₃<𝑐₄+𝑐₂+2𝜏_𝑑2or𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑3,5>𝑉_th,𝑐₃<𝑐₅+𝑐₂+2𝜏_𝑑22𝐸−𝑔₊𝑡_1,1,𝑡_1,2𝜏₊1+𝑔₊𝑡_1,1,𝑡_1,2+𝑔₋𝑡_1,2,𝑡_1,3𝜏₋1+𝑔₋𝑡_1,2,𝑡_1,3+𝑔₋𝑡_1,2,𝑡_1,4𝜏₋1+𝑔₋𝑡_1,2,𝑡_1,4+𝑔₋𝑡_1,2,𝑡_1,5𝜏₋1+𝑔₋𝑡_1,2,𝑡_1,5if𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑3,5<𝑉_th𝑔2𝐸⁺𝑡_1,2,𝑡_1,4𝜏₊1+𝑔₊𝑡_1,2,𝑡_1,4if𝑣𝑡_1,4+𝜏_𝑑1,4>𝑉_th.𝜕𝐸²𝜕𝑐₃=2𝐸−𝑔₊𝑡_1,1,𝑡_1,3𝜏₊1+𝑔₊𝑡_1,1,𝑡_1,3+−𝑔₊𝑡_1,2,𝑡_1,3𝜏₊1+𝑔₊𝑡_1,2,𝑡_1,3+𝑔₋𝑡_1,3,𝑡_1,4𝜏₋1+𝑔₋𝑡_1,3,𝑡_1,4if𝑣𝑡_1,3+𝜏_𝑑1,3>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑4,5>𝑉_th,𝑐₃>2𝑐₄−𝑐₅−2𝜏_𝑑2𝐸−𝑔₊𝑡_1,1,𝑡_1,3𝜏₊1+𝑔₊𝑡_1,1,𝑡_1,3+−𝑔₊𝑡_1,2,𝑡_1,3𝜏₊1+𝑔₊𝑡_1,2,𝑡_1,3if𝑣𝑡_1,3+𝜏_𝑑1,3>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑4,5>𝑉_th,𝑐₃<2𝑐₄−𝑐₅−2𝜏_𝑑𝑔2𝐸⁺𝑡_1,3,𝑡_1,5𝜏₊1+𝑔₊𝑡_1,3,𝑡_1,5if𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑3,5>𝑉_th,𝑐₃>𝑐₅+𝑐₂+2𝜏_𝑑2or𝑣𝑡_1,5+𝜏_𝑑1,5>𝑉_th2𝐸−𝑔₋𝑡_1,2,𝑡_1,3𝜏₋1+𝑔₋𝑡_1,2,𝑡_1,3if𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑3,5>𝑉_th,𝑐₃<𝑐₅+𝑐₂+2𝜏_𝑑2or𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,4+𝜏_𝑑3,4>𝑉_th,𝑐₃<𝑐₄+𝑐₂+2𝜏_𝑑2or𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑3,5<𝑉_th𝑔2𝐸⁺𝑡_1,3,𝑡_1,4𝜏₊1+𝑔₊𝑡_1,3,𝑡_1,4if𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,4+𝜏_𝑑3,4>𝑉_th,𝑐₃>𝑐₄+𝑐₂+2𝜏_𝑑2or𝑣𝑡_1,4+𝜏_𝑑1,4>𝑉_th2𝐸−𝑔₊𝑡_1,1,𝑡_1,3𝜏₊1+𝑔₊𝑡_1,1,𝑡_1,3+−𝑔₊𝑡_1,2,𝑡_1,3𝜏₊1+𝑔₊𝑡_1,2,𝑡_1,3+𝑔₋𝑡_1,3,𝑡_1,4𝜏₋1+𝑔₋𝑡_1,3,𝑡_1,4+𝑔₋𝑡_1,3,𝑡_1,5𝜏₋1+𝑔₋𝑡_1,3,𝑡_1,5if𝑣𝑡_1,3+𝜏_𝑑1,3>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑4,5<𝑉_th.𝜕𝐸²𝜕𝑐₄=2𝐸−𝑔₋𝑡_1,3,𝑡_1,4𝜏₋1+𝑔₋𝑡_1,3,𝑡_1,4if𝑣𝑡_1,3+𝜏_𝑑1,3>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑4,5>𝑉_th,𝑐₄<𝑐₃+𝑐₅+2𝜏_𝑑2or𝑣𝑡_1,3+𝜏_𝑑1,3>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑4,5<𝑉_th𝑔2𝐸⁺𝑡_1,4,𝑡_1,5𝜏₊1+𝑔₊𝑡_1,4,𝑡_1,5if𝑣𝑡_1,3+𝜏_𝑑1,3>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑4,5>𝑉_th,𝑐₄>𝑐₃+𝑐₅+2𝜏_𝑑2or𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑3,5>𝑉_th,𝑐₃>𝑐₂+𝑐₅+2𝜏_𝑑2or𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑3,5>𝑉_th,𝑐₃<𝑐₂+𝑐₅+2𝜏_𝑑22𝐸−𝑔₊𝑡_1,3,𝑡_1,4𝜏₊1+𝑔₊𝑡_1,3,𝑡_1,4+𝑔₋𝑡_1,4,𝑡_1,5𝜏₋1+𝑔₋𝑡_1,4,𝑡_1,5if𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,4+𝜏_𝑑3,4>𝑉_th,𝑐₄<2𝑐₃−𝑐₂−2𝜏_𝑑𝑔2𝐸⁻𝑡_1,4,𝑡_1,5𝜏₋1+𝑔₋𝑡_1,4,𝑡_1,5if𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,4+𝜏_𝑑3,4>𝑉_th,𝑐₄>2𝑐₃−𝑐₂−2𝜏_𝑑2𝐸−𝑔₋𝑡_1,2,𝑡_1,4𝜏₋1+𝑔₋𝑡_1,2,𝑡_1,4if𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑3,5<𝑉_th2𝐸−𝑔₊𝑡_1,1,𝑡_1,4𝜏₊1+𝑔₊𝑡_1,1,𝑡_1,4+−𝑔₊𝑡_1,2,𝑡_1,4𝜏₊1+𝑔₊𝑡_1,2,𝑡_1,4+−𝑔₊𝑡_1,3,𝑡_1,4𝜏₊1+𝑔₊𝑡_1,3,𝑡_1,4+𝑔₋𝑡_1,4,𝑡_1,5𝜏₋1+𝑔₋𝑡_1,4,𝑡_1,5if𝑣𝑡_1,4+𝜏_𝑑,1,4>𝑉_th𝑔2𝐸⁺𝑡_1,4,𝑡_1,5𝜏₊1+𝑔₊𝑡_1,4,𝑡_1,5if𝑣𝑡_1,5+𝜏_𝑑,1,5>𝑉_th.𝜕𝐸²𝜕𝑐₅=0if𝑣𝑡_1,3+𝜏_𝑑1,3>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑4,5>𝑉_th,𝑐₅>2𝑐₄−𝑐₃−2𝜏_𝑑2𝐸−𝑔₊𝑡_1,4,𝑡_1,5𝜏₊1+𝑔₊𝑡_1,4,𝑡_1,5if𝑣𝑡_1,3+𝜏_𝑑1,3>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑4,5>𝑉_th,𝑐₅<2𝑐₄−𝑐₃−2𝜏_𝑑2𝐸−𝑔₊𝑡_1,3,𝑡_1,5𝜏₊1+𝑔₊𝑡_1,3,𝑡_1,5+−𝑔₊𝑡_1,4,𝑡_1,5𝜏₊1+𝑔₊𝑡_1,4,𝑡_1,5if𝑣(𝑡_1,2+𝜏_𝑑)_1,2>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑3,5>𝑉_th,𝑐₅<2𝑐₃−𝑐₂−2𝜏_𝑑2𝐸−𝑔₊𝑡_1,4,𝑡_1,5𝜏₊1+𝑔₊𝑡_1,4,𝑡_1,5if𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑3,5>𝑉_th,𝑐₅>2𝑐₃−𝑐₂−2𝜏_𝑑2𝐸−𝑔₋t_1,4,𝑡_1,5𝜏₋1+𝑔₋𝑡_1,4,𝑡_1,5if𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,4+𝜏_𝑑3,4>𝑉_th2𝐸−𝑔₋𝑡_1,2,𝑡_1,5𝜏₋1+𝑔₋𝑡_1,2,𝑡_1,5if𝑣𝑡_1,2+𝜏_𝑑1,2>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑3,5<𝑉_th2𝐸−𝑔₋𝑡_1,3,𝑡_1,5𝜏₋1+𝑔₋𝑡_1,3,𝑡_1,5if𝑣𝑡_1,3+𝜏_𝑑1,3>𝑉_th,𝑣𝑡_1,5+𝜏_𝑑4,5<𝑉_th2𝐸−𝑔₋𝑡_1,4,𝑡_1,5𝜏₋1+𝑔₋𝑡_1,4,𝑡_1,5if𝑣𝑡_1,4+𝜏_𝑑1,4>𝑉_th2𝐸−𝑔₊𝑡_1,1,𝑡_1,5𝜏₊1+𝑔₊𝑡_1,1,𝑡_1,5+−𝑔₊𝑡_1,2,𝑡_1,5𝜏₊1+𝑔₊𝑡_1,2,𝑡_1,5+−𝑔₊𝑡_1,3,𝑡_1,5𝜏₊1+𝑔₊𝑡_1,3,𝑡_1,5+−𝑔₊𝑡_1,4,𝑡_1,5𝜏₊1+𝑔₊𝑡_1,4,𝑡_1,5if𝑣𝑡_1,5+𝜏_𝑑1,5>𝑉_th. (A.1)

Acknowledgment

This work was supported by the National Science Foundation, under ECCS Grant 0925407.