Pre-Synaptic Pool Modification (PSPM): A supervised learning procedure for recurrent spiking neural networks

Learning synaptic weights of spiking neural network (SNN) models that can reproduce target spike trains from provided neural firing data is a central problem in computational neuroscience and spike-based computing. The discovery of the optimal weight values can be posed as a supervised learning task wherein the weights of the model network are chosen to maximize the similarity between the target spike trains and the model outputs. It is still largely unknown whether optimizing spike train similarity of highly recurrent SNNs produces weight matrices similar to those of the ground truth model. To this end, we propose flexible heuristic supervised learning rules, termed Pre-Synaptic Pool Modification (PSPM), that rely on stochastic weight updates in order to produce spikes within a short window of the desired times and eliminate spikes outside of this window. PSPM improves spike train similarity for all-to-all SNNs and makes no assumption about the post-synaptic potential of the neurons or the structure of the network since no gradients are required. We test whether optimizing for spike train similarity entails the discovery of accurate weights and explore the relative contributions of local and homeostatic weight updates. Although PSPM improves similarity between spike trains, the learned weights often differ from the weights of the ground truth model, implying that connectome inference from spike data may require additional constraints on connectivity statistics. We also find that spike train similarity is sensitive to local updates, but other measures of network activity such as avalanche distributions, can be learned through synaptic homeostasis.


Introduction:
Optimizing the synaptic weights of artificial spiking neural networks (SNNs) to reproduce a desired input/output function is a particularly complex supervised learning problem.SNN neuron models, like the LIF neuron, usually have membrane potentials which obey a differential equation while the potential is less than a threshold, but reset upon surpassing the threshold value.Each reset is recognized as a spike and is stored as a 1 in the binary output spike train for the corresponding neuron and time step.Synapses between neurons allow the firing of neuron j to influence the membrane potential of neuron i, either in an excitatory or an inhibitory manner.To optimize the SNN, one minimizes a numerical distance between the output spike train of the network and the reference spike trains by changing the synapses of the network.Given this problem's large solution space and the non-analyticity of it's loss function, gradient-based procedures are inapplicable and approximations are required.A given spike may be caused by the strength of the synapses from other spiking neurons or due to external stimulating current.This potential over-determination of spike events complicates any attempt to reproduce a desired input/output relation.In addition, the space of possible fully recurrent neural networks with weight matrix values ranging over a discrete set (with, say d elements) expands exponentially (like d N 2 ) with the number of neurons N in the network, rendering the solution to these input/output transformations in a network with many neurons a significant challenge.
While past work (Bohte et al., 2002) has successfully developed an analogue for spiking neural networks of the traditional backpropagation algorithm, it requires a layered, feed-forward structure.Layered network models are inspired by the large-scale structure of the cortex, but fail to account for the highly recurrent nature of individual layers in a biological neural network.Though this layered structure is useful for engineering applications, it lacks the feedback central to biological neural networks.A number of efforts have been made to develop algorithms for optimizing recurrent SNNs (Gutig and Sompolinsky, 2006), (Pavlidis et al., 2005), (Russell et al., 2010), (Jin et al., 2007), (Booij and tat Nguyen, 2005), (Schrauwen and Campenhout, 2004), (Schrauwen and Campenhout, 2006), but these efforts have been limited in the scope of the problems they attempt to solve.In particular, previous algorithms have imposed a limit of one spike per neuron or were successful only with networks with less than 50 neurons.Each limit relaxes the problem of exact spike time matching on an arbitrary set of desired spike trains.One other study (Gardner and Grüning, 2016) developed a more general-purpose learning rule, but it was primarily intended for machine learning purposes, whereas here we aim primarily to produce model spiking neural networks that capture the inputoutput relationship observed in data from other networks.Motivated by the stochastic nature of biological neural computations, our relaxation, to which we refer as Pre-Synaptic Pool Modification, removes the requirement of exactly matching spike timings.Instead we seek to generate artificial SNNs whose spiking behavior matches the statistical features of some desired set of spike trains, in particular a van-Rossum like distance (van Rossum, 2001), agreement between inter-spike-interval distributions, and critical behavior (Plenz and Niebur, 2014).This relaxation holds particular and immediate usefulness in computational neuroscience by allowing the creation of generative models for neural data.Pre-Synaptic Pool Modification (PSPM) may also be extended to applications with artificial neural networks.
Our PSPM algorithm is based on the concept of synaptic scaling (Turrigiano, 2017)- (Keck et al., 2017).In synaptic scaling, neurons respond to excess or inadequate network-level activity by modifying the strength of their synapses.In response to strong and weak external inputs, this modification prevents, respectively, oversaturation and extinction of spiking activity.Synaptic scaling is crucial in recurrent networks, which are otherwise at risk of runaway activity resulting from feedback loops within the network.In addition to this network-wide synaptic scaling, our algorithm also changes individual synaptic weights in the network, modulating the activity neuron by neuron to eliminate unwanted spikes and induce desired spikes.
We implemented the proposed PSPM learning rules and used them to optimize the weight matrices of fully recurrent LIF neural networks.Application our proposed technique improved the similarity between the optimized and reference spike trains on a variety of measures including integrated differences of convolved activity signals, similarity of inter-spike interval distributions and a power-law distribution of neuronal avalanches.

Methods:
Leaky Integrate-and-Fire (LIF) Model Our neural networks consisted of Leaky Integrate-and-Fire (LIF) neurons.For all neurons i in the network, the membrane potential V i (t) had a time dependence given by the following equations: Where τ is membrane time constant, V th is the voltage threshold, R m is the membrane resistance, and I i (t) is the external input current into neuron i at time t.Further, W ij is the weight of the network connection from neuron j to neuron i, and s j (t) is the binary value denoting whether or not neuron j spiked at time t.V th represents the threshold membrane potential above which a neuron will spike.In our simulations, we imposed the condition that 20% of the neurons in our network are inhibitory with W ij < 0 for all inhibitory neurons j.Network connectivity was all-to-all but we allowed W ij = 0, corresponding to an absence of synaptic strength between neuron j and neuron i.For our model, τ = 30 ms, R m = 100 M Ω, and V th = 30 mV , which are biologically plausible parameters for a real neuron.The differential equation is evaluated discretely using Euler's method with a step size of 3 ms.

Initialized Network Parameters
In our simulations, we used one LIF network, termed the reference network, to represent the goal input/output map, from input stimulating currents to output spike trains.A distinct network called the naive network was also generated and then optimized with our learning rules to produce output spike trains similar to those of the reference.All networks consisted of N = 400 LIF neurons with 20% of the neurons j randomly assigned as inhibitory (W ij < 0 for all inhibitory neurons j) and the remaining 80% set as excitatory (W ij > 0 for all excitatory neurons j).
Initial weight matrices of the reference and naive networks had magnitudes that were drawn from 4 distinct distributions: (1) uniform: reference and naive network weight matrix magnitudes were both drawn from a uniform distribution on [0 V, 5×10 −3 V ], (2) gaussian: reference and naive weight values were both drawn from a gaussian distribution with (µ = 4×10 −4 V , σ = 4×10 −4 V ), (3) sparse: reference and naive weight matrix values first drawn from a uniform distribution on [0 V, 5 × 10 −3 V ] with 50% of the synapses in the naive network subsequently set to 0 V , (4) naive-half-max : reference weight magnitudes drawn from a uniform distribution on [0 V, 5 × 10 −3 V ] but naive network weight magnitudes drawn from a uniform distribution on [0 V, 2.5 × 10 −3 V ].As described above, the 20% of neurons j designated as inhibitory had weight matrix values set to negative (W ij < 0 for all inhibitory neurons j).The excitatory neurons j had positive weight values (W ij > 0 for all excitatory neurons j).The same neurons j that were inhibitory in the reference network were inhibitory in the naive network.The variety of initial weight distributions was chosen to demonstrate the versatility of the algorithm for uniform and gaussian networks.The sparse and naive-half-max initial configurations demonstrate that the algorithm can improve the performance of a naive network with weights drawn from an initially distinct distribution from that of the reference.
The external input current values at each time step were drawn from a Gaussian distribution with mean µ = 2.5 × 10 −10 A and standard deviation σ = 1 × 10 −10 A. The values of the input currents were tuned to be strong enough to produce a network-wide mean spike frequency 0 < f ≤ 1 Hz when all weights were set to zero, but could produce mean spike rates greater than 10 Hz when the weights were nonzero.This restriction guaranteed that the external input currents remained smaller than the currents generated by synaptic connections within the network, and thus the network weights significantly affected the output spike trains.The distribution described above met these requirements.

Learning Paradigm
We conducted 30 trials for each of the 4 sets of network parameters described above (Initial Network Parameters), creating a reference network and a naive network by drawing weights from the specified distributions.The reference and naive networks are simulated with the same input currents to generate reference outputs R and naive outputs N respectively.The naive network's weights are then changed according to the PSPM rules (PSPM Learning Rules).This modified network is simulated again with the same inputs to produce new observed spike trains O, whose dissimilarity with R determines a new set of PSPM prescribed weight changes.This process is repeated for 150 epochs resulting in a network we define as optimized (Figure 1) with output spike trains O defined as the optimized spike trains.
To evaluate whether improvements in similarity between R and O (Spike Train Similarity Measures) are not solely due to network-wide homeostatic adjustments to the weights but rather are due to the precise, neuron-by-neuron synaptic changes specified by our PSPM method, a control condition was established.At the beginning of the learning procedure, a copy of the naive network is instantiated.During the course of optimization, changes were made to this copy's weights according to the following rule: for every synaptic weight change in the naive network, an identical change is made in the control network but at a random synapse (Figure 1).After the 150 epochs, this network is defined to be the control and its output spike trains defined as C. Because the control network benefited from the same number and magnitudes of additive weight changes that the optimized network did, any difference in performance between the optimized and control networks is due entirely to the precisely specified weight alterations called for by our PSPM learning rules.

Pre-Synaptic Pool Modification (PSPM) Learning Rules
At a high level, our PSPM algorithm functions like the mechanisms of synaptic homeostasis (Turrigiano, 2017)- (Keck et al., 2017).In response to differences between desired, reference spike trains R and the spike trains observed in a simulation O, individual weight values are additively modified (Figure 2).
To assess how well the spikes in the reference outputs match up with those of the observed outputs, we optimally pair individual spikes in the ith observed spike train o i (t) with their counterparts in the ith reference spike train r i (t) using a string-matching dynamic program.This dynamic program begins by creating a cost matrix C ∈ R n×m with elements that represent the cost of matching each possible spike pair where n is the number of spikes in r i (t) and m is the number of spikes in o i (t).C kl , therefore, represents the cost of matching the kth spike of r i (t) to the lth spike of o i (t).Cost is computed as the number of timesteps between two matched spikes, or simply approximated as a cap = 15 if a spike remains unmatched.The cost for unmatched spikes a cap caps the number of timesteps allowed between two spikes for them to be considered a match.Therefore, any spike in r i (t) that is more than 15 timesteps or 45 ms away from a spike in o i (t) will not be paired with the spike in o i (t).The dynamic program then finds the lowest cost path through the matrix from position [n, m] to position [1,1].This corresponds to identifying the optimal set of spike pairs P where P = {(k, l): spike k in r i (t) and spike l in o i (t) are paired to minimize total cost}.
Our algorithm then targets the spikes in either spike train that are not paired after application of the dynamic program.For any spike that remains unmatched, there are two cases to be addressed: the case of an extra spike in the observed spike train o i (t) at time t extra , and the case of a missing spike in o i (t) at time t missing .In the case of an extra spike, the algorithm determines which inbound synaptic connections W ij are likely to have most strongly contributed to that spike by identifying which neurons j in the presynaptic pool spiked within z = 10 timesteps of t extra .These strongly contributing synaptic weights W ij between significant neurons j and the neuron with the extra spike i are then stochastically diminished by a small amount by subtracting a small random number δ from W ij .In the case of a missing spike, a similar procedure is followed, save that relevant weights W ij are stochastically increased (Figure 2).The magnitudes of these modifications δ are drawn from a uniform distribution on [0 V, 10 −7 V ].Note that weight matrix values are not allowed to change sign, which would correspond to an excitatory synapse becoming inhibitory or vice-versa.Instead, if an excitatory weight is decreased below zero or an inhibitory weight is increased above zero, the weight is simply set to zero.
In addition, homeostatic weight modifications are made to keep the entire network at a desired level of activity.In the event of excess spiking throughout the observed spike trains O, all weights (not just the excitatory weights) of the observed network are stochastically diminished by subtracting a small random number from each of the weights.If the observed network produces an inadequate number of spikes, all weights are stochastically increased by adding small random numbers to each of the weights.Further, changes to weights are additive increases and decreases, rather than multiplicative increases or decreases in the magnitude of the weight, though we do not allow inhibitory connections to become excitatory, nor vice versa.If the reference spike trains R contain x spikes and the observed spike trains O contain y spikes, changes to each weight value W ij are determined by drawing from a uniform distribution over

Spike Train Similarity Measures
To assess agreement between spike trains after the algorithm was run, we used a modification of the van Rossum distance metric (van Rossum, 2001).In the original van Rossum paper, spike trains are convolved with an exponential window to smooth the signal.To generate a distance measure between the spike trains, absolute or squared differences between convolved signals are then integrated over the duration of the signals.In this project it was desirable to have tolerance on both sides of a spike (before and after), so instead of using an exponential windowing function, we convolved binary spike trains with a Gaussian window as described in Schreiber et.al (Shreiber et al., 2003).Our Gaussian filter had a mean µ = 0 and a standard deviation of σ = 5 √ 2 time-steps or roughly σ ≈ 21 ms.The kernel, or windowing function, K(t) = e −t 2 /2σ 2 = e −t 2 /(100 timesteps 2 ) is convolved with the ith binary spike train s i (t) ∈ {0, 1} T , which is defined for the time interval of the simulation [0, T ] where T is the total number of time steps in a simulation.The result of this convolution is defined as the activity signal for neuron i A total network activity signal can be found by summing individual activity signals a(s i (t), t) over each neuron index i with N the total number of neurons.If S = [s 1 (t), s 2 (t)..., s i (t), ..., s N (t)] ∈ {0, 1} T ×N is the matrix containing each of the output spike trains of a network, then the total activity signal for spike trains S is From these two types of activity signals, two distance metrics are established.First a pairwise distance metric D P (S, R) is constructed by comparing the individual activity signals of two sets of spike trains (with the same number of neurons N ) S, R ∈ {0, 1} T ×N , neuron by neuron over the time interval [0, T ].If s i (t) and r i (t) are the ith binary spike trains for distinct network outputs S and R, then the pairwise distance measure would be We also define an aggregate distance measure between the two sets of spike trains S and R as where A(S, t) and A(R, t) are total activity signals of spike trains S and R respectively.
In addition to these spike train distance measures, we also assessed agreement between spike trains by comparing inter-spike-interval (ISI) distributions, which contain information about the regularity of spiking in a network.For a given spike train s i (t), the inter-spike interval distribution is a series of observations of the number of time steps between adjacent spikes.For instance, if at time t, s i (t) = 1 and the time of the next spike is t + a so s i (t + a) = 1 and s i (t + b) = 0 for 0 < b < a, then the value a is appended to the ISI distribution.For a network wide measure ISI distribution, observations from each neuron are concatenated to produce a distribution with the entire set of observations.To compare ISI distributions between two output spike trains, we conduct a two-sample Kolmogorov-Smirnoff test.This non-parametric statistical test measures the probability that two sets of observations are drawn from the same distribution by evaluating the largest disparity in the respective cumulative probability distributions.A small p-value warrants rejection of the null hypothesis which is that the two observations are drawn from the same distribution.Thus a high p-value indicates a high goodness of fit between two ISI distributions.

Probabalistic Integrate-and-Fire (PIF) Model
In our second set of simulations, reference binary spike trains R were generated from a probabalistic integrate and fire (PIF) model with N = 400 neurons (Karimipanah et al., 2017b)- (Larremore et al., 2012).Like the synapses in the LIF, the strength of PIF synapses are represented by a weight matrix W ∈ R N ×N .To impose sparsity, we set the probability that a synapse is created between any two neurons i and j to p = 10%.In our PIF network, all synapses were excitatory.The weight values were drawn from a uniform distribution on [0, 0.02] so that the mean (nonzero) weight value was 0.01.The maximum eigenvalue λ was calculated and each element of the weight matrix W was subsequently divided by λ so that the maximum eigenvalue of W is 1 and the network operates at criticality (Karimipanah et al., 2017b).Input values I i (t) were drawn from a Poisson distribution with mean and variance µ = σ 2 = 1 which were then multiplied by 0.001 to produce a network wide average spiking frequency greater than 0.003 spikes per time-step or 1 Hz if each time-step is 3 ms, consistent with the overall spike rate of our LIF model.The binary state of neuron i at time t + 1, s i (t + 1), is given in terms of the state of all other neurons at time t: where ξ i (t) are random numbers drawn from a uniform distribution on [0, 1] and Θ is the unit step function.Whereas the LIF network parameters have dimensions corresponding to voltage and current, the PIF model parameters are dimensionless.This PIF network was simulated for T = 50, 000 time-steps to produce the reference spike trains R = [r 1 (t), ..., r i (t), ..., r N (t)] ∈ {0, 1} T ×N .This simulation time was chosen so that there were 5000 or more avalanches (Avalanche Analysis).The PIF outputs R were then broken up into 5 spike trains of 10,000 timesteps each of which served as a set of reference spike trains for the PSPM learning procedure.These spike trains were split so as to reduce the algorithm's input size, improving the speed at which a solution could be obtained.Naive LIF networks were initialized with weight magnitudes drawn from a uniform distribution on [0, 1 × 10 −3 ] with all synapses designated as excitatory.

Avalanche Analysis
To test whether PSPM successfully induced criticality in the optimized LIF networks, avalanche statistics were calculated for each of the output spike trains R, N , O, and C. From a set of spike trains s i (t) for neurons i, a summed network spiking F (t) = i s i (t) was evaluated.Avalanches were defined as events where the summed network spiking F (t) exceeded the 20th percentile of the all summed activity values over the simulation interval [0, T ].An avalanche persists from the time step F (t) first passes above the 20th percentile threshold until F (t) sinks below this threshold.This percentile threshold was chosen so that in a simulation of T = 50, 000 time steps, the number of avalanches exceeded 5000 for the PIF network outputs R. For each avalanche, the size S and duration D were recorded.Size is defined as the number of spikes within an avalanche while duration is the number of time steps for which the avalanche persists.As described in (Karimipanah et al., 2017a), a maximum likelihood method was employed to fit power law probability distributions for the sizes and durations of avalanches.Namely that the size distribution follows P (S) ∼ S −τ while the duration distribution follows P (D) ∼ D −α .At criticality, the average avalanche size < S > and duration D of avalanches also obeys a power law < S >∼ D −β with critical exponent β.By using the calculated critical exponents α and τ from the size and duration probability densities one can calculate a predicted β p = (α − 1)/(τ − 1) which can then be compared to an empirically observed critical exponent β o obtained by fitting average size < S > and duration D observations for the avalanches in the output spike train.

Results: PSPM Learning Improves Spike Train Distance Measures
In simulations with LIF reference and naive networks generated from the distributions described in Initial Network Parameters, optimization improved agreement between the optimized network's output spike trains O and the outputs R of the reference network.Figure 3 (A)-(D) show sample raster plots of reference R, naive N , optimized O, and control C spike trains for a representative trial with uniform initial network parameters (Methods).The reference network outputs R and the naive network outputs N are similar in the overall level of network activity, but close inspection reveals significant differences in the timing of spikes for individual neurons in the network.Because both networks were stimulated with the same input currents, the difference in the reference and reference weight matrices account for the observed disparity in the output spike trains.The optimized spike trains O in Figure 3 (C), however, demonstrate considerably better qualitative agreement with R than N does, indicating improvement due to PSPM learning.The control network C, despite ben-efiting from similar homeostatic adjustments to its weight matrix as the optimized weight matrix (Methods), shows disagreement with R due to the weight changes made at random synapses during learning.This indicates that network performance is sensitive to the prescribed weight changes called for by our PSPM learning rules.
To corroborate these qualitative assessments of spike train similarity for this trial, we generate network activity signals with the Van-Rossum like method described in Spike Train Similarity Measures.Figure 4 (A)-(D) shows the single neuron activity signals for an example neuron, in this case, the first neuron in the network.Again, the optimized activity signal a(o 1 (t), t) shows better agreement with the reference a(r 1 (t), t) than either the naive a(n 1 (t), t) or control a(c 1 (t), t) activity signals do.In addition to single neuron signals, we show the total network activity signals for this trial in Figure 4 (E)-(H) to visually capture network-wide spike train behavior.As before, the optimized network activity A(O, t) shows the greatest agreement with the total activity signal of the reference A(R, t).The pairwise and aggregate distance measures (Methods) for this trial were D P (N , R) = 3.2436 × 10 5 and D A (N , R) = 6.8155×10 7 for the naive network, D P (O, R) = 7.172×10 4 and D A (O, R) = 1.0015×10 7 for the optimized network, and D P (C, R) = 2.8335 × 10 5 and D A (C, R) = 5.3108 × 10 7 for the control network, verifying the qualitative inspection of the spike trains.PSPM learning, therefore, improved agreement with the reference network's outputs for this trial.
Distance measure results for each of the initial network configurations described in Initial Network Configurations are shown in Figure 5.For each initial configuration, 30 trials were conducted.Tables 1 and 2 show these pairwise and aggregate distance measures respectively.For each of the initial network configurations, the distance measures of the optimized networks are lower than the distance measures for the naive and control networks.The fact that PSPM improves performance of the optimized network but not the control demonstrates the importance of neuron-by-neuron precision in synaptic weight modifications during learning (Methods).

Improved agreement in ISI distributions and basic spike statistics dependent on Initial Network Configuration
In addition to reducing distance measures between the optimized spike trains O and the reference spike trains R, the PSPM procedure also improves goodness of fit between inter-spike interval (ISI) distributions.For example, the ISI distributions for each of the four spike trains R, N , O, and C for a sample trial are shown in Figure 6 (A)-(D).Qualitatively, the ISI distribution of the optimized network (C) shows best agreement with that of the reference network (A).The ISI distributions of the naive (B) and control (D) networks differ somewhat from that of the reference network.To quantify similarity between ISI distributions, two-sample Kolmogorov-Smirnov (KS) tests were conducted between the set of reference ISI observations and the observed ISI values from the naive, optimized, and control outputs (Methods).Figure 7 shows the KS test p-values recorded for each of the 30 trials along with spike number mean and variance (computed over neurons in the network).The ISI distribution for O demonstrates the closest fit with that of R as indicated by the high mean Initial Configuration D P (N , R) for sparse, and p KS = 0 ± 0 for naive-half-max.For the control outputs C, p KS = .0002± .0009 for uniform, p KS = .0009± .0042 for gaussian, p KS = 0 ± 0 for sparse, and p KS = .00003± .0002 for naive-half-max.The mean p KS values are highest for the ISI distribution of the optimized network demonstrating best average agreement between the reference ISI distribution and that of the optimized network.This indicates that the PSPM algorithm also improves agreement between ISI distributions of the reference and optimized network.However, the large standard deviations in p KS for each network and initial configuration indicates a large trial to trial variability.Learning the ISI distribution appears sensitive to the initial weights of the naive and reference networks in each trial.
In addition, spike number mean and variance show sensitivity for the initial network conditions (Figure 7).Mean spike numbers were calculated for a collection of spike trains S = [s 1 (t), s 2 (t)..., s i (t), ..., s N (t)] ∈ {0, 1} T ×N by first counting the number of spikes in each individual spike train s i (t) and then taking an average over neurons i. Variance was similarly calculated over the spike numbers in each individual neuron's spike train.While the mean spike number varies sporadically throughout optimization, spike number variance tends to increase with application of PSPM.Interestingly, the gaussian and sparse networks show a larger spike number variance for the optimized outputs O than for the control outputs C, indicating that for these distributions, network-wide homeostatic adjustments do not con- siderably influence the spike number variance but the targeted spike-matching weight adjustments due to PSPM increase spike number variance.

PSPM does not induce agreement between weight matrices, Possible I/O Degeneracy
With the PSPM algorithm's demonstrated success in reducing spike train distance measures and improving goodness of fit of ISI distributions, a plausible expectation may be that the algorithm improves similarity between the reference and optimized weight matrices R, O ∈ R N ×N but an analysis of the resulting weight matrices shows that this is not the case.Summed squared differences were calculated between the elements of naive N , optimized O, and control C weight matrices and the reference weight matrix R.These errors are shown in Table 3 averaged over 30 trials for each of the initial network configurations.The difference between the optimized weight matrices O and the naive weight matrices R is consistently the largest for uniform, gaussian, and sparse and only shows slight improvement for the naive-half-max condition where the elements of the naive and reference weight matrices, N and R respectively, are drawn from distinct probability distributions.A possible explanation for this disagreement between optimized and reference weight matrices despite improvement in distance measures and ISI statistics is that certain input/output maps (from input currents to output spike trains) are degenerate with respect to network weights.That is to say that many possible network weight matrices could accomplish the same input-output mapping.The degeneracy likely increases as the objective of precise spike timing is relaxed and one simply seeks improvement in the accuracy of the input-output map.It remains to be seen whether using PSPM to train a network with a variety of reference input-output observations (with inputs drawn from different distributions) could eliminate the effect of degeneracy by reducing the set of possible networks that could accomplish all of the input-output mappings instead of the single observation used in this study (Discussion).
PSPM used to generate critical behavior in a LIF network.
In the second set of simulations, a probabilistic integrate-and-fire (PIF) networks was generated with the maximum eigenvalue of the weight matrix tuned to 1 (Methods).The PIF network was driven with Poisson inputs and the outputs were evaluated with avalanche analysis.LIF networks were then optimized to reproduce the critical outputs of the PIF, but in a deterministic model net-work.The spike trains of this PIF network were simulated and subsequently used as the reference data for runs of our algorithm.Namely, the 50,000 time steps of the PIF simulation were split into five sets of 10,000 time steps, each of which was used for a run of the PSPM algorithm.The resulting output spike trains were concatenated to provide adequate data for critical avalanche analysis.While the PIF network received external, scaled Poisson inputs (Methods), the LIF network was stimulated using a Gaussian input current distribution identical to that used in our first set of simulations.Poisson inputs were only used for the PIF network on account of its different dynamics.Input currents in the PIF represent probabilities of firing due to external input, whereas input currents in the LIF change the membrane potential.Unlike the first set of simulations, the LIF network synapses were exclusively excitatory, as was the case for the PIF networks described by Karimipanah (Karimipanah et al., 2017b) We investigated the avalanche statistics for the critical PIF reference network and compared it with the avalanche statistics for the naive, optimized, and control LIF network from ten runs of the algorithm (Figure 8).A necessary condition for criticality is that the probability distributions of the size and duration of neural avalanches follow power laws, namely that P (S) ∼ S −τ and P (D) ∼ D −α .In addition to power law distributions, criticality requires a relationship between the size and duration of avalanches < S >∼ D −β , where β = α−1 τ −1 .For a given set of spike trains, one can thus determine whether the predicted and empirically fit values of β agree.As expected, the critical PIF network demonstrates the marked agreement between the predicted and fit β values with absolute difference |β o − β p | = .006(Figure 8 A).Because the PIF network was already tuned to criticality, this agreement is unsurprising.In contrast, the naive LIF networks with weights drawn from a uniform distribution on [0, 1 × 10 −3 ] had fewer avalanches and the avalanches that were observed failed to follow a power law with size or duration as shown in (B).After running our algorithm, the concatenated output spike trains of the optimized and control LIF networks were also subjected to avalanche analysis.Of the three LIF networks, the optimized network demonstrated the best agreement between its predicted and observed values of β with absolute difference |β o − β p | = .064as is evident in Figure 8 (C).This indicates that the weight alterations made during the application of PSPM successfully induced critical activity in the optimized LIF network.Interestingly, however, the control condition, shown in Figure 8 (D) exhibits decent agreement between the predicted and observed β values with |β o − β p | = .128.This surprising agreement is of potential theoretical interest, as it could indicate that non-targeted alterations to the weight matrix of a sort similar to that of biological synaptic scaling are sufficient to induce criticality, providing a potential explanation of the emergence of criticality in biological neural networks which should be explored in future work.

Discussion:
In our attempt to optimize a spiking neural network (without limits being placed on the spike counts or network structure) to match desired spiking behavior under comparable external stimulus conditions, we successfully reproduced statistical features of biological neural networks, including spike train distance measures and ISI distributions.Further, our successful induction of criticality into a LIF neuron network demonstrates the possibility of critical behavior emerging in a purely deterministic, biologically plausible system through simple synaptic scaling procedures.
While past work (Bohte et al., 2002) (Gutig and Sompolinsky, 2006)- (Schrauwen and Campenhout, 2006) has developed numerous approximations of solutions to the exact spike matching problem, each was restricted either in allowable spike number, neuron count, or network structure.PSPM, by contrast, removes the requirement of exactly matching spike timings.This relaxation of the problem allowed us to generate artificial SNNs whose spiking behavior matches the statistical features of a desired set of spike trains.Other learning rules have successfully reproduced spike trains by modifying the weight matrix of a LIF network by comparing VR signals directly (Gardner and Grüning, 2016).Instead of making weight updates based on a convolved van-Rossum distance, PSPM compares binary spike trains directly with a dynamic program.Only after we have trained the network do we quantify statistical matching via a van-Rossum like distance (van Rossum, 2001), agreement.Our measure also improved inter-spike-interval distributions and generated near-critical behavior when trained on PIF outputs.The most immediate application of PSPM is to produce generative models of network activity within the cortex.However, future work may also extend the algorithm's use to applications with artificial neural networks and machine learning.
The largest limitation of the present work is that each of our networks received only a single input-output observation and were trained to generate optimized outputs O as similar as possible to reference outputs R when exposed to identical inputs.Drawing an analogy to computer vision, our procedure in this work was like testing a supervised learning algorithm by presenting a neural network a single image and its corresponding label.While this can lead to potential over-fitting, it also provides less information about the structure of the reference network and its input-output mapping (Abu-Mostafa et al., 2012).Further work could elucidate the performance of PSPM when a network is trained on input data drawn from a variety of distributions.One of our secondary goals was to infer the weight matrix of an artificial neural network given its external inputs and output spike trains.This proved more challenging than anticipated, as numerous (if not potentially infinite) networks may produce spike trains with the same statistical features.Given this degeneracy (Caudill et al., 2009), our algorithm was unable to faithfully replicate the weight matrix of a reference network.By incorporating data sets with a variety of inputs into training, better agreement between resulting weight matrices may be possible.
Further approximations and modifications of this model generation problem could be developed to suit the demands of various theoretical problems in neuroscience and machine learning.Examples of such problems include: (i) Predicting and modeling the neural correlates of a given external stimulus.(ii) Maximizing the efficiency of recurrent spiking networks, particularly in the form of neuromorphic hardware where there is a desire to achieve more processing power per spike (Nahmias et al., 2013).(iii) Mapping mathematical functions (i.e.generalized input-output relationships) onto recurrent spiking neural networks for problems in machine learning and artificial intelligence.
In addition, future efforts may produce heuristics better suited to the above problems.Given the problem's resistance to approximation algorithms, heuristics are the most viable direction for future approaches.In addition, problems related to features of the network structure itself may be of interest, including the development of algorithms for replicating statistical properties of a reference network given only the inputs and spiking output of that network.Such features include in-degree and out-degree distributions, clustering coefficients, and the size distribution and number of cliques present in the network.
A further useful comparison is that of our algorithm's generative models and those of the various types of restricted Boltzmann machines (RBMs).In comparison to traditional RBMs which require the network be organized into distinct layers, our network places no restrictions on the structure of the network, allowing full recurrence.In addition, even modern RBMs such as Sutskever et al's recurrent temporal restricted Boltzmann machine use stochastic neurons which are not intended to be biologically plausible.By contrast, our algorithm works with the substantially more realistic LIF neuron.While both RBMs and our pre-synaptic pool modified (PSPM) networks successfully capture a statistical transformation from input distribution to output distribution, we offer improved correspondence with biological neural networks, and thus provide a better model of such networks' behavior.

Reference Naive
Copy of Naive  In particular, these are the activity signals for neuron 1 in this trial.If r 1 (t), n 1 (t), o 1 (t), c 1 (t) represent the spike trains of the first neuron in the reference, naive, optimized, and control networks respectively then the activity signals a(r 1 (t), t), a(n 1 (t), t), a(o 1 (t), t), a(c 1 (t), t) are plotted in (A), (B), (C), and (D) respectively (Methods).Transitioning from the naive state to the optimized state there is clear improvement in the agreement with the reference spike trains total activity signal a(r 1 (t), t).Similarly, the post-optimization signal a(o 1 (t), t) exhibits greater agreement with the reference signal a(r 1 (t), t) than does the control total activity signal A(C, t) with the reference signal.The fact that

Figure 1 :Figure 2 :Figure 4 :
Figure1: Learning Paradigm (A) Three networks are initialized: the reference network, a naive network and its copy.In addition, a set of input currents is generated and will be used throughout the trial.When simulated, the reference network produces goal spike trains R that the learning rule aims to reproduce in the optimized network.A naive network and its identical copy are created with weights distinct from the reference.(B) The reference and naive networks are simulated with the same input currents to generate outputs R and N and the PSPM learning rules are followed to change the weights of the naive network, henceforth denoted as the pre-optimized network (middle column).For every change in the pre-optimized network, there is an identical change in the copy, except the change is made at a random synapse.(C) This process is repeated for the specified number of epochs (150) of the algorithm.(D) Upon completion, the network in the middle column is considered optimized, the reference (left column) is unchanged and the remaining network (right column) is considered the control.

Figure 5 :Figure 6 :Figure 7 :
Figure4: PSPM improves qualitative agreement between output activity signals and reference activity signals(Methods).Subplots (A)-(D) show single neuron activity signals of network outputs for a representative trial.In particular, these are the activity signals for neuron 1 in this trial.If r 1 (t), n 1 (t), o 1 (t), c 1 (t) represent the spike trains of the first neuron in the reference, naive, optimized, and control networks respectively then the activity signals a(r 1 (t), t), a(n 1 (t), t), a(o 1 (t), t), a(c 1 (t), t) are plotted in (A), (B), (C), and (D) respectively (Methods).Transitioning from the naive state to the optimized state there is clear improvement in the agreement with the reference spike trains total activity signal a(r 1 (t), t).Similarly, the post-optimization signal a(o 1 (t), t) exhibits greater agreement with the reference signal a(r 1 (t), t) than does the control total activity signal A(C, t) with the reference signal.The fact that (A) and (C) show closest agreement indicates successful learning due to PSPM.Total activity signals are plotted in (E)-(H).If the total activity signals for the reference spike trains R, naive spike trains N , optimized spike trains O, and control spike trains C are represented by A(R, t), A(N , t), A(O, t), A(C, t), then these are respectively plotted under labels (E), (F), (C), and (D) (Methods).Again, the optimized total activity signal most closely matches the reference, demonstrating successful learning.Indeed for this trial the aggregate distance measures were D A (R, N ) = 3.24 × 10 5 , D A (R, O) = 7.17 × 10 4 , and D A (R, C) = 2.83 × 10 5 .

Table 1 :
Pairwise distance measures for each of the 4 initial network configurations.Mean and standard deviation reported for 30 trials.Actual distance measure values are 10 5 times those reported above.

Table 2 :
Aggregate distance measures for each of the 4 initial network configurations.Mean and standard deviation reported for 30 trials.Actual distance measure values are 10 7 times those reported above.
p-value for each of the initial network configurations.KS test p-values from the optimized outputs O were p KS = .0756± .1206 for uniform, p KS = .0379± .1007 for gaussian, p KS = .00057± .0028 for sparse, and p KS = .2398± .3217 for naive-half-max.For comparison, KS tests conducted on the naive outputs N resulted in p KS = .00005±.0001 for uniform, p KS = .0007±.0321 for gaussian, p KS = 0±0

Table 3 :
Weight Matrix Differences for each initial network configuration.Mean and standard deviation reported for 30 trials.