Hi folks, I am experimenting with implementing the non-Markovian Gillespie Algorithm (nMGA) in Brian based on Vestergaard and Genois (PLoS Comput. Biol. 2015, p9). Is there a way to retrieve real-time evaluations of the threshold/reset conditions?
At the moment my method is implemented by reaction=NeuronGroup(1,code=nMGA,...)
where
nMGA = '''
tau : 1
Q : 1
transition = tau<Q*dt : boolean
...
...
'''
Say I have a reaction in the nMGA scheme every time tau<Q*dt
where tau
and Q
are set inside my nMGA implementation (within reaction=NeuronGroup(...)
. This reaction is detected passing threshold='tau<Q*dt
to the NeuronGroup
.
Because the nMGA requires some tricks for multiple transitions within the same time step, I am using reaction.run_regularly(code=code_after_transition)
as well, where
code_after_transition = '''
if not transition:
...
else
...
'''
I want to avoid defining transition
in the reaction
equations and have it instead directly passed from the reaction.thresholder
or whatever evaluates the threshold condition. Something in the spirit of lastspike
.
Hi Maurizio,
we don’t have a variable that is explicitly set by the thresholder, spikes are stored in a dedicated data structure. I have to admit that I do not exactly understand the issue: why do you not want to define transtion
in the equations? If you wanted to set a variable after each spike (or transition), you could set it in the reset
. Or the same for a custom event and run_on_event
.
@mstimberg You are totally right. Sorry. I was caught in a classic loop-of-thoughts. To answer your question. The idea is that I have a transition = tau<Q*dt : boolean
in the equations. And then if (not transition): <do some stuff> else: <do some other stuff upon transition>
. This is the case of the nMGA, with the complication, that upon a transition, the algorithm requires to do some additional stuff after the reset and before updating to the next time step.
On a side note, following up on to your reply, do the reset
or threshold
code strings allow specification of complex code in the form of if (...) : else : ...
or for
or while
loops (the latter for example in the reset
code)? Or do I need to resort run_on_event
in these cases?
You cannot use control structures like if
/else
or for
/while
in Brian code, the only direct way to implement it would be in C++/Cython code as part of user-defined functions. If you want to assign two different values to a variable depending on a condition, you can use the workaround var = int(condition)*value_1 + int(not condition)*value_2
. Otherwise you’ll need to do everything with custom events and run_on_event
, e.g. you can do stuff like:
group = NeuronGroup(..., events={'transitioned': 'transition',
'no_transition': 'not transition'})
group.run_on_event('transitioned', 'do something')
group.run_on_event('no_transition', 'do something else')
1 Like