Problem with biexponential synaptic currents: no currents recorded

[Hedyeh]

Dear all,
I have problem with implementing biexponential synaptic currents. Actually, no synaptic currents has been recorded. In case I use s_GABA_rec_tot_post, the error I get is: s_GABA_rec_tot_post cannot be used as a variable name, the -pre and -post suffixes are used to refer to …
While when I changed the name to s_GABA_rec_tot_in the error did not came up, but no currents was recorded. I wonder what causes such a problem.
My other question is whether ,summed, should be there? And what is the difference between AMPA and NMDA currents from modeling point of view?
A part of code regarding the equations is attached in the following.

Thank you in advance for your help and time.
Hedyeh

eqs_E = ‘’’
dv / dt = (- g_m_E * (v - V_L) - I_syn) / C_m_E : volt (unless refractory)

I_syn = I_AMPA_ext + I_AMPA_rec + I_NMDA_rec + I_GABA_rec : amp

I_AMPA_ext = g_AMPA_ext_E * (v - V_E) * s_AMPA_ext : amp
ds_AMPA_ext / dt = (- s_AMPA_ext / tau_AMPA_decay) : 1

I_NMDA_rec = g_NMDA_E * (v - V_E) / (1 + Mg2 * exp(-0.062 * v / mV) / 3.57) * s_NMDA_tot : amp
s_NMDA_tot : 1

I_AMPA_rec = g_AMPA_rec_E * (v - V_E) * s_AMPA_rec_tot : amp
s_AMPA_rec_tot : 1

I_GABA_rec = g_GABA_E * (v - V_I) * s_GABA_rec_tot : amp
s_GABA_rec_tot : 1
‘’’

eqs_I = ‘’’
dv / dt = (- g_m_I * (v - V_L) - I_syn) / C_m_I : volt (unless refractory)

I_syn = I_AMPA_ext + I_AMPA_rec + I_NMDA_rec + I_GABA_rec : amp

I_AMPA_ext = g_AMPA_ext_I * (v - V_E) * s_AMPA_ext : amp
ds_AMPA_ext / dt = (- s_AMPA_ext / tau_AMPA_decay) : 1

I_NMDA_rec = g_NMDA_I * (v - V_E) / (1 + Mg2 * exp(-0.062 * v / mV) / 3.57) * s_NMDA_tot : amp
s_NMDA_tot : 1

I_AMPA_rec = g_AMPA_rec_I * (v - V_E) * s_AMPA_rec_tot : amp
s_AMPA_rec_tot : 1

I_GABA_rec = g_GABA_I * (v - V_I) * s_GABA_rec_tot : amp
s_GABA_rec_tot : 1
‘’’

populations

E = NeuronGroup(N_E, eqs_E, threshold=‘v > V_thr’, reset=‘v = V_reset’, refractory=tau_rp_E, method=‘euler’)
E.v = V_L
I = NeuronGroup(N_I, eqs_I, threshold=‘v > V_thr’, reset=‘v = V_reset’, refractory=tau_rp_I, method=‘euler’)
I.v = V_L

Synaptic equations

eqs_NMDA = ‘’’
s_NMDA_tot_post = w * s_NMDA : 1 (summed)
ds_NMDA / dt = - s_NMDA / tau_NMDA_decay + alpha * x * (1 - s_NMDA) : 1 (clock-driven)
dx / dt = - x / tau_NMDA_rise : 1 (clock-driven)
w : 1
‘’’

eqs_AMPA_rec=‘’’
s_AMPA_rec_tot_post = w * s_AMPA_rec : 1 (summed)

ds_AMPA_rec / dt = ((tau_AMPA_decay / tau_AMPA_rise) ** (tau_AMPA_rise / (tau_AMPA_decay - tau_AMPA_rise)) * x_AMPA_rec - s_AMPA_rec ) / tau_AMPA_rise : 1 (clock-driven)

dx_AMPA_rec / dt = (- x_AMPA_rec / tau_AMPA_decay) : 1 (clock-driven)
w : 1
‘’’

eqs_GABA_rec=‘’’
s_GABA_rec_tot_post = w * s_GABA : 1
ds_GABA / dt = ((tau_GABA_decay / tau_GABA_rise) ** (tau_GABA_rise / (tau_GABA_decay - tau_GABA_rise)) * x_GABA - s_GABA ) / tau_GABA_rise : 1 (clock-driven)

dx_GABA/ dt = (- x_GABA / tau_GABA_decay) : 1 (clock-driven)
w : 1
‘’’

Receptors

eqs_pre_NMDA = ‘’’
x += 1
‘’’

eqs_pre_rec = ‘’’
x_AMPA_rec += 1
‘’’

eqs_pre_ext = ‘’’
s_AMPA_ext += 1
‘’’

eqs_pre_gaba = ‘’’
x_GABA += 1
‘’’

Connections involving NMDA synapses

C_E_E = Synapses(E, E, model=eqs_NMDA, on_pre=eqs_pre_NMDA, delay = 1.0 * ms, method=‘euler’)
C_E_E.connect(p=sparsness)
C_E_E.connect(‘i != j’)

C_E_I = Synapses(E, I, model=eqs_NMDA, on_pre=eqs_pre_NMDA, delay = 1.0 * ms, method=‘euler’)
C_E_I.connect(p=sparsness)

Connections involving AMPA synapses

C_E_E = Synapses(E, E, model=eqs_AMPA_rec, on_pre=eqs_pre_rec, delay = 1.0 * ms, method=‘euler’)
C_E_E.connect(p=sparsness)
C_E_E.connect(‘i != j’)

C_E_I = Synapses(E, I, model=eqs_AMPA_rec, on_pre=eqs_pre_rec, delay = 1.0 * ms, method=‘euler’)
C_E_I.connect(p=sparsness)

Connections involving GABA synapses

C_I_E = Synapses(I, E, model=eqs_GABA_rec, on_pre=eqs_pre_gaba, delay = 1.0 * ms, method=‘euler’)
C_I_E.connect(p=sparsness)

C_I_I = Synapses(I, I, model=eqs_GABA_rec, on_pre=eqs_pre_gaba, delay = 1.0 * ms, method=‘euler’)
C_I_I.connect(p=sparsness)
C_I_I.connect(‘i != j’)

External inputs

C_P_E = PoissonInput(E, ‘s_AMPA_ext’, N=N_X, rate=rate_exc, weight=1.52)
C_P_I = PoissonInput(I, ‘s_AMPA_ext’, N=N_X, rate=rate_inh, weight=1.52)

[mstimberg]

Hi @Hedyeh,

@Hedyeh

I have problem with implementing biexponential synaptic currents. Actually, no synaptic currents has been recorded. In case I use s_GABA_rec_tot_post, the error I get is: s_GABA_rec_tot_post cannot be used as a variable name, the -pre and -post suffixes are used to refer to …
While when I changed the name to s_GABA_rec_tot_in the error did not came up, but no currents was recorded. I wonder what causes such a problem.
My other question is whether ,summed, should be there?

An expression like s_GABA_rec_tot_post = w * s_GABA : 1 (summed) in the synapses means: at every time step, set the s_GABA_rec_tot variable of each post-synaptic cell to the sum of w * s_GABA, summing over all synapses that connect to the cell. You can record the value of s_GABA for each synapse with something like:

synapse_mon = StateMonitor(C_I_E, 's_GABA', record=True)

or you can record the total value of s_GABA for each neuron with:

neuron_mon  = StateMonitor(E, 's_GABA_rec_tot', record=True)

or the same for the total current (... = StateMonitor(E, 'I_GABA_rec', record=True)).

You cannot write s_GABA_rec_tot_post = ... without (summed) in the synaptic equations, since there is more than one synapse targetting each post-synaptic cell. Each of these synapses would have a different right-hand-side of the expression, it would therefore be undefined what the value for the post-synaptic cell should be.

If you use a name that does not end in post (e.g. s_GABA_rec_tot_in), then you cannot use summed. What you are doing in this case is defining an expression for internal use in the synapse. You can record this expression, but it does not have any effect on the post-synaptic cell, it is therefore not of much use.

@Hedyeh

And what is the difference between AMPA and NMDA currents from modeling point of view?

The equations that you wrote for AMPA and GABA are correct, but very inefficient. Since the equations use the same time constants for all synapses, you can have a much more efficient implementation. What you are currently calculating is transforming the incoming spikes into a biexponential current for each synapse, and then you add up all these currents for each target neuron. Mathematically, this is equivalent (since the differential equations are linear, and the time constants are the same) to transforming the incoming spikes directly into a current of the target neuron. I.e., instead of

# synapse equations
eqs_AMPA_rec='''
s_AMPA_rec_tot_post = w * s_AMPA_rec : 1 (summed)

ds_AMPA_rec / dt = ((tau_AMPA_decay / tau_AMPA_rise) ** (tau_AMPA_rise / (tau_AMPA_decay - tau_AMPA_rise)) * x_AMPA_rec - s_AMPA_rec ) / tau_AMPA_rise : 1 (clock-driven)

dx_AMPA_rec / dt = (- x_AMPA_rec / tau_AMPA_decay) : 1 (clock-driven)
w : 1
'''
eqs_pre_rec = '''
x_AMPA_rec += 1
'''
# neuron equations 
eqs = '''...
I_AMPA_rec = g_AMPA_rec_E * (v - V_E) * s_AMPA_rec_tot : amp
s_AMPA_rec_tot : 1
'''

you’d use:

# synapse equations
eqs_AMPA_rec='w : 1'
eqs_pre_rec = '''
x_AMPA_rec += w
'''
# neuron equations 
eqs = '''...
I_AMPA_rec = g_AMPA_rec_E * (v - V_E) * s_AMPA_rec : amp
ds_AMPA_rec / dt = ((tau_AMPA_decay / tau_AMPA_rise) ** (tau_AMPA_rise / (tau_AMPA_decay - tau_AMPA_rise)) * x_AMPA_rec - s_AMPA_rec ) / tau_AMPA_rise : 1

dx_AMPA_rec / dt = (- x_AMPA_rec / tau_AMPA_decay) : 1
'''

This is much more efficient, since it does the integration of the differential equations only once per neuron and not once per synapse.

However, you cannot do this for the NMDA equations, since they are non-linear functions of v. You’d therefore actually have to do the opposite, i.e. move the definition of the current into the synapse. I think it should be something like this:

# synapse equations
eqs_NMDA = '''
I_NMDA_rec_post = g_NMDA_I * (v_post - V_E) / (1 + Mg2 * exp(-0.062 * v_post / mV) / 3.57) * (w * s_NMDA) : amp (summed)
ds_NMDA / dt = - s_NMDA / tau_NMDA_decay + alpha * x * (1 - s_NMDA) : 1 (clock-driven)
dx / dt = - x / tau_NMDA_rise : 1 (clock-driven)
w : 1
'''
eqs_pre_NMDA = '''
x += 1
'''
# neuron equations 
eqs = '''...
I_NMDA_rec : amp
'''

Hope that answers your questions. If you are interested why NMDA is modelled differently from AMPA/GABA, have a look here: 3.1 Synapses | Neuronal Dynamics online book

PS:
You can wrap your code like this to be displayed nicely:

```python
print('my code')
```

[Hedyeh]

Hi Marcel,

Many thanks for your detailed reply. Very useful and informative, indeed!
I still have questions: By recording synaptic currents, in fact the incoming currents to each neuron is recorded, right?
I tried to plot s_AMPA_rec as well as synaptic currents, but they all stick to zero with no change. However, neurons are spiking due to external poissonian currents. Why should it be so?
you said I am using the same time constants, but that is not the case as both synaptic and membrane time constants are different across neurons and AMPA and GABA and NMDA synapses. Am I mistaken?
Thanks a lot.
Best,
Hedyeh

[mstimberg]

Hi,

@Hedyeh

I still have questions: By recording synaptic currents, in fact the incoming currents to each neuron is recorded, right?

Yes, there is one “total current” for each type of synapse and each neuron.

@Hedyeh

I tried to plot s_AMPA_rec as well as synaptic currents, but they all stick to zero with no change.

For s_AMPA_rec (recorded for each synapse in this case, right?), I do not see how stay at zero if the synapse receives spikes. For the current, it could stay at zero if the conductance (g_AMPA_rec) is zero, and/or when v is exactly equal to V_E, but I guess this is not the case… Are you sure that the neurons are actually spiking? They get Poissonian input, but maybe the input is too weak to make them spike?

@Hedyeh

you said I am using the same time constants, but that is not the case as both synaptic and membrane time constants are different across neurons and AMPA and GABA and NMDA synapses.

In the model you posted, the equations refer to external constants such as tau_AMPA_decay, etc., so these time constants are the same across all neurons/synapses. Or are you using different equations in your actual model?

[mstimberg]

@mstimberg

Are you sure that the neurons are actually spiking? They get Poissonian input, but maybe the input is too weak to make them spike?

Oh, or you could have forgot to set the weight of the synapses, which means that they are zero.

BTW, I just saw that you are connecting with two connect calls:

C_E_E.connect(p=sparsness)
C_E_E.connect('i != j')

I doubt that this is what you really want. This means to create synapses between all possible pairs with the given probability, and additionally to connect all neurons to all neurons (except to themselves). You probably want to use

C_E_E.connect('i != j', p=sparsness)

instead.

[Hedyeh]

Thanks again, I double checked:
v is not equal to V_E, neurons are spiking, I merged the two connection calls into one call containing both self-connection exclusion and sparsness; but as you can see from the plot, nothing changed, currents are still zero.
Regarding the time constants, if you mean considering a distribution, then that is not the case for me. You are right,

I use a specific time constant for each synaptic current.

[Hedyeh]

BTW, connection strengths are not zero…

[mstimberg]

@Hedyeh

I use a specific time constant for each synaptic current.

Right, what I meant is that all synapses of the same type use the same time constants.

@Hedyeh

Thanks again, I double checked:
v is not equal to V_E, neurons are spiking, I merged the two connection calls into one call containing both self-connection exclusion and sparsness; but as you can see from the plot, nothing changed, currents are still zero.

Did you also make sure that the synaptic weights w are not zero?

[Hedyeh]

Sorry for such stupid question: are the weights defined only inside neuron and synapse equations? or somewhere else as well?

[mstimberg]

Your synaptic equations define w : 1, which means that each synapse has a value for w. You’ll therefore have to set it after creating the synapses:

C_E_E.connect(...)
C_E_E.w = 'rand()'  # set a random weight

You can also set the same weight for all synapses with something like

C_E_E.w = 1

but in that case you don’t really need a weight per synapse, i.e. you could also remove w : 1 from the equations, and instead of w refer to something like w_e_e, a single value that you set outside of your equations (i.e. like tau_AMPA_rise and other constants).

[Hedyeh]

Dear Marcel,

Thanks to the points you advised, currents and s_AMPA_rec etc. are recorded successfully! Not zero anymore.
Thanks for creating such an environment so that everyone can raise their questions.

Regards,
Hedyeh