Creating synapses with recurrent connections within a single neuron group

[daruwal...@gmail.com]

Hi everyone

I have a neuron group that takes in spikes from a spike generator input stimuli. These spikes are excitatory, and they should increase the membrane potential. I also want to feed the output spikes from the group back in as an input for the next time tick, which should be an inhibitory connection. Below is some sample code to illustrate what I am trying to do:

weights = '''

w : volt

w_fb : volt

'''

synapse = Synapses(input_group, decorrelator_group, weights, on_pre='v_post += w', on_post='v_post += w_fb')

synapse.connect(i = input_idx, j = 0)

synapse.w = '16 * volt'

synapse.w_fb = '-16 * volt'

synapse.post.delay = tau # tau is defined earlier in the code as 1 * ms

Currently, the behavior I am seeing has two sequential output spikes (within defaultclock.dt = tau = 1 * ms apart), and the second inhibitory spike is causing the membrane potential to decrease by more than -16. I have attached a screen capture of my plot that illustrates this effect.

I feel like there is a timing issue related to the evaluation of groups and synapses happening, but I can't figure it out.

Thanks,

Kyle

[daruwal...@gmail.com]

Here is another figure showing this behavior.

[Marcel Stimberg]

Hi Kyle,

does your synapse connect several cells to one target cell? This would
explain it: a single pre-synaptic spike will increase the post-synaptic
once, but if the post-synaptic cell spikes this affects all the synapses
connecting to this cell. I.e., if you have n inputs to one cell, then
the on_post statement will be executed n times. This might seem
counter-intuitive but if you look at it from the point of view of a
synapse, then there is really no other way to implement it that would
make sense.

Having said all that, if I understand correctly your ultimate goal is
that an incoming spike increases the post-synaptic potential for 1ms,
and then reduces it by the same amount, i.e. you want to have a
step-like post-synaptic potential. For this, you should rather use two
pre-synaptic pathways, with different delays. You'd write it like this:

synapse = Synapses(input_group, decorrelator_group, weights,

on_pre={'step_up': 'v_post += w',

'step_down': 'v_post += w_fb'},

delay={'step_down': 1*ms})


If you want to have different delays (i.e. "step length") for individual
connections, then don't include the delay definition in the initializer,
but instead define it afterwards:

synapse.step_down.delay = ...


Best,

Marcel

[daruwal...@gmail.com]

I only have one input neuron connect to this neuron, so there should only be one synapse. Thanks for explaining about the n input case though; it will be very helpful down the road.

My goal is to have one input that is excitatory. The neuron's membrane should rise as input spikes come in. The threshold is random, so the neuron will fire if the random threshold for that time step is below the current membrane potential. At that point, I want to feedback the output spike that has been produced as an inhibitory input that decreases the membrane potential. If there are no input spikes coming in and no output spikes being produced, the potential should stay the same (i.e. there is no leak term).

Thanks,

Kyle

[Marcel Stimberg]

Hi Kyle,

ok I'd need to know more about the model to diagnose this further, then. Is it maybe possible that the membrane potential stays above the threshold for more than one time step and therefore you get more than one inhibitory input? And if you simply want to lower the membrane potential after a spike, couldn't you use the neuron's reset statement for this?

Best,
  Marcel

--
http://www.facebook.com/briansimulator
https://twitter.com/briansimulator
 
New paper about Brian 2: Stimberg M, Goodman DFM, Benichoux V, Brette R (2014).Equation-oriented specification of neural models for simulations. Frontiers Neuroinf, doi: 10.3389/fninf.2014.00006.
---
You received this message because you are subscribed to the Google Groups "Brian" group.
To unsubscribe from this group and stop receiving emails from it, send an email to briansupport...@googlegroups.com.
To post to this group, send email to brians...@googlegroups.com.
Visit this group at https://groups.google.com/group/briansupport.
For more options, visit https://groups.google.com/d/optout.

[Kyle Daruwalla]

Hi Marcel,

Sorry for the long delay, but I wanted to figure out the most convenient way to present you with the necessary code.

Attached are three files — decorrelator_sim.py, ibm_signal_processing.py, and ibm_base.py. Below is a description of each file:

• ibm_base.py: Implements custom functions for equations and includes the create_ibm_neuron function that returns a neuron group that models an IBM TrueNorth neuron given the input parameters to the function
• ibm_signal_processing.py: Includes the create_decorrelator function that creates a decorrelator IBM neuron and the necessary synapses
• decorrelator_sim.py: Generates an input spike train and simulates a decorrelator neuron. The input spike train is a stochastic representation of a number in [0, 1] (i.e. if I want the number 0.2, then the spike train is a series of 1s (spikes) and 0s (lack of spike) as long as num_ticks where the probability of a spike occurring at tick is 0.2)

In the create_ibm_neuron function, you will see that we have some code that thresholds the neuron to a max (up_bound) and min (low_bound) value. These are extremes that the membrane potential of the neuron can reach. Previously, we did not have that code in the model, and the decorrelator neuron exhibited the behavior in the previously attached figures. However, once we added that code in, the neuron behaved correctly. The strange part is that we can look at the plots before we added that code in, and we can see that the membrane potential never comes close to the up_bound or the low_bound values. So, we are baffled trying to figure out why the change fixed things. Moreover, we noticed that changing defaultclock.dt to 0.01 ms would result in the correct behavior (though the input and output streams were not sufficiently decorrelated) even without the up_bound and low_bound code. Currently, we suspect there is a detail related to the timing or order of how Brian performs computations each time step that is eluding us.

FYI, our neuron model should have the same behavior regardless of the length of each time step. The neuron parameter tau = defaultclock.dt globally.

If you would like to run our code, it is available on Github at https://github.com/UW-PHARM/brian2.

Thanks for all your help,

Kyle

[Marcel Stimberg]

Hi Kyle,

I had a look at the code you sent and the code on github and I am still
not entirely clear about the issue. If I run your code with a fixed seed
to have the same input spikes, then I do not see any change in behaviour
regardless of whether I include the up_bound/low_bound code or not.

Re-reading your earlier mails after having had a look at your model, I
also do not quite get why having two subsequent decreases in membrane
potential is problematic/suprising: since your threshold is random, the
membrane potential can still be above threshold after one 16V decrease,
therefore triggering another spike and another 16V decrease.

> FYI, our neuron model should have the same behavior regardless of the
> length of each time step. The neuron parameter tau = defaultclock.dt
> globally.

You mean when you look at things over time steps, not over time, right?
Because you'll get 10 times more spikes when you decrease your dt by a
factor of 10 (you'll compare against the random threshold 10 times as
often). If you wanted to have Poisson spiking with an identical rate
over time, you'd need your threshold to depend on the time step (e.g.
for a PoissonGroup, the threshold condition is "rand() < rate*dt").

Oh, a small thing I noticed: the C++ code for lfsr_rand is incorrect
(and therefore gives different results than Cython), bit22 is defined in
the same way as bit32.

Finally, one general question: did you have a specific issue that made
you fork the entire Brian2 repository instead of developing your
"ibm_lib" as an independent package? We had extensibility in mind when
we designed Brian 2, so if there is something that you wanted to do that
wasn't possible without modifying Brian2 itself, then this would
probably be considered a bug!

Best,
Marcel

[Kyle Daruwalla]

Hi Marcel,

Firstly, thank you for putting so much effort into this.

> If I run your code with a fixed seed to have the same input spikes, then I do not see any change in behaviour regardless of whether I include the up_bound/low_bound code or not.

This is an interesting behavior. As you may be able to tell from the code, the up_bound/low_bound code is part of the positive and negative threshold events. Perhaps we had an error in our standard thresholding code, and when we added the up_bound/low_bound code, we inadvertently fixed the error with the non-up_bound/low_bound thresholds. I will roll back to the older version of the model equations when I get a chance, and I will test that code with a fixed seed like you have.

Just to clarify, when you say the behavior is unchanged, do you mean that across multiple runs of the simulation file, the output is the same. That would be expected, since the seeds are fixed. We are concerned with only a single run of the simulation script. The input stream of spikes should be at some frequency X. We would like the output stream to be at X frequency (with a tight error bound). However, the two streams should not be correlated. So, if I take the total number of spikes in the input stream and divide it by the total number of ticks for the simulation, then that result should be the same as when I divide the total number of spikes in the output stream by the total number of ticks. However, the spikes in the input stream and output stream should not occur during the same time tick. Before we added the up_bound/low_bound code, our output stream did not have the same frequency as our input stream.

> Re-reading your earlier mails after having had a look at your model, I also do not quite get why having two subsequent decreases in membrane potential is problematic/suprising: since your threshold is random, the membrane potential can still be above threshold after one 16V decrease, therefore triggering another spike and another 16V decrease.

That’s correct, but the behavior we saw in the figures was a 16V decrease without a spike occurring. Only the first two 16V decreases in a chain of 3 or 4 could be accounted for bt a output spike occurring.

> You mean when you look at things over time steps, not over time, right?

Yes, we look over time ticks. For us, the length of each time tick is a function of the clock rate we supply to all the computational units (since this is all modeled after a digital CMOS neural substrate).

> Oh, a small thing I noticed: the C++ code for lfsr_rand is incorrect (and therefore gives different results than Cython), bit22 is defined in the same way as bit32.

Thanks for catching this! It will save me a lot of frustrated debugging time :).

> Finally, one general question: did you have a specific issue that made you fork the entire Brian2 repository instead of developing your "ibm_lib" as an independent package? We had extensibility in mind when we designed Brian 2, so if there is something that you wanted to do that wasn't possible without modifying Brian2 itself, then this would probably be considered a bug!

Our group selected Brian because its frontend is straightforward and easy to understand. Furthermore, with numpy readily available, we can use almost all of the array operations we would use in MATLAB which was our previous environment. Most importantly, Brian is a framework that takes models of neurons and their connections defined by mathematical equations, then generates code for simulation. This is important to us, because our end goal is to translate a model of a neural system into synthesizable Verilog code for an FPGA. Specifically, our field is computer architecture, so our interests are to take biologically relevant learning mechanisms and map them to novel hardware platforms. With this in mind, it makes sense to fork our own repository, since we expect to be adding and modifying the rendering and generation backend of Brian quite a lot.

Again, thanks for all your help with our issues.

- Kyle

[Marcel Stimberg]

Hi Kyle,

> Just to clarify, when you say the behavior is unchanged, do you mean
> that across multiple runs of the simulation file, the output is the
> same. That would be expected, since the seeds are fixed.

what I meant was that when I run the same simulation, i.e. a simulation
where I expect identical outcomes due to the fixed seeds, with and
without the up_bound/low_bound code in the thresholds, I do not see any
difference in the results. So unless you can provide a seed value for
which you need to add the up_bound/low_bound workaround, it's hard to
debug this any further.

> Our group selected Brian because its frontend is straightforward and easy to understand. Furthermore, with numpy readily available, we can use almost all of the array operations we would use in MATLAB which was our previous environment. Most importantly, Brian is a framework that takes models of neurons and their connections defined by mathematical equations, then generates code for simulation. This is important to us, because our end goal is to translate a model of a neural system into synthesizable Verilog code for an FPGA. Specifically, our field is computer architecture, so our interests are to take biologically relevant learning mechanisms and map them to novel hardware platforms.

I wasn't questioning your use of Brian, I was just curious why you think
you needed to for the whole repository :)

> With this in mind, it makes sense to fork our own repository, since we expect to be adding and modifying the rendering and generation backend of Brian quite a lot.

This is what I am curious about. We designed Brian 2 so that other
projects could extend it in various ways such as adding new functions,
but even adding support for completely new code generation target
languages or platforms, *without* the need to touch the core of Brian
itself. Therefore, if you come across an issue that you can't solve
without modifying something in Brian, please file a bug in our issue
tracker. From what I see in your project so far, it does not seem to be
necessary to change the Brian core, so you could develop/package your
project independently of Brian (of course depending on it), which should
hopefully make it more sustainable in the long term. Feel free to mail
me directly about this let's discuss it in a bug on our issue tracker.

Best,
Marcel

[Kyle Daruwalla]

Hi Marcel,

> what I meant was that when I run the same simulation, i.e. a simulation where I expect identical outcomes due to the fixed seeds, with and without the up_bound/low_bound code in the thresholds, I do not see any difference in the results. So unless you can provide a seed value for which you need to add the up_bound/low_bound workaround, it's hard to debug this any further.

Ah, thanks for the clarification. I can’t give a you a specific seed for which the results will be different, but I don’t think it is necessary to debug this further. The information and guidance that you have already provided us has been valuable, and I expect we will discover the behavior of the bug (if there even is one!) in more detail as we continue to use that code segment.

> This is what I am curious about. We designed Brian 2 so that other projects could extend it in various ways such as adding new functions, but even adding support for completely new code generation target languages or platforms, *without* the need to touch the core of Brian itself. Therefore, if you come across an issue that you can't solve without modifying something in Brian, please file a bug in our issue tracker. From what I see in your project so far, it does not seem to be necessary to change the Brian core, so you could develop/package your project independently of Brian (of course depending on it), which should hopefully make it more sustainable in the long term. Feel free to mail me directly about this let's discuss it in a bug on our issue tracker.

You are probably right — forking the entire repository was probably unnecessary. Right now, we are very much in the exploratory phase of this project. I expect that as we have a solid understanding of Brian and our goals, then we will be in touch to make sure the produced package makes Brian and our own work as sustainable as possible. Ideally, we will build on top of Brian without modifying it, and if that remains the case, then we will create an independent package that depends on Brian (as you described). Right now, there are no issues with Brian that make that outcome infeasible.

Thanks,
Kyle