Problem with explicitly assigning current variables in model

Description of problem

hi there, I am new to Brian2 so apologies if this is an obvious mistake.

I am trying to run an LIF networks where I can explicitly add current to two variables I_e (external current) and I_i (internal current) which are then added onto the voltage in the model at each step.

As a test to see that this is working correctly I have simply added the inputs from I_e and I_i to the current voltage, using the equations provided from Brian2 LIF examples and compared the spikes. But weirdly when I do this the results look quite different - I was wondering if I am incorrectly assigning my variables?

Minimal code to reproduce problem

here is my model:
lif ="""
dv/dt = -(V-v_rest) / t_tau : 1 (unless refractory)
V = v + I_e + I_i : 1 #combine current voltage with incoming currents
I_i : 1 #internal current
I_e : 1 #external current"""

here is the model i am trying to reproduce - if my equation is doing what it should be doing the output should be the same
lif = “”"
dv/dt = -(v-v_rest) / t_tau : 1 (unless refractory)"""

What you have already tried

i think the error may be to do with how I am adding current to my I_i and I_e variables.

In the online examples i can see that internal weights are accumulated to the voltage variable, via on_pre as ‘v+=w’.
Given that I am adding I_i onto v in my equation I have changed this to on_pre: ‘I_i = w’

Similarly for the Poisson input group the target_var from the online examples was ‘target_var = v’, which I have changed to ‘target_var = I_e’

One thing in particular I have noticed is that it seems like the I_e variable accumlates current from previous time steps when it should just add the current from the input from that time step alone - I am not sure how exactly to fix this?

Any thoughts on why the equations are not doing the same thing would be much appreciated!

Thanks!

Hi. I think there is some confusion around the various types of synaptic interactions/inputs. You are certainly not the first or only user with this problem, we should definitely address this more concretely in some kind of tutorial/documentation/blog post… For your specific question, I am having a bit of a trouble understanding the details, a runnable example would be great!

The equations in lif look ok in principal, but maybe a bit uncommon in the way you define V. A more common way would be to write:

lif = '''
dv/dt = -((v - v_rest) - I_e - I_i ) / tau : 1 (unless refractory)
I_i : 1
I_e : 1
'''

(The parentheses around (v - v_rest) are of course not strictly necessary, but it makes it clearer that your equation has three currents: the leak current given by (v - v_rest) and the currents I_i and I_e). As a really minor note: I think most people would think that I_e and I_i stand for excitatory and inhibitory current. If this is not the case for you (i.e. these are really “internal” and “external” currents), then I’d probably change their names a bit.

Now, I think the confusion stems from the fact that I_i and I_e are currents that only change when you assign to them, i.e. they don’t have any dynamics in the absence of input. If you used a summed variable, this would be ok, since this mechanism sets the target variable at each time step. If you use the on_pre statement of a Synapses object, however, this only triggers a change in the target variable for each incoming spike. Setting on_pre='I_i = w' is therefore almost certainly not what you want to do, since it will fix the value of the I_i current for eternity. Similarly, a PoissonInput models the additive effect of a large number of impinging spikes for every timestep, and is equivalent to Poisson-distributed spikes followed by on_pre='I_e += w'. So in this case it does not “clamp” the current to a fixed value, but instead continues to increase it (as you observed).

In the examples you cite where v += w or target_var = 'v' is used, the synaptic effect is modelled as a so-called “delta synapse”, i.e. an instantaneous jump in the membrane potential, without any explicitly modelled current. Each of the incoming spikes (or summed activity in the case of PoissonInput) makes the membrane potential go up, but its leak current then makes the membrane potential go back to zero again.

If you want to actually model a current, then you need to specify its dynamics. The typical way of doing it is to say that it decays exponentially back to 0 in absence of input and jumps up with each spike. I.e. in your model:

lif = '''...
dI_i/dt = -I_i/tau_i : 1
dI_e/dt = -I_e/tau_e : 1
...'''
synapses = Synapses(..., on_pre='I_i += w')
poisson_input = PoissonInput(..., target_var='I_e')

The time constants tau_e and tau_i will determine the duration of the effect of each incoming spike. Hope that makes things a bit clearer?

1 Like

Hi Marcel - thanks for the prompt response!

I think I have been using incorrect terminology so apologies for not being clear - what I am looking to do is model my synaptic effects as ‘delta synapses’ as in the example, but instead of having the membrane potential be updated in the v variable, have it update two separate variables that i can then divide by a leak term, before adding to v.

I_int: At each time step, for each neuron add the summed weights (as defined in my synapses group) of all incoming spikes onto the membrane potential (see below).

i_ext: At each time step for each neuron add the Poisson input onto the membrane potential

These inputs could then be divided by a leak term (g_L) to get the leak current which is added onto the voltage at each time step. When I run this my dynamics are drastically different from what I see when I run the LIF example (’-(v-v_rest)/tau’ with ‘on_pre: v+=w’ and ‘target_var = v’) and that’s where my confusion has arisen from.

Judging from your response, the problem may be arising from the fact that values are only updated when a spike occurs, and then is fixed until another spike occurs? So in instances when there are 0 spikes the current would still carry over from before. In that case the answer may be setting I_int and i_ext as summed variables instead so they are updated at each time step rather than only when a spike occurs?

Please see my full code below:

#Network dynamics
N = nodes.shape[0] #Define number of cells
v_rest= 0 #Resting potential
v_th = 1 #Spike threshold
t_ref= 2.0 * b2.ms #absolute refractory period
t_syn_del = 1.5 * b2.ms #delay between presynaptic spike and postsynaptic increase
t_tau = 20. * b2.ms #time taken for voltage to reach 63% of final value
g_L = 1 #leak current

#define dynamics for each cell
lif = ‘’’
dv/dt = -(v - v_rest) + ((I_ext + I_int)/g_L) / t_tau : 1 (unless refractory)
I_int : 1
I_ext : 1
‘’’

net_dyn = b2.NeuronGroup(
N, model=lif,
threshold=“v>v_th”, reset=“v = v_rest”, refractory=t_ref,
method=“euler”)

net_dyn.v = v_rest #set starting value for voltage

#Network connectivity + weights - ignore the details of this as this is a separate class I am using to define network connectivity
s = i
k = 15
p = 0.3
curr = bap_netsim(dist).adjmat_generate(s, p, k, divisor, soften)
A = curr.A
W = curr.adj_mat

#Build synapses
net_syn = b2.Synapses(net_dyn, net_dyn, ‘w:1’, on_pre=“I_int=w”, delay=t_syn_del)
rows, cols = np.nonzero(A)
net_syn.connect(i = rows, j = cols)
net_syn.w = W[rows, cols]

#External input
N_extern=400 #number of presynaptic excitatory poisson neurons
poisson_input_rate= 10. * b2.Hz #poisson rate of external population
w_external= 0.01 #synaptic weight of excitatory external possion neurons onto all neurons
external_poisson_input = b2.PoissonInput(target=net_dyn, target_var=“I_ext”, N=N_extern,
rate=poisson_input_rate, weight=w_external)

On a side note, the reason i am explicitly defining the internal and external values rather than using the simple ‘-(v-v_rest)/tau’ example is that the formulae I am familiar with for leaky integrators explicitly separates out the incoming current and divides by a leak term - I am not sure why/how we can do away with that in the examples you provide.

Thanks so much for the help!

Hi again. I have to admit I am now more confused than before. If you want to use “delta synapses”, then they have to update the membrane potential. Can you say why you want to have “two separate variables”? Do you want to keep track of updates from different sources? If you don’t need that, you can have everything target v directly.
Regarding the division by the leak conductance (which is 1 in your example :slight_smile: ). Since you are using a unitless formulation of the integrate-and-fire model, this is normally already considered to be taken into account when defining the strength of the increase. In general, I’d say that it is unnecessary to link back the increase of the membrane potential to a current divided by a conductance if you are using such a simple delta synapse model. It’s an infinitely short current pulse, but its effect is really straightforward: it increases the membrane potential by a given amount. So if you don’t have some kind of experimental measurement that makes it necessary to link your model values back to the experimental values, I would not complicate things further. If this is your situation, then probably you do not want to use this kind of neuron/synapse model in the first place.

All that said: if you need to separate out the contributions from the two groups for analysis purposes, then there is a solution to do things similar to your current approach, but it is a bit hacky and non-standard. I can show you how to do this if necessary.

This makes sense - I wasn’t totally sure why we could do away with explicitly defining the incoming current + leak. That was the only reason I thought I needed to individually define them, so it sounds like sticking with simple delta synapses that directly update v is all I need.

Thanks!

1 Like