STDP weights converging to zero

[touches]

Description of problem

I am working on implementing the paper: Pattern recognition using Spiking Neural Networks. This is also posted under the Projects category. As per the paper, STDP is the learning rule applied. I have used the Brian2 STDP example from the official Brian2 tutorials. The problem is on running the script, the weight converges to 0 despite getting the correct firing pattern of neurons as given in the paper

Minimal code to reproduce problem

taupre = taupost = 20*ms

wmax = 1

Apre = 0.01

Apost = -Apre*taupre/taupost*1.05

stdp_eqs = '''
         w : 1
         dapre/dt = -apre/taupre : 1 (event-driven)
         dapost/dt = -apost/taupost : 1 (event-driven)
       '''

stdp_pre = '''
      v_post += w * mV
      apre += Apre
      w = clip(w+apost, 0, wmax * mV)
'''

stdp_post = '''
        apost += Apost
        w = clip(w+apre, 0, wmax * mV)
'''
# Inp -> Hidden
Syn_inp_hid = b2.Synapses(izh , 
                          izh_hidden , 
                          model = stdp_eqs , 
                          on_pre = stdp_pre , 
                          on_post = stdp_post , 
                          method = 'euler' , 
                          name = 'Syn_inp_hid')

# Hidden -> Output
Syn_hid_out = b2.Synapses(izh_hidden , 
                          izh_output , 
                          ''' w : 1''' , on_pre = 'v += w * mV' , 
                          method = 'euler' , 
                          name = 'Syn_hid_out')



# MAKE CONNECTIONS
Syn_inp_hid.connect()
Syn_hid_out.connect()

Syn_inp_hid.w = 'rand() * 0.8'
Syn_hid_out.w = 'rand() * 0.8'

========

What you have aready tried

On setting wmax to very high values like say 500, I get the following actiivty from the weights

Expected output (if relevant)

This is the graph taken from the above paper, the weights can be seen climbing to higher values and 0.8 is considered to be the max value specified by the authors.

Actual output (if relevant)

The STDP plot is at the bottom. This is the graph generated when I run the script.

The architecture is feedforward (i.e.) Input -> Hidden Layer -> Output. I have made the Synapses connections accordingly. Any suggestion on what to look for would be appreciated.

[assergazina]

Good afternoon,

I was reading the article and I also saw the code, I do not know how much it affects to get the same value as in a paper. I found out that the parameters are used different values than actual values which were in the article.

And I am stack with how you got “dv/dt” from the (1)-(2) equations.

[assergazina]

My previous question was about parameters. And I tried to change some parameters and got this output. It seems the same figure on paper (Fig.1)

But I’m curious about even I tried those values but the result gives different outputs. Why?
Even though, here I also want to mention “Voltage range”, on the article range lies between -60 and 40

[mstimberg]

Not quite sure this is really an issue for the support category – I think Brian is doing what it is supposed to do, it is just a question of setting parameters, etc.? Maybe more appropriate to discuss this in the Pattern recognition in Spiking Neural Nets using Brian2 topic.

Or maybe I misunderstood the problem? The code posted in the first post is not complete, so it is hard to see what is going on. If weights go to zero that normally means that the post-synaptic neurons does not spike (or very little), in a typical STDP rule the weights only increase in the on_post rule.

[touches]

@mstimberg: For posting this in the support category, I read through some of the posts in the support category and it had questions on implementations, so I thought this fits in as an incorrect implementation of the STDP rule which includes wrong parameters. The support category also reads:

“This category is for asking questions of the community and developers when something isn’t working right.”

I am not sure of where the boundary lies on what type of questions to post in support. Let me know if this is a mistake to post here and I will refrain from making such posts

And for the code, the format of the support category wanted minimal code for reproducing the error, I didn’t want to paste the full code risking a response saying not to post lots of code The code snippet given actually implements the Synapses and STDP, the rest of the code is the NeuronGroup definition.

This is the extended code if any help can be offered:

# Izhikevich Equation
izh_eqs = '''

dv/dt = (0.04/ms/mV) * v**2 + (5/ms) * v + 140*mV/ms - u + I : volt
du/dt = a * (b * v - u)                                : volt/second

a                                                : 1/second
d                                                : volt/second
I  = digit_ta_input(t , i): volt/second


'''

c = -65 * mV

reset_iz = '''

v = c
u = u + d

'''


stdp_eqs = '''
             w : 1
             dapre/dt = -apre/taupre : 1 (event-driven)
             dapost/dt = -apost/taupost : 1 (event-driven)
           '''

stdp_pre = '''
          v_post += w * mV
          apre += Apre
          w = clip(w+apost, 0, wmax * mV)
'''

stdp_post = '''
            apost += Apost
            w = clip(w+apre, 0, wmax * mV)
'''



# INPUT LAYER
izh = b2.NeuronGroup(inp_neu ,
                     izh_eqs , 
                     threshold = 'v > 30 * mV' , 
                     reset = reset_iz , 
                     method = 'euler' , 
                     name = 'Izh')


# izh.I = np.mean(inp_data) * mV/ms
izh.a = 0.02/ms
izh.d = 8 * mV/ms
izh.v = c


# HIDDEN LAYER
izh_hidden = b2.NeuronGroup(hidden_neu ,
                            izh_eqs , 
                            threshold = 'v > 30 * mV' , 
                            reset = reset_iz ,
                            method = 'euler' ,
                            name = 'HiddenLayer')


# OUTPUT LAYER
izh_output = b2.NeuronGroup(output_neu ,
                            izh_eqs , 
                            threshold = 'v > 30 * mV' , 
                            reset = reset_iz ,
                            method = 'euler' ,
                            name = 'OutputLayer')



                          


# Inp -> Hidden
Syn_inp_hid = b2.Synapses(izh , 
                          izh_hidden , 
                          model = stdp_eqs , 
                          on_pre = stdp_pre , 
                          on_post = stdp_post , 
                          method = 'euler' , 
                          name = 'Syn_inp_hid')

# Hidden -> Output
Syn_hid_out = b2.Synapses(izh_hidden , 
                          izh_output , 
                          ''' w : 1''' , on_pre = 'v += w * mV' , 
                          method = 'euler' , 
                          name = 'Syn_hid_out')



# MAKE CONNECTIONS
Syn_inp_hid.connect()
Syn_hid_out.connect()

Syn_inp_hid.w = 'rand() * 0.8'
Syn_hid_out.w = 'rand() * 0.8'



spike_mon_izh = b2.SpikeMonitor(izh)
state_mon_inp = b2.StateMonitor(Syn_inp_hid , ['w'] , record = True)

plt.plot(state_mon_inp.t/ms , state_mon_inp.w[0 : 25].T / wmax)
ax.set_xlabel('Time')
ax.set_ylabel('Weight/wmax')
ax.set_title("STDP Weights")
plt.tight_layout()

I can see the Neuron spikes from the raster plot and when I set the wmax value high, there is activity so am not really sure which parameter value to investigate.

[touches]

@assergazina: Thank you for offering to read through and giving some insight into the problem. The parameters that are in the paper is for the Izhikevich neuron model. They have not specified the parameter values for STDP.

The dv/dt seen in the code is derived by rearranging and expanding the equation given in the paper and the constant values are empirically obtained via experiments by E Izhikevich himself in his paper on the Izhikevich neurons. This is the standard form of writing this framework of neurons

[kernfel]

@touches

stdp_pre = '''
      v_post += w * mV
      apre += Apre
      w = clip(w+apost, 0, wmax * mV)
'''

You have a unit problem. Your w is dimensionless, but you’re clipping between 0 and 1 * mV which presumably evaluates to 1e-3… which is exactly the range you’re getting, as far as I can see from your w/wmax plot.

@mstimberg, should clip() do unit checks?

[mstimberg]

@touches

I am not sure of where the boundary lies on what type of questions to post in support. Let me know if this is a mistake to post here and I will refrain from making such posts

Yea, I don’t really know either yet For me, the general idea is: Support for questions where you have model/equations/parameters, but it is not working in Brian as you expect. If there are doubts about the underlying model/equations/parameters, then this should rather go into Science, Projects, and Showcases. I thought this rather fell into the latter category, but @kernfel’s comment points to a Brian issue as well, so maybe it should indeed be a Support question. Either way, let’s continue the discussion here for now.

@touches

And for the code, the format of the support category wanted minimal code for reproducing the error, I didn’t want to paste the full code risking a response saying not to post lots of code

This is indeed a tricky issue as well We certainly appreciate not having a full dump of code, but the “minimal code” should be runnable and somehow reproduce the problem. In your concrete example, I guess you could for example show your issue with just the hidden layer, removing the output layer?

@kernfel

You have a unit problem. Your w is dimensionless, but you’re clipping between 0 and 1 * mV which presumably evaluates to 1e-3 … which is exactly the range you’re getting, as far as I can see from your w/wmax plot.

well spotted.

@kernfel

@mstimberg, should clip() do unit checks?

It should and I am a bit confused why it does not! I’ll investigate, this seems to be a genuine

[mstimberg]

@mstimberg

It should and I am a bit confused why it does not! I’ll investigate, this seems to be a genuine

It is indeed a bug, I opened a github issue.

[mstimberg]

@assergazina

Can I send to your email address one file about Units(only 16 pages)?

I’d rather prefer to discuss things here (can be a new topic), since other users can benefit from it as well. Also, “only 16 pages” sounds a bit scary

Just to make sure: you already had a look at the documentation and the tutorial that explain Brian’s unit mechanism? Units are indeed essential, but you can also define variables as dimensionless (except for things like time constants that always need units of time). This is quite common for simple models, several examples in the Brian documentation (e.g. Brette 2004) use this approach.

[touches]

Thank you for spotting that, after assigning the units for the weights, I am able to see some activity from the synapses.

This has wmax set to 0.8.

[touches]

@mstimberg

This is indeed a tricky issue as well We certainly appreciate not having a full dump of code, but the “minimal code” should be runnable and somehow reproduce the problem. In your concrete example, I guess you could for example show your issue with just the hidden layer, removing the output layer?

The hidden layer is a NeuronGroup definition and am not sure how it would have represented the problem in the STDP parameters . The example required the STDP equations for kernfel to spot the issue and the Synaptic connections gives some context to it. I believe to my credit, my post was concrete in terms of explaining the concern I was facing while still following the guidelines. In the future any posts related to modelling or parameters will be in Science, Projects, and Showcases. I did see a conducatance based synapse implementation post in Support, hence the motivation to post it here

[mstimberg]

@touches

The hidden layer is a NeuronGroup definition and am not sure how it would have represented the problem in the STDP parameters .

What I meant that there was no need for the output layer and the connections to it. But let’s not overanalyze this code example

[touches]

Yeah seeing that posts by other users are not overanalyzed or rather not analyzed at all as much as this post . TBH this is kind of upsetting when the post is made to learn and contribute as a way of expressing genuine thanks and gratitude for this wonderful Brian simulator library. I hope different users dont have different standards. Sadly, for Brian there is not much user support available in any outside forums like Stackexchange for example and hope my future posts here is not criticized and commented as to whether it is filled under the right category or has the proper code representation etc. Before posting I make sure to thoroughly check the standards and guidelines and only after that I post here or in any forums.

[mstimberg]

Hi @touches. I am sorry to hear that you found our interaction upsetting… As I said in my earlier posts, we just created this forum and we are still figuring things out, this was the main reason to propose a category change (one of the big difference to our previous forum is that we now have categories in the first place). I also did not mean to criticize your code, I only pointed out that the initial code was not runnable, and therefore making it hard for others to understand the problem (even though @kernfel had a keen eye and spotted the unit issue). When I suggested that you don’t need the output layer to show your issue, I was just trying to be helpful and give a suggestion for future issues. Coming up with good minimal examples is hard and often half the answer to a problem! I did not intend to discuss the code question any further, but you apparently misunderstood my comment so I tried to clarify things. When I said “but let’s not overanalyze this code example” I only wanted to say (apparently in a clumsy way): don’t worry about it, there’s no need to polish the code example any further.

But I also had the impression that from your reply to @kernfel that your basic issue why you posted this support question was now solved and now you are only trying to find the parameters that they used in the paper (which they unfortunately did not include)? So I am genuinely confused what kind of feedback you are currently looking for.