Problem with the ‘+=’ function

Support
1

I have constructed a network of excitatory pyramidal (PYR) neurons and inhibitory (PVN) neurons. For the PYR neurons, there is an adaptation conductance ‘u’ with a -65mV reversal potential that has an exponential decay and to which a fixed value is added at each spike threshold and voltage reset. However, in simulations ‘u’ is not receiving the increment upon voltage reset. Is there someone out there more experienced than me at spotting my coding error? See my code below.

Thanks very much, Mitch

#this is only some of the code, but should be enough. I monitored voltage and adaptive conductance ‘u’ in neuron 1. ‘u’ started at 1 nS, as specified in code and decayed exponentially. But when neuron fired action potential and reset, ‘u’ did not receive an increment.

### Define the PYR cell equations

eqs_pyr = “”"
Iext : amp
Isyn_e_pyr = g_e*(v-Ee)s_sum_e_pyr : amp
Isyn_i_pyr = g_i(v-Ei)s_sum_i_pyr : amp
s_sum_e_pyr : 1
s_sum_i_pyr : 1
Iapp : amp
k = (int(v<vt) * klow) + (int(v>=vt) * khigh) : siemens/volt
du/dt = a(-u) : siemens
dv/dt = (k*(v-vr)(v-vt)+Ishift+Iext-Isyn_e_pyr-Isyn_i_pyr-u(v-va))/C_pyr : volt
“”"
reset = ‘ ‘ ‘
v = c

u += 1
‘ ‘ ‘

### Define the PYR neuron group

PYR = NeuronGroup(N1, model=eqs_pyr, reset= reset, threshold=‘v >= vpeak’, method = ‘euler’)
PYR.Iapp = 0.00000000000015 * amp
PYR.Iext = (randn(len(PYR))PYR.Iapp10 + PYR.Iapp)
print(PYR)

### set initial conditions
PYR.v = (rand(len(PYR))*0.002 * volt) - (0.052 * volt)
PYR.u = 1 * nS

(there were no relevant error messages)

2

Hi @MitchG_HunterCollege, could you please provide a “minimum working example”, i.e. code that can be run ? I don’t see any reason why an increase in u should not be applied in the reset, but your code above does not (should not) run – it should raise a unit mismatch error: u += 1 should not work, since u has units of siemens.

3

Hello, Marcel
Thanks for helping.
Here is complete code. There are PYR and PVN groups and four synaptic models. Running the code successfully executes, but gives a warning related to the problematic variable ‘u’:
“WARNING ‘u’ is an internal variable …. But also exists in the run namespace with the value 1. * nsiemens. The internal variable will be used:
As you can see from the plots, after a time delay caused by the initial value u=1, PYR[0] receives multiple EPSCs and IPSCs and fires action potentials, but u does not undergo +1 reset.
My guess is there is a simple solution that I am not seeing.
Mitch

import matplotlib.pyplot as plt
from brian2 import *
defaultclock.dt = 0.02*ms
seed(3)

from brian2.codegen import runtime

C_pyr = 115 * pF
C_pvn = 30 * pF
vr = -52.0 * mV
vpeak = 30.0 * mV
c = -53.0 * mV
N1 = 5000 # number of PYR neurons in network
N2 = 50 # number of PVN neurons in network
klow = 0.1 * nS/mV
khigh = 3.3  * nS/mV
a = 0.008 /ms
va = -65 * mV
vt = -49.0 * mV
Ishift= 0 * pA

# synaptic parameters
g_e = 0.285 * nS
g_i = 0.285 * nS
alpha = 2000 /second
Beta = 333.33 /second
s_inf =alpha/(alpha+Beta)
tau_s =1/(alpha+Beta)
Ee=-15 *mV
Ei=-68 *mV

eqs_pyr = """
Iext  : amp
Isyn_e_pyr = g_e*(v-Ee)*s_sum_e_pyr : amp
Isyn_i_pyr = g_i*(v-Ei)*s_sum_i_pyr : amp
s_sum_e_pyr : 1
s_sum_i_pyr : 1
Iapp : amp
k = (int(v<vt) * klow) + (int(v>=vt) * khigh) : siemens/volt
du/dt = a*(-u) : siemens
dv/dt = (k*(v-vr)*(v-vt)+Ishift+Iext-Isyn_e_pyr-Isyn_i_pyr-u*(v-va))/C_pyr : volt
"""

reset = '''
u += 1
v = c
'''

eqs_pvn = """
Isyn_e_pvn = g_e*(v-Ee)*s_sum_e_pvn : amp
Isyn_i_pvn = g_i*(v-Ei)*s_sum_i_pvn : amp
s_sum_e_pvn : 1
s_sum_i_pvn : 1
k = (int(v<vt) * klow) + (int(v>=vt) * khigh) : siemens/volt
dv/dt = (k*(v-vr)*(v-vt)+Ishift-Isyn_e_pvn-Isyn_i_pvn)/C_pvn : volt
"""

reset = '''
v = c
'''

PYR = NeuronGroup(N1, model=eqs_pyr, reset= reset, threshold='v >= vpeak', method = 'euler')
# create subgroups
PYR1 = PYR[:4950]
PYR2 = PYR[50:]
PYR1.Iapp = 0.00000000000015 * amp
PYR2.Iapp = 0.00000000000015 * amp
PYR1.Iext = (randn(len(PYR1))*PYR1.Iapp*10 + PYR1.Iapp)
PYR2.Iext = (randn(len(PYR2))*PYR2.Iapp*10 + PYR2.Iapp)
print(PYR)

PVN = NeuronGroup(N2, model=eqs_pvn, reset= reset, threshold='v >= vpeak', method = 'euler')
print(PVN)

S_pyr_pyr=Synapses(PYR,PYR,
       model="""
              ds0/dt=-(s0-s_inf)/(tau_s) :1 (clock-driven)
              ds1/dt=-s1*Beta :1 (clock-driven)
              s_tot=clip(s0,0,s1) :1 
              s_sum_e_pyr_post = s_tot : 1 (summed)
       """,
       on_pre="""
              s0=s1 
              s1=(s_inf*(1-exp(-0.001*second/tau_s))+s1*exp(-0.001*second/tau_s))*exp(0.001*Beta*second)
              delay=0.001*second
       """, 
          method = 'exact')
S_pvn_pyr=Synapses(PVN,PYR,
       model="""
              ds0/dt=-(s0-s_inf)/(tau_s) :1 (clock-driven)
              ds1/dt=-s1*Beta :1 (clock-driven)
              s_tot=clip(s0,0,s1) :1
              s_sum_i_pyr_post = s_tot : 1 (summed)
       """,
       on_pre="""
              s0=s1
              s1=(s_inf*(1-exp(-0.001*second/tau_s))+s1*exp(-0.001*second/tau_s))*exp(0.001*Beta*second)
              delay=0.001*second
       """,
          method = 'exact')

S_pyr_pvn=Synapses(PYR,PVN,
       model="""
              ds0/dt=-(s0-s_inf)/(tau_s) :1 (clock-driven)
              ds1/dt=-s1*Beta :1 (clock-driven)
              s_tot=clip(s0,0,s1) :1
              s_sum_e_pvn_post = s_tot : 1 (summed)
       """,
       on_pre="""
              s0=s1
              s1=(s_inf*(1-exp(-0.001*second/tau_s))+s1*exp(-0.001*second/tau_s))*exp(0.001*Beta*second)
              delay=0.001*second
       """,
          method = 'exact')

S_pvn_pvn=Synapses(PVN,PVN,
       model="""
              ds0/dt=-(s0-s_inf)/(tau_s) :1 (clock-driven)
              ds1/dt=-s1*Beta :1 (clock-driven)
              s_tot=clip(s0,0,s1) :1
              s_sum_i_pvn_post = s_tot : 1 (summed)
       """,
       on_pre="""
              s0=s1
              s1=(s_inf*(1-exp(-0.001*second/tau_s))+s1*exp(-0.001*second/tau_s))*exp(0.001*Beta*second)
              delay=0.001*second
       """,
          method = 'exact')
# set PYR to PYR connectivity to 0.2%
S_pyr_pyr.connect(condition='rand() < 0.002')

# set PYR to PVN connectivity to 5%
S_pyr_pvn.connect(condition='rand() < 0.05')

# set PVN to PVN connectivity to 5%
S_pvn_pvn.connect(condition='rand() < 0.05')

# set PVN to PYR connectivity to 5%
S_pvn_pyr.connect(condition='rand() < 0.05')

# set initial conditions
S_pyr_pyr.s0=0
S_pyr_pyr.s1=0
S_pyr_pvn.s0=0
S_pyr_pvn.s1=0
S_pvn_pyr.s0=0
S_pvn_pyr.s1=0
S_pvn_pvn.s0=0
S_pvn_pvn.s1=0
PYR.v = (rand(len(PYR))*0.002 * volt) - (0.052 * volt)
PVN.v = (rand(len(PVN))*0.002 * volt) - (0.052 * volt)
PYR.u = 1 * nS

#record PYR membrane potential, adaptive conductance, synaptic currents of cells and spike times in network
PYR_v = StateMonitor(PYR,'v',record = [0,1,2,3,4,995,996,997,998,999])
PYR_u = StateMonitor(PYR,'u',record = [0,1,2,3,4,995,996,997,998,999])
PYR_Isyn_e_pyr = StateMonitor(PYR,'Isyn_e_pyr',record= [4])
PYR_Isyn_i_pyr = StateMonitor(PYR,'Isyn_i_pyr',record= [4])
PYR_spktimes = SpikeMonitor(PYR, record=[0,1,2,3,4,995,996,997,998,999])
#record PVN membrane potential of cells and spike times in network
PVN_v = StateMonitor(PVN,'v',record = [0,1,2,3,4,5])
PVN_spktimes = SpikeMonitor(PVN, record=[0,1,2,3,4,5,6,7,8,9])
#run for x seconds of simulated time
duration = 1 * second

#include neuron groups, synapse groups, and monitors in our simulation
net = Network(PYR,PYR1,PYR2,PVN,S_pyr_pyr,S_pyr_pvn,S_pvn_pyr,S_pvn_pvn,PYR_v,PYR_spktimes,PYR_u,PYR_Isyn_e_pyr,PYR_Isyn_i_pyr,PVN_v,PVN_spktimes)
net.run(duration)

# plot PYR membrane potential
plt.figure(0)
plt.plot(PYR_v.t/ms,PYR_v.v[0]/mV)
plt.xlabel("Time (ms)")
plt.ylabel("PYR Membrane Potential (mV)")
plt.title('PYR network Cell 1 Vm')
plt.show()

# plot PYR Adaptive Conductance u
plt.figure(1)
plt.plot(PYR_u.t/ms,PYR_u.u[0]/nsiemens)
plt.xlabel("Time (ms)")
plt.ylabel("PYR Adaptive Conductance u (nS)")
plt.title('PYR network Cell 1 Adaptive Conductance')
plt.show()

# plot PYR excitatory synaptic current Isyn_e_pyr
plt.figure(2)
plt.plot(PYR_Isyn_e_pyr.t/ms,PYR_Isyn_e_pyr.Isyn_e_pyr[0]/pA)
plt.xlabel("Time (ms)")
plt.ylabel("EPSCs (pA)")
plt.title('PYR network Cell 1 EPSCs')
plt.show()

# plot PYR inhibitory synaptic current Isyn_i_pyr
plt.figure(3)
plt.plot(PYR_Isyn_i_pyr.t/ms,PYR_Isyn_i_pyr.Isyn_i_pyr[0]/pA)
plt.xlabel("Time (ms)")
plt.ylabel("IPSCs (pA)")
plt.title('PYR network Cell 1 IPSCs')
plt.show()
4

Hi @MitchG_HunterCollege – I wouldn’t necessarily call this a minimal example :wink: Also, please wrap code examples in triple backticks like this

```
print('some python code')
```

This will make it more readable and easy to copy & paste :pray:

Regarding your code, in the beginning you write:

eqs_pyr = """..."""

reset = '''
u += 1
v = c
'''

eqs_pvn = """..."""

reset = '''
v = c
'''

PYR = NeuronGroup(N1, model=eqs_pyr, reset= reset, threshold='v >= vpeak', method = 'euler')

The second assignment to reset overwrites the first one, before you hand it over to NeuronGroup, therefore both of your groups use the reset v = c, without updating u! When you fix this, you will get an error, though, since the reset for PYR should include something like u += 1*nS instead of u += 1.

Hope that gets you going!

5

Hi Marcel. Thank you for identifying the likely flaw in my code. I will work on this and let you know how it resolves. Regards, Mitch

6

Hi @MitchG_HunterCollege just a quick check whether your issue has been solved? If yes, please mark it as solved by pressing the :ballot_box_with_check: button under my reply :pray:

7

Hi Marcel. Problem solved (for now!). Thanks, Mitch

1 Like