Resonances induced by Spiking Time Dependent Plasticity
RResonances induced by Spiking Time DependentPlasticity
Pau Vilimelis Aceituno
Abstract
Neural populations exposed to a certain stimulus learn to represent it better.However, the process that leads local, self-organized rules to do so is unclear. Weaddress the question of how can a neural periodic input be learned and use theDifferential Hebbian Learning framework, coupled with a homeostatic mechanismto derive two self-consistency equations that lead to increased responses to thesame stimulus. Although all our simulations are done with simple Leaky-Integrateand Fire neurons and standard Spiking Time Dependent Plasticity learning rules,our results can be easily interpreted in terms of rates and population codes.
In a more general setting, we would think that forcing a network of neurons to aperiodic repetitive signal should make the network represent the signal better. Further-more, implemented with the well-known rules presented before using simulations andan analytical approach based on two self-consistency equations relying neural activityand its weights.As a introductory example, we can consider a single neuron with two autapses –synapses to the neuron itself – with different delays presented in Fig. 1. One of theautapses has a delay of ms while the other has a delay of ms , on the orders ofmagnitude od non-myelinated central axons [1]. If the neuron is externally forced tofire every ms , the autapse with delay two has its presynaptic spikes coincide withpostsynaptic ones, and thus will get reinforced, while the one with delay of one secondwould not be affected, on account that the pre- and postsynaptic spikes are too far apartto undergo plasticity. This leads to a neuron with a very strong autapse with delay two,which is effectively a resonant dynamical system with frequency Hz .The previous schema is simple to understand but biologically implausible; synapsesare typically very fast and autapses are extremely rare, so we will instead use largenetworks of neurons where the spikes will travel through a large network with shortsynapses instead of few large synapses. This will require having a network whereneurons are active at different points in time, and thus the receptive fields of the neuronsmust correspond to different phases of the input. To be more precise, the model that we will use for simulations and some of our analysiswill be the Leaky Integrate and Fire (LIF) model with a refractory period [2]. In this1 a r X i v : . [ q - b i o . N C ] J un = 200ms t = 100 msd = 200ms d = 200ms d = 200msd = 100ms d = 100ms d = 100ms d = 100mst = 100 mst = 100 ms Figure 1:
Schema of the emergence of resonances : A minimal example of a neu-ral network that induces resonances. The three leftmost schemas show the evolutionmovement of a spike through the autapses and the rightmost plot the resulting structureof the network. A neuronmodel, the state of a neuron at a given time is described by its membrane potential v ( t ) ,which evolves according to the equation τ m dv ( t ) dt = − ( v ( t ) − v ) + u ( t ) , (1)where τ m = 10 ms , v = − mV . u ( t ) is the input to the neuron at time t . When themembrane potential reaches a certain threshold v th = − mV , the neuron ”fires” or”spikes”, meaning that it emits a pulse of current – the spike – that will be sent to otherneurons in the form of a delta function. After firing, the membrane potential is resetto its resting state v and kept frozen at this value for a fixed period of time called therefractory period t ref = 1 ms .The firing of a neuron generates pulses of current that arrive at other neurons, whichin turn update their membrane potentials. If neuron a receives the spikes of neuron b wewill say that there is a synapse going from the second to the first. The receiving neuronis called postsynaptic and the sending neuron is the presynaptic one. This synapse ischaracterized by a weight w ab and a delay d ab which correspond, respectively, to thegain and the latency that the pulse of neuron a goes through before arriving at b .At this point it is useful to discuss the input u ( t ) . As mentioned before, the inputto a neuron is often given by spikes from other neurons, which would give us u ( t ) = (cid:88) k w k δ ( t − t k ) , (2)where w k is the weight of the k th spike and t k its arrival time to the neuron. Note thatwe can also aggregate small contributions from many other neurons forming a smoothfunction of time or a stochastic variable. In this case, u ( t ) might be better described bythe statistics of this variable, such as the expected value and the variance.Finally, we shall note that the activity of neuron populations does not need to be de-scribed by a set of spike times. In many cases we are interested in a coarser descriptionand then it is better to use variables such as the instantaneous firing rate – the averagenumber of spikes per time unit – instead the exact spike times [3]. Furthermore, ifwe are considering time scales much larger than τ m , the exact evolution of the inner2euron state is irrelevant and we can focus on the coarse neuron activity which is givenby x ( t ) = f ( u ( t )) , (3)where f is a smooth function with lower bounded by zero – meaning that the neuroncannot fire less than zero spikes – and upper bounded by t ref , because the neuron needs ms between spikes.Figure 2: Leaky Integrate-and-Fire as a smooth function : Here we plot the firingrate of a LIF neuron as function of the parameters of a Gaussian input. The left plotpresents the average of firing rate for an input presented for 5 seconds and with the twovariables being the mean and standard deviation of a Gaussian input. The middle andright plots correspond to cuts through that surface on the respective directions, with thecenter one being the mean firing rate for a standard deviation of V /s and the rightone corresponding to varying levels of noise for a mean input of − V /s . To ease our computations we will use periodic signals with symmetries among neuronsand only excitatory synapses. Specifically, we will have the external input being u n ( t ) = u m ( t + ∆ t mn ) = u n ( t + T ) (4)where T is the period of the signal, and ∆ t mn simply encodes the difference of timingbetween neurons n and m . Alternatively, we can think of this external input to everyneuron as a periodic signal with neurons having different receptive fields.This external input modifies the activity of the neurons with the equation x n ( t ) = f (cid:32) u n ( t ) + N (cid:88) m =1 W mn x m ( t ) (cid:33) (5)where f is a monotonically increasing function with upper and lower bounds as shownin Fig. 2. This is just a coarse representation of the LIF neurons that represents theprobability of a spike at time t –sometimes called the instantaneous firing rate– whenthe time scale of the input is much larger than the time scale of the neuron.3 Derivation of self-consistency equations
Given the activities of two neurons we can then apply the STDP equation, which givesus ∆ w mn = − γ ( r n ) w mn + (cid:90) ∞−∞ (cid:90) T x m ( t ) x n ( t + ∆ t ) A ± ( w mn ) e ∆ tτS dtd ∆ t. (6)Note that this equation has some parameters that change as the weights evolve, such asthe rate r n or the activities of the neurons. Here we are not interested in the evolutionof those, but rather in their final values, so we will treat this equation together withEq. 5 as a system at equilibrium. This implies that r n is constant, x n ( t ) is fixed and ∆ w mn is zero, thus W mn = 1 γ ( r n ) (cid:90) ∞−∞ (cid:90) T x m ( t + ∆ t ) x n ( t ) A ± ( w mn ) e − | ∆ t | τS dtd ∆ t. (7)Furthermore, since we have the γ ( r n ) homeostatic term we do not need to have boundson the weights, so we will drop that dependency on A ± .If we notice that the timescale of the STDP is very small so that τ S → we cansimplify the integral (cid:90) ∞−∞ x n ( t + ∆ t ) A ± e ∆ tτS d ∆ t = (cid:90) −∞ x n ( t + ∆ t ) A − e ∆ tτS d ∆ t + (cid:90) ∞ x n ( t + ∆ t ) A + e − ∆ tτS d ∆ t ≈ κ S x n ( t ) + κ A ˙ x n ( t ) (8)where ˙ x n is the derivative of x n and the constants κ A and κ S are constants that accountfor the symmetric and asymmetric parts of the STDP kernel, κ A = (cid:90) ∞ A + ( w mn ) e − ∆ tτS d ∆ t + (cid:90) −∞ A − ( w mn ) e ∆ tτS d ∆ tκ S = (cid:90) ∞ A + ( w mn ) e − ∆ tτS d ∆ t − (cid:90) −∞ A − ( w mn ) e ∆ tτS d ∆ t, (9)and therefore we can compute the final weights by w mn = (cid:98) γ ( r n ) (cid:90) T x m ( t ) [ κ S x n ( t ) + κ A ˙ x n ( t )( t )] dt (cid:99) = (cid:98) β S (cid:104) x m , x n (cid:105) + β A (cid:104) x m , ˙ x n (cid:105)(cid:99) (10)where β S , β A correspond to the values of κ S , κ A normalized by γ ( r n ) , and (cid:104) , (cid:105) cor-responds to the correlation over the interval [0 , T ] . The brackets (cid:98)·(cid:99) denote the factthat the weights cannot be negative; even if we were to use balanced networks withinhibitory synapses, the weights would simply decay to zero.4t this point it is convenient to modify the notation and associate to every neuron aspatial variable θ such that we can write the input as u n ( t ) = u ( t, θ n ) = u ( t − cθ n ) (11)where u ( t − cθ n ) is a periodic function that corresponds to the input with differentdelays due to the neurons position.Naturally this leads to a similar formulation in terms of the neuron activity, x n ( t ) = x ( t, θ n ) = x ( t − cθ n ) . (12)Here it is important to note that by using this type of notation where the position is onlyrelevant in terms of time we are imposing a spatial symmetry that must be justified.In a network where the neurons have different connections, we cannot claim thatthe activity of one neuron is equal to the activity of another neuron with a shift. This issimply because their inputs from other neurons might differ. Thus, what we are statinghere is that not only do neurons obtain the same input with a phase change, but alsothat the matrix W is built so that W mn = W m (cid:48) n (cid:48) ∀ m − n = m (cid:48) − n (cid:48) mod N. (13)Here we shall note that most biological networks of neurons are sparse, meaningthat most pairs of neurons are not connected. In small networks, this would imply thatthe input to a single neuron is not given by its expectation, as we will assume here. Themain justification here is the sheer number of connections that every neuron has, whichis typically assumed to be on the order of . [4]. This implies that we can be fairlysure the input to a neuron follows the law of large numbers.Furthermore, the number of neurons can also be assumed to be large, implying thateven if the phases θ n are randomly drawn from a uniform probability, the whole phasespace will be uniformly covered, again by the law of large numbers.Another side effect of the symmetries imposed is that we can now describe theactivity of the whole network by a single equation x ( t ) = f (cid:32) u ( t ) + (cid:90) T w (∆ θ ) x ( t + c ∆ θ ) d ∆ θ (cid:33) . (14)Notice that we are talking about a periodic input, and the phase of each neuron is onlygiven by the phase of that neuron in temporal units that are arbitrary. Thus, we will set c = 1 and now we can exchange the integral by a convolution, so that x ( t ) = f ( u ( t ) + [ x ∗ w ] ( t )) (15)where ∗ is the convolution operator.We can also use our new notation and the normalization of c to rewrite Eq. 10 w (∆ θ ) = β A (cid:90) T x ( t ) ˙ x ( t + ∆ θ ) dt + β S (cid:90) T x ( t ) x ( t + ∆ θ ) dt = β A [ x ∗ ˙ x ] (∆ θ ) + β S [ x ∗ x ] (∆ θ ) . (16)5q. 15 and Eq. 16 determine the final weights and activity of a neural networksubject to the input u ( t ) after STDP has modified the synapses.Now we would like to know if they can implement a network adaptation to induceresonance to a given input signal u ( t ) . Neural Activity Before STDP
50 100 150 200 250 300
Time (ms) N eu r on I nde x Neural Activity After STDP
50 100 150 200 250 300
Time (ms) N eu r on I nde x -40 -20 0 20 40 Relative Time (ms) A v e r age R a t e Before STDPAfter STDP
Connectivity Matrix
100 200 300 400 500
Neuron Index N eu r on I nde x Connectivity Matrix
100 200 300 400 500
Neuron Index N eu r on I nde x -1000 0 1000 Input Phase Distance A v e r age C onne c t i on W e i gh t Synaptic Weight Distribution
Before STDPAfter STDP
Figure 3:
Evolution of a circulant peak of activity : We simulated 400 neurons andfeed each one with a sinusoid with phases uniformly distributed across the full circlefor all neurons. The upper left and upper center plot correspond to raster plots of neuralspikes, and the upper right is the neural activity averaged over ten periods and over allneurons with the phase being centered with respect to their inputs. The lower left andcenter plots are the matrix weights where brighter color represents higher weights, andthe right lower plot is the average vector w . After repeating the input many times, theSTDP changes the weights which take the shape of a sinusoid (lower row), and theactivity is shifted to be slightly advanced and higher (upper row). The low effect of theweight is simply a result of the parameters of the neurons and weights, which do notallow large weights without generating instabilities through the STDPThe first problem in Eq. 15 is that we have the non-linear function f ( · ) , renderingclosed-form solutions complicated. The first approach to deal with such situations is to6inearize Eq. 15, x ( t ) = (cid:37) (cid:34) u ( t ) + (cid:90) T w (∆ θ ) x ( t + c ∆ θ ) d ∆ θ (cid:35) (17)where (cid:37) is the derivative of f around the mean input to a single neuron, (cid:37) = ∂f∂u (cid:32) ¯ u + ¯ x (cid:90) T w (∆ θ ) d ∆ θ (cid:33) (18)with ¯ x, ¯ u are the mean activity and input respectively. Note that the derivative of f is taken with respect u , although it can be taken with respect to any of the inputs to aneuron. The means can be used because the linearization implies that the fluctuationsare small, but in any case the main point is that (cid:37) is constant.Now that the system consists only of linear equations with convolutions whosesolutions are periodic, we can notice that all the operations involved can be written inthe Fourier domain, F [ x ] ( ω ) = (cid:37) ( F [ u ] ( ω ) + F [ x ] ( ω ) F [ w ] ( ω )) F [ w ] ( ω ) = ( β A ω + β S ) F [ x ] ( ω ) , (19)which we can merge to obtain F [ x ] ( ω ) = (cid:37) (cid:16) F [ u ] ( ω ) + F [ x ] ( ω ) ( β A ω + β S ) (cid:17) (20)This last equation has the advantage that it gives us a direct connection between theinput and the activity. More specifically, if we put a single sinusoid with period T asan input, the Fourier transform consists of two deltas at ω = ± π . Thus, the previousequation is zero everywhere except at the position of the deltas, and there it yields adepressed cubic equation that can be solved analytically.Even though we have an analytical solution, it is worth noticing that the valuesof (cid:37), β A , β S must be estimated numerically as they depend nonlinearly on the inputstatistics to each neuron. Furthermore, there are severe constraints on the parametersof the STDP that arise from the stability of the plasticity and depend on the numberof neurons and synapses. Thus, in the simulations presented in Fig. 3 we observe thatthere is only a small modification of the input activity.That limitation does not imply that our previous results were useless. First, becausethe results relate the general shape of the STDP kernel to qualitative modifications onthe activity, namely the shift in signal phase linked to the strength and sign of β A thatarises from the two solutions for ± ω . Second, because we can obtain the shape of w as a sinusoid which in our case leans towards the asymmetry and thus gives us weightsthat are sinusoids with a phase difference in the interval (cid:2) , π (cid:3) , mostly toward the endof it as shown in Fig. 3. The second approach to this problem to deal with Eq. 15 is to take advantage of the par-ticular the nonlinearity of f ( · ) – with upper and lower saturation – and take a signal –7 eural Activity Before STDP
100 200 300 400 500 600
Time (ms) N eu r on I nde x Neural Activity After STDP
100 200 300 400 500 600
Time (ms) N eu r on I nde x -100 -50 0 50 Relative Time (ms) A v e r age R a t e Before STDPAfter STDP
Connectivity Matrix
50 100 150 200
Neuron Index N eu r on I nde x Connectivity Matrix
50 100 150 200
Neuron Index N eu r on I nde x -600 -400 -200 0 200 400 600 Input Phase Distance A v e r age C onne c t i on W e i gh t Synaptic Weight Distribution
Before STDPAfter STDP
Figure 4:
Evolution of a circulant peak of activity : We simulated 500 neurons withrandom connections and imposed an activity u ( t ) ∝ exp (cid:16) − t τ u (cid:17) where the time con-stant τ u is on the order of magnitude of the leak time constant so that the activity ofneurons with phase t − cδθ affects the activity of neurons at t . This activity modifiesthe initially random network into a circulant-like network, where every neuron con-nects only to its left neighbors (lower panels), which in turn reduces the noise in thenetwork by decreasing the amount of input at non-active times and increasing the inputat the time when the neurons should be active (upper panels). u ( t ) – that is encoded as a sequence of very sharp symbols with small overlapping, andthen let STDP converge into a neural network where each of those symbols ”prepares”the activation of the next.The intuition here is that we have a signal consisting in two levels of activity: thepeak and the background activity. To improve the signal we only need to push thosetwo levels apart, meaning that the neurons that would be active should have a rate ashigh possible and those that are not meant to be active should be as silent as possi-ble. This can be seen as a long sequence of symbols where only a small populationof neurons should be active for every symbol and with the populations having verylittle overlap; then, the causal associations promoted by STDP will create a network ofneural populations where the populations associated to each symbol have strong con-nections towards the subsequent population and thus prepare it to receive an input andfire.In more practical terms, we are using a simple input which consists of a very nar-8ow peak of current concentrated on very few neurons, and a constant backgroundactivity corresponding to the probability of a neuron firing without input. This peakcorresponds to a mollified Dirac of activity, which we can plug into Eq. 16, w ( θ ) = (cid:106) β A (cid:104) ˜ δ ( θ ) ∗ a ˜ δ ( θ − (cid:15) ) − ˜ δ ( θ ) ∗ a ˜ δ ( θ + (cid:15) ) (cid:105) + β S (cid:104) ˜ δ ( θ ) ∗ ˜ δ ( θ ) + c (cid:105)(cid:107) (21)where ˜ δ is the mollified Dirac function, whose derivative is approximated by a ˜ δ ( θ − (cid:15) ) − a ˜ δ ( θ + (cid:15) ) , with (cid:15) being a very small number and a being a scaling that dependson the exact mollification, and the c corresponds to the probability of two neuronsrandomly firing at similar times thus strengthening their connections by β S . Assumingthat the mollified version behaves similarly to a Dirac, we can apply the convolutionsand obtain w ( θ ) ≈ (cid:106) β A a ˜ δ ( θ − (cid:15) ) − β A ˜ δ ( θ + (cid:15) ) + β S ˜ δ ( θ ) + cβ S (cid:107) . (22)In our case β A (cid:29) β S > , and the weights cannot be negative, the shape of theweights remains constant with value β S c for most θ , then a peak with height aβ A at − (cid:15) and a drop to zero at (cid:15) . Furthermore, with those weights the activity gets reinforcedbecause the neurons at phase θ promote the activity of neurons at phase θ + δθ , thus theactivity increases due to neighboring neurons is mostly arriving at the time when theneuron is getting the high level input. This leads to an increased signal to noise ratio,which is shown in Fig. 4 In this chapter we derived solutions for the evolution of synaptic activity under periodicinputs that use differential Hebbian learning under a self-consistency set-up. In thisway we have shown that the features of Echo State Networks that were important inthe first part of the thesis can emerge out of synaptic plasticity rules that are well knownin the field of neuroscience. In particular, we saw that increased signal-to-noise ratiosappear in both the linear and the sparse solution, showing that the signal processingnotions can provide neuroscience insight.The networks proposed here differ from the Echo State Networks used in earlierchapters in the sense that the two solutions presented require very specific receptivefields. Therefore, we cannot use the input as we did in Echo State Networks, but mustinstead project it into the network in such a way that different neurons are active atdifferent times.Although the self-consistency equations used before can in principle be solved nu-merically, the fact that the parameters β A , β S are unknown makes such numerical ap-proaches problematic. However, the limitations outline here are general to theoreticalneuroscience, as biological neural networks are typically heterogeneous and estimatingthe parameters of its components is an arduous task. The results here should thereforebe taken qualitatively, in the sense that we would expect STDP to adapt networks ofneurons to their input in such a way as to increase the signal to noise-ratio and maybeinduce an advancement of the neural activity with respect to the input.9nother remark to this chapter is that we have only proposed the two simplest solu-tions compatible with Eq. 15 and Eq. 16, but naturally there are other possible options.A straightforward approach would be to increase the multiplicity of the solution, forinstance by having multiple isolated sinusoids, or multiple peaks of activity per neuron,in such a way that the general shape of the solution stays the same but neurons haveactivities composed by multiple receptive fields. Acknowledgement
P.V.A is supported by the Supported by BMBF and Max Planck Society
References [1] Harvey A Swadlow and Stephen G Waxman. Axonal conduction delays.
Scholar-pedia , 7(6):1451, 2012.[2] L. Lapique. Recherches quantitatives sur lexcitation electrique des nerfs traitecomme une polarization.
Journal de Physiologie Pathologie Genetique , 1907.[3] Wulfram Gerstner, Werner M Kistler, Richard Naud, and Liam Paninski.
Neuronaldynamics: From single neurons to networks and models of cognition . CambridgeUniversity Press, 2014.[4] Jeff Hawkins and Subutai Ahmad. Why neurons have thousands of synapses, atheory of sequence memory in neocortex.