PoissonGroup spike timing

[Daniel Bliss]

Hi all,

I'm replicating using brian2 a model that was originally written in C++.  Recently my brian2 version has produced some results that don't match with the C++ version, despite my using identical parameter settings to that version.  I've been discussing this with my collaborator, who wrote the original code, and he had this to say:

"I have thought that another potential source for discrepancies could be the Poisson implementation. In my code I used an implementation based on time events (ISIs follow an exponential distribution, which means that Poisson events are randomly distributed within each dt). This compared to the implementation based on number of events in each dt (if it's what you're using), can have an impact on final results. This is because one vs the other implementation has a distinct effect on Poisson-dependent synaptic summation, which is decreased by randomizing the time events (otherwise, they happen totally in sync in each dt)."

I'm using the default PoissonGroup settings.  Does brian2's PoissonGroup use an implementation based on time events or number of events?  If the latter, is there a way for me to achieve the former in brian2?  I'm not sure what I can do with the 'when' argument to PoissonGroup.

Many thanks,
Dan

[Dan Goodman]

Hi,

Yes it's quite likely that differences in the Poisson implementation
could lead to differences in results. The way PoissonGroup works is that
at each timestep for each neuron it computes the condition rand()<p for
a certain value of p, and if this condition is true it fires a spike.
This is a good approximation when dt is sufficiently small that dt*dt is
negligible, i.e. when the probability of getting two spikes from a
single neuron within the same window dt is small enough not to worry about.

Also note that all Brian simulations are run on a fixed width clock, so
for example all spikes happen at a time that is a multiple of dt. If you
are comparing to code that is event driven (i.e. spikes can occur at
arbitrary times) then you could get substantially different results in
some cases.

If you want to see if this sort of thing is causing your differences,
try running a few simulations with a very small value of dt.

Also note that you can use PoissonInput in some cases, and this will
replicate the effect of having multiple spikes within a single time
window dt. Note that this is a new class, not available in the released
versions of Brian, only by downloading it directly from the github
repository. Documentation is also a little sparse:

http://brian2.readthedocs.org/en/latest/reference/brian2.input.poissoninput.PoissonInput.html

Dan

> --
>
> ---
> You received this message because you are subscribed to the Google
> Groups "Brian Development" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to brian-developm...@googlegroups.com
> <mailto:brian-developm...@googlegroups.com>.
> For more options, visit https://groups.google.com/d/optout.

[Daniel Bliss]

Dan,

Many thanks for the reply.  Right now I have the following:

simulation_clock = brian2.Clock(dt=0.02 * brian2.ms)
fext = 1800.0 * brian2.hertz
pe_pfc = brian2.NeuronGroup(NE, . . . , clock=simulation_clock)
pg_e_pfc = brian2.PoissonGroup(NE, fext, clock=simulation_clock)
cpe_pfc = brian2.Synapses(pg_e_pfc, pe_pfc, . . . , connect='i == j', clock=simulation_clock)

So, each neuron receives input from one poisson source, which is firing at 1800 Hz.  This firing rate implies that the neuron receives, on average, one spike every 0.56 ms.  With dt set to 0.02 ms, do I have to worry about the fact that this simulation, unlike the original, is not event-driven?

Regardless of the answer to that question, it seems that my code could benefit considerably from a switch from PoissonGroup to PoissonInput.

Thanks again,
Dan

> <mailto:brian-development+unsub...@googlegroups.com>.

[Daniel Bliss]

Also, for what it's worth, I found this in the Brian2 docs:

In the future, the implementation of PoissonGroup will change to a more efficient spike generation mechanism, based on the calculation of inter-spike intervals. Note that, as can be seen in its equivalent NeuronGroup formulation, a PoissonGroup does not work for high rates where more than one spike might fall into a single timestep. Use several units with lower rates in this case (e.g. use PoissonGroup(10, 1000*Hz) instead ofPoissonGroup(1, 10000*Hz)).

https://brian2.readthedocs.org/en/2.0b2/user/input.html

[Dan Goodman]

Hi,

Yeah, PoissonInput is probably going to be better for you than
PoissonGroup. If you do want to use PoissonGroup, then it would be
better to use more neurons each with a lower firing rate. 1800 Hz is
unrealistically high for a single neuron, you might want to break that
into 100 neurons each firing at 18 Hz say, and connect each of those 100
neurons to the same output neuron. Note though that still, PoissonInput
is better than PoissonGroup for this.

0.02*ms is very small so that should cover you. I would try PoissonInput
first though and see if that makes the difference.

Incidentally, with such high input firing rates, you could probably
approximate it very well with a stochastic differential equation. See
Fourcaud and Brunel 2002 (in Neural Computation).

Dan

On 10/04/2015 23:02, Daniel Bliss wrote:
> Also, for what it's worth, I found this in the Brian2 docs:
>
> In the future, the implementation of PoissonGroup

> <https://brian2.readthedocs.org/en/2.0b2/reference/brian2.groups.poissongroup.PoissonGroup.html#brian2.groups.poissongroup.PoissonGroup> will

> change to a more efficient spike generation mechanism, based on the
> calculation of inter-spike intervals. Note that, as can be seen in its
> equivalent NeuronGroup

> <https://brian2.readthedocs.org/en/2.0b2/reference/brian2.groups.neurongroup.NeuronGroup.html#brian2.groups.neurongroup.NeuronGroup>formulation,
> a PoissonGroup
> <https://brian2.readthedocs.org/en/2.0b2/reference/brian2.groups.poissongroup.PoissonGroup.html#brian2.groups.poissongroup.PoissonGroup> does

> > <mailto:brian-developm...@googlegroups.com>.

> > For more options, visit https://groups.google.com/d/optout

> <https://groups.google.com/d/optout>.

>
> --
>
> ---
> You received this message because you are subscribed to the Google
> Groups "Brian Development" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to brian-developm...@googlegroups.com

> <mailto:brian-developm...@googlegroups.com>.

[Daniel Bliss]

"1800 Hz is unrealistically high for a single neuron" -- Of course it's unrealistically high in terms of the biology, but, computationally, does it make a difference whether I split this into 100 neurons or keep it as one (as in, does it make a difference to the neurons that receive the spikes)?  I don't care about the biological realism of this "neuron" -- I don't record from it and it doesn't receive any input; it's meant to represent generalized input to my network from "other cortical areas."

         > <mailto:brian-development+unsub...@googlegroups.com>.

         > For more options, visit https://groups.google.com/d/optout

        <https://groups.google.com/d/optout>.

--

---
You received this message because you are subscribed to the Google
Groups "Brian Development" group.
To unsubscribe from this group and stop receiving emails from it, send

an email to brian-development+unsubscribe@googlegroups.com
<mailto:brian-development+unsub...@googlegroups.com>.

For more options, visit https://groups.google.com/d/optout.

--

--- You received this message because you are subscribed to the Google Groups "Brian Development" group.

To unsubscribe from this group and stop receiving emails from it, send an email to brian-development+unsubscribe@googlegroups.com.

[Dan Goodman]

The only difference it makes splitting it into 100 neurons is that it
reduces the probability of two spikes ending up in the same bin, and so
improves the approximation we use. I only mentioned biological realism
because that explains why we used the approximation we did: for
realistic neurons it's a good approximation. However, your use case of
summed activity of many neurons is precisely why we made PoissonInput so
I'd strongly recommend using that.

Dan

> <https://brian2.readthedocs.__org/en/2.0b2/reference/brian2.__groups.poissongroup.__PoissonGroup.html#brian2.__groups.poissongroup.__PoissonGroup

> <https://brian2.readthedocs.org/en/2.0b2/reference/brian2.groups.poissongroup.PoissonGroup.html#brian2.groups.poissongroup.PoissonGroup>>
> will
> change to a more efficient spike generation mechanism, based on the
> calculation of inter-spike intervals. Note that, as can be seen
> in its
> equivalent NeuronGroup

> <https://brian2.readthedocs.__org/en/2.0b2/reference/brian2.__groups.neurongroup.__NeuronGroup.html#brian2.__groups.neurongroup.NeuronGroup
> <https://brian2.readthedocs.org/en/2.0b2/reference/brian2.groups.neurongroup.NeuronGroup.html#brian2.groups.neurongroup.NeuronGroup>__>formulation,
> a PoissonGroup
> <https://brian2.readthedocs.__org/en/2.0b2/reference/brian2.__groups.poissongroup.__PoissonGroup.html#brian2.__groups.poissongroup.__PoissonGroup

> <https://brian2.readthedocs.org/en/2.0b2/reference/brian2.groups.poissongroup.PoissonGroup.html#brian2.groups.poissongroup.PoissonGroup>>
> does
>
> not work for high rates where more than one spike might fall into a
> single timestep. Use several units with lower rates in this case
> (e.g.
> use PoissonGroup(10, 1000*Hz)instead ofPoissonGroup(1, 10000*Hz)).
>

> https://brian2.readthedocs.__org/en/2.0b2/user/input.html
> <https://brian2.readthedocs.org/en/2.0b2/user/input.html>
> <https://brian2.readthedocs.__org/en/2.0b2/user/input.html

> <https://brian2.readthedocs.org/en/2.0b2/user/input.html>>
>
> On Friday, April 10, 2015 at 2:51:22 PM UTC-7, Daniel Bliss wrote:
>
> Dan,
>
> Many thanks for the reply. Right now I have the following:
>
> |

> simulation_clock =brian2.Clock(dt=0.02*brian2.__ms
> <http://brian2.ms>)
> fext =1800.0*brian2.hertz
> pe_pfc =brian2.NeuronGroup(NE,...,__clock=simulation_clock)
> pg_e_pfc =brian2.PoissonGroup(NE,fext,__clock=simulation_clock)
> cpe_pfc =brian2.Synapses(pg_e_pfc,pe___pfc,...,connect='i ==

> http://brian2.readthedocs.org/__en/latest/reference/brian2.__input.poissoninput.__PoissonInput.html
> <http://brian2.readthedocs.org/en/latest/reference/brian2.input.poissoninput.PoissonInput.html>
>
> <http://brian2.readthedocs.__org/en/latest/reference/__brian2.input.poissoninput.__PoissonInput.html

> > an email to brian-developm...@__googlegroups.com
> <mailto:brian-developm...@googlegroups.com>
> >
> <mailto:brian-developme...@googlegroups.com
> <mailto:brian-development%2Bunsu...@googlegroups.com>>.
> > For more options, visit
> https://groups.google.com/d/__optout
> <https://groups.google.com/d/optout>
> <https://groups.google.com/d/__optout

> <https://groups.google.com/d/optout>>.
>
> --
>
> ---
> You received this message because you are subscribed to the Google
> Groups "Brian Development" group.
> To unsubscribe from this group and stop receiving emails from
> it, send

> an email to brian-development+unsubscribe@__googlegroups.com
> <mailto:brian-development%2Bunsu...@googlegroups.com>
> <mailto:brian-developme...@googlegroups.com
> <mailto:brian-development%2Bunsu...@googlegroups.com>>.
> For more options, visit https://groups.google.com/d/__optout

> <https://groups.google.com/d/optout>.
>
>
> --
>
> --- You received this message because you are subscribed to the
> Google Groups "Brian Development" group.
> To unsubscribe from this group and stop receiving emails from it,

> send an email to brian-development+unsubscribe@__googlegroups.com
> <mailto:brian-development%2Bunsu...@googlegroups.com>.
> For more options, visit https://groups.google.com/d/__optout

> <https://groups.google.com/d/optout>.
>
>
> --
>
> ---
> You received this message because you are subscribed to the Google
> Groups "Brian Development" group.
> To unsubscribe from this group and stop receiving emails from it, send

> an email to brian-developm...@googlegroups.com
> <mailto:brian-developm...@googlegroups.com>.

[Daniel Bliss]

Hi Dan,

I've finally gotten around to making the switch from PoissonGroup to PoissonInput.  Unfortunately, after making this switch, I've seen drastic changes in my results that imply I must have done something wrong.  Here is a git diff that shows my PoissonGroup (the old) and PoissonInput (the new) implementations:

-fext = 1800.0 * brian2.hertz
+n_ext = 1000
+f_ext = 1.8 * brian2.hertz
 
 sim_dt = 0.02 * brian2.ms
 
@@ -152,11 +153,6 @@ if __name__ == '__main__':
                                 clock=simulation_clock, method='rk2',
                                 namespace=name_space['mt']['inh'])
 
-    pg_e_pfc = brian2.PoissonGroup(NE, fext, clock=simulation_clock)
-    pg_i_pfc = brian2.PoissonGroup(NI, fext, clock=simulation_clock)
-    pg_e_mt = brian2.PoissonGroup(NE, fext, clock=simulation_clock)
-    pg_i_mt = brian2.PoissonGroup(NI, fext, clock=simulation_clock)
-
     # -------------------------------------------------------------------------
     # AMPA Connections
     # -------------------------------------------------------------------------
@@ -214,15 +210,10 @@ if __name__ == '__main__':
     # External Connections (from Unmodeled Cortical Areas)
     # -------------------------------------------------------------------------
 
-    cpe_pfc = brian2.Synapses(pg_e_pfc, pe_pfc, pre='s_ext += 1',
-                              connect='i == j', clock=simulation_clock)
-    cpi_pfc = brian2.Synapses(pg_i_pfc, pi_pfc, pre='s_ext += 1',
-                              connect='i == j', clock=simulation_clock)
-
-    cpe_mt = brian2.Synapses(pg_e_mt, pe_mt, pre='s_ext += 1',
-                             connect='i == j', clock=simulation_clock)
-    cpi_mt = brian2.Synapses(pg_i_mt, pi_mt, pre='s_ext += 1',
-                             connect='i == j', clock=simulation_clock)
+    pg_e_pfc = brian2.PoissonInput(pe_pfc, 's_ext', n_ext, f_ext, weight=1)
+    pg_i_pfc = brian2.PoissonInput(pi_pfc, 's_ext', n_ext, f_ext, weight=1)
+    pg_e_mt = brian2.PoissonInput(pe_mt, 's_ext', n_ext, f_ext, weight=1)
+    pg_i_mt = brian2.PoissonInput(pi_mt, 's_ext', n_ext, f_ext, weight=1)

My first guess is that it was a mistake to set weight to an int.  Perhaps it should be '1' (a string)?  (Another difference between the old and new implementations is that the new one uses 1000 units at 1.8 Hz, rather than 1 unit at 1800 Hz -- though I believe these should be equivalent.)

(As a side note, I read the Fourcaud & Brunel paper you recommended.  It's not clear to me whether they provide a definition for Poisson input that could help me here.)

Many thanks, as always,

Dan B

To unsubscribe from this group and stop receiving emails from it, send an email to brian-developm...@googlegroups.com.

[Daniel Bliss]

Update: The PoissonInput code in my previous email is giving me an effective firing rate of 360 Hz (the same for all neurons).  I'm not sure whether this has to do with the fact that PoissonInput is interpreting my weight as having units of radians.  (This passes the units check even though my variable 's_ext' has units of '1' -- i.e., no units.)

-- Dan B

[Daniel Bliss]

I think the problem may be that PoissonInput is not passing CodeRunner a clock or dt, so it's using the default clock, not my simulation clock.

I'll make a fix, see whether it works, and then submit a pull request.

[Marcel Stimberg]

Hi Dan,

> I think the problem may be that PoissonInput is not passing CodeRunner
> a clock or dt, so it's using the default clock, not my simulation clock.

well spotted, PoissonInput should use the clock of the group it connects
to, but apparently it does not.

Note that if your "simulation clock" is the same for all objects, then
you should rather not create a new Clock object with the dt you want,
but instead set defaultclock.dt. Even if there are a few objects for
which you'd like to use a different dt, use their dt keyword argument
instead of creating a new Clock (this is more of a stylistic issue,
behind the scenes it should do the same thing).

Best,
Marcel

[Daniel Bliss]

Many thanks for the style tips.  I've adopted them, and, as you say, the resultant code functions in the same way as the code it replaced.

(Pull request for PoissonInput submitted, as you know.)

[Daniel Bliss]

Hi Marcel,

I'm digging up an old email thread between us here (see below).  Is setting defaultclock.dt at the top of my script still the best way to go?  For some reason I'm getting errors in a new script I've just written, saying that my dt is changing after defining a PoissonInput.  This seems to have something to do with PoissonInput's before_run method.  I'm looking into the cause of the error now, but it'd be great to hear your thoughts.  I'll report back if I find a solution.

Thanks,

Dan

On Wed, May 6, 2015 at 8:43 AM, Marcel Stimberg <marcel....@inserm.fr> wrote:

[Daniel Bliss]

P.S. Setting dt=brian2.defaultclock.dt in all my calls to NeuronGroup (after defining brian2.defaultclock.dt at the top of my script) seems to solve the problem.  Still, seems inelegant.  Am I doing something wrong?

-- Dan

[Marcel Stimberg]

Hi Daniel,

could you give some detail on how you create the relevant objects? You first set the dt of the defaultclock and then you create a NeuronGroup and a PoissonInput group without any dt/clock information set? When you get an error about a changing dt then this should mean that the dt of the group that PoissonInput connects to changed between the creation of the PoissonInput and the run call.

Best,
  Marcel

[Daniel Bliss]

Sorry, I missed your reply somehow.  Here's my code:

brian2.defaultclock.dt = 0.02 * brian2.ms

# . . . 

pe = brian2.NeuronGroup(NE, eqs_e, threshold='v > Vth_e', reset=reset_e,

                            refractory='tr_e', method='rk2',

                            namespace=name_space['exc'],

                            dt=brian2.defaultclock.dt)

pi = brian2.NeuronGroup(NI, eqs_i, threshold='v > Vth_i', reset=reset_i,

                            refractory='tr_i', method='rk2',

                            namespace=name_space['inh'],

                            dt=brian2.defaultclock.dt)

# . . . 

pg_e = brian2.PoissonInput(pe, 's_ext', n_ext, f_ext, weight=1)

pg_i = brian2.PoissonInput(pi, 's_ext', n_ext, f_ext, weight=1)

When I remove the dt argumen to my NeuronGroup calls, I get the error.  Let me know if more code would be helpful here.

Best,

Dan

[Marcel Stimberg]

Hi Dan,

thanks, I could reproduce your issue, it is an actual bug and I filed an
issue for it here: https://github.com/brian-team/brian2/issues/684

It should not be very difficult to fix, but until then please continue
to use your workaround.

Thanks for reporting this, best
Marcel