Brian1 to brian2 conversion; issue could be with update_nmda

Support
1

Description of problem

I am trying to convert the code from brian1 to brian2, related to the paper
“A computational model of Major Depression: The role of glutamate dysfunction on cingulo-frontal network dynamics” Ramirez-Mahaluf J.P., Roxin A., Mayberg H.S. and Compte A. Cerebral Cortex, 2015

Minimal code to reproduce problem

NE=800                    
NI=200                     

#Biophysical parameters
tauav = 2*ms            
tauad = 2*ms            # tau AMPA decay on dlPFC
taun = 100*ms          # tau NMDA decay
taux = 2*ms              # tau NMDA rise
taug = 10*ms            # tau GABA decay
Vt  =-50*mvolt           # spike threshold
Vr  =-55*mvolt           # reset value
Elv  =-70*mvolt         
El  =-70*mvolt            
refE= 2*ms # refractory periods piramidal cell
refI= 1*ms                 # refractory period inhibitory cell
cmE= 500*pF             #capacitance piramidal cel   
cmI= 200*pF              #capacitance interneuron  
tauE =20*ms              #tau piramidal cel
tauI =10*ms               #tau interneuron 
alpha =0.5*kHz
S=1	
N=800/NE                  # Factor for rescaling the weights according to the number of neurons 

#Connection parameters
wgEEN = 0.001761*(1/S)*N     #weight  excitatory to excitatory through NMDA
wgEEA = 0.0009454*(1/S)*N    #weight  excitatory to excitatory through AMPA
wgEIN = 0.0012*(1/S)*N         #weight  excitatory to inhibitory through NMDA
wgEIA = 0.0004*(1/S)*N         #weight  excitatory to inhibitory through AMPA
wgIEG = 0.005*(1/S)*N          #weight inhibitory to excitatory through GABA
wgIIG = 0.004865*(1/S)*N     #weight inhibitory to  inhibitory through GABA
wgEIA1 = 0.0004*(1/S)*N      #weight vACC excitatory to dlPFC inhibitory through NMDA
wgEIA2 = 0.0004*(1/S)*N      #weight  dlPFC excitatory to vACC excitatory through NMDA

#equations excitatory cell vACC 
eqsE1 = '''
dV/dt = (-gea*V-gen*V/(1+exp(-V/Ven)/3.57)-gi*(V+70*mV)-(V-Elv))/(tauE)  : volt
dgea/dt = -gea/(tauav)           : 1
dgi/dt = -gi/(taug)              : 1
dspre/dt = -spre/(taun)+alpha*xpre*(1-spre) : 1
dxpre/dt = -xpre/(taux)                     : 1
gen: 1

'''
#equations inhibitory cell vACC 
eqsI1 = '''
dV/dt = (-gea*V-gen*V/(1+exp(-V/Ven)/3.57)-gi*(V+70*mV)-(V-El))/(tauI)    : volt
dgea/dt = -gea/(tauav)           : 1
dgi/dt = -gi/(taug)             : 1
dspre/dt = -spre/(taun)+alpha*xpre*(1-spre) : 1
dxpre/dt = -xpre/(taux)                     : 1
gen: 1


'''
#equations excitatory cell dlPFC 
eqsE2 = '''
dV/dt = (-gea*V-gen*V/(1+exp(-V/Ven)/3.57)-gi*(V+70*mV)-(V-El))/(tauE)   : volt
dgea/dt = -gea/(tauad)           : 1
dgi/dt = -gi/(taug)             : 1
dspre/dt = -spre/(taun)+alpha*xpre*(1-spre) : 1
dxpre/dt = -xpre/(taux)                     : 1
gen: 1


'''
#equations inhibitory cell dlPFC 
eqsI2 = '''
dV/dt = (-gea*V-gen*V/(1+exp(-V/Ven)/3.57)-gi*(V+70*mV)-(V-El))/(tauI)    : volt
dgea/dt = -gea/(tauad)           : 1
dgi/dt = -gi/(taug)             : 1
dspre/dt = -spre/(taun)+alpha*xpre*(1-spre) : 1
dxpre/dt = -xpre/(taux)                     : 1
gen: 1

'''

#Populations of neurons: 
Pev = NeuronGroup(NE, model= eqsE1, threshold=Vt, reset= Vr, refractory=refE)      #vACC excitatory neurons 
Piv = NeuronGroup(NI, model= eqsI1, threshold=Vt, reset= Vr, refractory=refI)          #vACC inhibitory neurons

Ped = NeuronGroup(NE, model= eqsE2, threshold=Vt, reset= Vr, refractory=refE)       #dlPFC excitatory neurons
Pid = NeuronGroup(NI, model= eqsI2, threshold=Vt, reset= Vr, refractory=refI)          #dlPFC inhibitory neurons

#Connection NMDA:
selfnmda_v = IdentityConnection(Pev, Pev, 'xpre', weight=1.0) #NMDA connections, excitatory to excitatory neurons in vACC
selfnmda_d = IdentityConnection(Ped, Ped, 'xpre', weight=1.0) #NMDA connections, excitatory to excitatory neurons in dlPF

#Connections AMPA and GABA:
Ceeav = Connection(Pev, Pev, 'gea', structure='dense') #AMPA connections, excitatory to excitatory neurons in vACC 
Ceiav = Connection(Pev, Piv, 'gea', structure='dense') #AMPA connections, excitatory to inhibitory neurons in vACC
Ciev = Connection(Piv, Pev, 'gi', structure='dense') # GABA connections, inhibitory to excitatory neurons in vACC
Ciiv = Connection(Piv, Piv, 'gi', structure='dense') # GABA connections, excitatory to excitatory neurons in vACC

Ceead = Connection(Ped, Ped, 'gea', structure='dense')#AMPA connections, excitatory to excitatory neurons in dlPFC 
Ceiad = Connection(Ped, Pid, 'gea', structure='dense') #AMPA connections, excitatory to inhibitory neurons in dlPFC
Cied = Connection(Pid, Ped, 'gi', structure='dense')# GABA connections, inhibitory to excitatory neurons in dlPFC
Ciid = Connection(Pid, Pid, 'gi', structure='dense')# GABA connections, excitatory to excitatory neurons in dlPFC

Ceiav1 = Connection(Pev, Pid, 'gea' )#AMPA connections, excitatory neurons in vACC target inhibitory neurons in dlPFC
Ceiad1 = Connection(Ped, Piv, 'gea' )#AMPA connections excitatory neurons in dlPFC target inhibitory neurons in vACC

Ceeav.connect_random(Pev, Pev, S, weight=wgEEA)  #AMPA connections, excitatory to excitatory neurons in vACC  
Ceiav.connect_random(Pev, Piv, S, weight=wgEIA)  #AMPA connections, excitatory to inhibitory neurons in vACC
Ciev.connect_random(Piv, Pev, S, weight=wgIEG)  # GABA connections, inhibitory to excitatory neurons in vACC
Ciiv.connect_random(Piv, Piv, S, weight=wgIEG)  # GABA connections, excitatory to excitatory neurons in vACC

Ceead.connect_random(Ped, Ped, S, weight=wgEEA) #AMPA connections, excitatory to excitatory neurons in dlPFC 
Ceiad.connect_random(Ped, Pid, S, weight=wgEIA) #AMPA connections, excitatory to inhibitory neurons in dlPFC
Cied.connect_random(Pid, Ped, S, weight=wgIEG) # GABA connections, inhibitory to excitatory neurons in dlPFC
Ciid.connect_random(Pid, Pid, S,weight=wgIIG) # GABA connections, excitatory to excitatory neurons in dlPFC

Ceiav1.connect_random(Pev, Pid, S, weight=wgEIA1) #AMPA connections, excitatory neurons in vACC target inhibitory neurons in dlPFC
Ceiad1.connect_random(Ped, Piv, S, weight=wgEIA2) #AMPA connections excitatory neurons in dlPFC target inhibitory neurons in vACC


#NMDA synapses
E_nmda_v = asarray(Pev.spre)
E_nmda_d = asarray(Ped.spre)
E_gen_v = asarray(Pev.gen)
E_gen_d = asarray(Ped.gen)
I_gen_v = asarray(Piv.gen)
I_gen_d = asarray(Pid.gen)

#Calculate NMDA contributions
@network_operation(when='start')
def update_nmda():
 E_gen_v[:] = wgEEN/wgEEA * numpy.dot(E_nmda_v,Ceeav.W)
 I_gen_v[:] = wgEIN/wgEIA * numpy.dot(E_nmda_v,Ceiav.W)
 E_gen_d[:] = wgEEN/wgEEA * numpy.dot(E_nmda_d,Ceead.W)
 I_gen_d[:] = wgEIN/wgEIA * numpy.dot(E_nmda_d,Ceiad.W)


#External noise:
extinput1E=PoissonGroup(NE,rates=1800*Hz)
extinput1I=PoissonGroup(NI,rates=1800*Hz)

input1_coE=IdentityConnection(extinput1E,Pev,'gea',weight=0.082708)
input1_coI=IdentityConnection(extinput1I,Piv,'gea',weight=0.081)

extinput2E=PoissonGroup(NE,rates=1800*Hz)
extinput2I=PoissonGroup(NI,rates=1800*Hz)

input2_coE=IdentityConnection(extinput2E,Ped,'gea',weight=0.082708)
input2_coI=IdentityConnection(extinput2I,Pid,'gea',weight=0.081)

#run
run(100*second)

What you have aready tried

Pev = NeuronGroup(NE, model= eqsE1, threshold=Vt, reset= Vr, refractory=refE)      #vACC excitatory neurons 
Piv = NeuronGroup(NI, model= eqsI1, threshold=Vt, reset= Vr, refractory=refI)          #vACC inhibitory neurons

Ped = NeuronGroup(NE, model= eqsE2, threshold=Vt, reset= Vr, refractory=refE)       #dlPFC excitatory neurons
Pid = NeuronGroup(NI, model= eqsI2, threshold=Vt, reset= Vr, refractory=refI)          #dlPFC inhibitory neurons


selfnmda_v = Synapses(Pev, Pev,   on_pre='xpre += 1.0') 
selfnmda_v.connect(j='i')
selfnmda_d = Synapses(Ped, Ped,  on_pre='xpre += 1.0')
selfnmda_d.connect(j='i')


Ceeav = Synapses(Pev, Pev, on_pre='gea += wgEEA') #AMPA connections, excitatory to excitatory neurons in vACC 
Ceiav = Synapses(Pev, Piv,  on_pre='gea += wgEIA') #AMPA connections, excitatory to inhibitory neurons in vACC
Ciev = Synapses(Piv, Pev, on_pre='gi += wgIEG') # GABA connections, inhibitory to excitatory neurons in vACC
Ciiv = Synapses(Piv, Piv,  on_pre='gi += wgIIG') # GABA connections, excitatory to excitatory neurons in vACC

Ceead = Synapses(Ped, Ped, on_pre='gea += wgEEA')#AMPA connections, excitatory to excitatory neurons in dlPFC 
Ceiad = Synapses(Ped, Pid, on_pre='gea += wgEIA') #AMPA connections, excitatory to inhibitory neurons in dlPFC
Cied = Synapses(Pid, Ped,   on_pre='gi += wgIEG')# GABA connections, inhibitory to excitatory neurons in dlPFC
Ciid = Synapses(Pid, Pid,   on_pre='gi += wgIIG')# GABA connections, excitatory to excitatory neurons in dlPFC

Ceiav1 = Synapses(Pev, Pid,  on_pre='gea += wgEIA1')#AMPA connections, excitatory neurons in vACC target inhibitory neurons in dlPFC
Ceiad1 = Synapses(Ped, Piv,  on_pre='gea += wgEIA2' )#AMPA connections excitatory neurons in dlPFC target inhibitory neurons in vACC

Ceeav.connect(p=S)
Ceiav.connect(p=S)
Ciev.connect(p=S)
Ciiv.connect(p=S)

Ceead.connect(p=S)
Ceiad.connect(p=S)
Cied.connect(p=S)
Ciid.connect(p=S)

Ceiav1.connect(p=S)
Ceiad1.connect(p=S)


#NMDA synapses

E_nmda_v = asarray(Pev.spre)
E_nmda_d = asarray(Ped.spre)

E_gen_v = asarray(Pev.gen)
E_gen_d = asarray(Ped.gen)
#I_gen_v = asarray(Piv.gen)
#I_gen_d = asarray(Pid.gen)
I_gen_v = asarray(Pev.gen)
I_gen_d = asarray(Ped.gen)

#Calculate NMDA contributions

@network_operation(when='start')
def update_nmda():
    E_gen_v[:] = wgEEN/wgEEA * numpy.dot(E_nmda_v,wgEEA)
    I_gen_v[:] = wgEIN/wgEIA * numpy.dot(E_nmda_v,wgEIA)
    E_gen_d[:] = wgEEN/wgEEA * numpy.dot(E_nmda_d,wgEEA)
    I_gen_d[:] = wgEIN/wgEIA * numpy.dot(E_nmda_d,wgEIA)

#External noise:
extinput1E=PoissonGroup(NE,rates=1800*Hz)
extinput1I=PoissonGroup(NI,rates=1800*Hz)

input1_coE=Synapses(extinput1E,Pev,on_pre='gea += 0.082708')
input1_coE.connect(j='i')
input1_coI=Synapses(extinput1I,Piv,on_pre='gea += 0.081')
input1_coI.connect(j='i')

extinput2E=PoissonGroup(NE,rates=1800*Hz)
extinput2I=PoissonGroup(NI,rates=1800*Hz)

input2_coE=Synapses(extinput2E,Ped,on_pre='gea += 0.082708')
input2_coE.connect(j='i')
input2_coI=Synapses(extinput2I,Pid,on_pre='gea += 0.081')
input2_coI.connect(j='i')

Expected output (if relevant)

The model is running. But it gives the similar kind of results when the NMDA synaptic strength is set to zero which shouldn’t happen. I want to be sure my update_nmda change is right becuase i think that might be the issue.
Thanks,
Selva

2

I copied your code as it is but I noticed that there are a lot of syntax issues and missing multiplication signs. Do you have a copy of the python code that you have run to get an output?

3

This is because the code was interpreted as markdown code, where multiplication signs need are used for italics and bold (written as *italics* and **bold**). I edited the original post so that it shows the Python code.
@selvamaran (and everyone else of course :blush: ), please post Python code like this:

```
# Here comes the code
print("I am Python code")
```

i.e. enclosed in triple backticks :pray:.

4

Thank you Marcel (@mstimberg ). Sure, in future i will enclose the code in triple brackets. Hai Genevieve (@genevievefahey ), let me know whether you were able to use the code now.
I did spend sometime exploring the issue and it seems as I mentioned the problem seems to be with update_nmda. When I track the gen state variable its negligible.
I looked at other codes that use update_nmda and seems each one use different style of coding and so i find it bit hard to understand.

I also wrote a simple single attractor network by simplifying the original code and was able to reproduce attractor state but even in it update_syntax syntax is not right. Basically not sure how to modify the code to sum all the NMDA conductance in update_NMDA in Brian 2.

5

Hi, I’ve had a look at your code and I am able to get it running. I wondered if it may have been to do with the IdentityConnection conversion or the update_nmda array dimensions but now I am unsure if this is the case and it is difficult for me to know what to look for to determine this.

Do you have a short simulation’s set of results I could compare against? What value you would expect for these .gen variables to obtain? If I can run() for only a couple of seconds that would be even better.

I can see that the original code has tasks written under the part that you provided. I am unsure if you are supposed to see non-negligible results from the update_nmda function until after the PoissonGroup rates have been adjusted.

As an example, if I run(1*second) the first few lines of Pev.gen variables look like this:

 3.54648573e-08 3.54648573e-08 3.54648573e-08 3.54648573e-08
 3.54648573e-08 3.54648573e-08 3.54648573e-08 3.41180360e-04
 2.55225101e-06 3.54648573e-08 7.50304089e-08 3.54648573e-08
 3.54648573e-08 3.54648573e-08 3.54648573e-08 3.54648573e-08

They are small, but one is only 3.4e-4, so that might not be negligible to your model. I’m not sure.

(setting np.random.seed(0) will fix the RNG so it produces the same results each time it runs, otherwise they may be slightly different)

6

Thanks @genevievefahey.
The issue is i am not sure whether i am writing the proper syntax for update_nmda (For example lines (sE=sum(E_nmda_v), E_gen_v[:] = wgEEN/wgEEA * sE) - not sure is this right). I have stripped the code to just single excitatory-inhibitory network (Even separatly reduced further with single excitatory network alone). I am able to reproduce the attractor state, NMDA being activated etc, but issue is that i need to scale down the impact of NMDA current.

Below is the code for excitatory-inhibitory network:

from brian2 import *
from numpy.fft import rfft,irfft
from scipy.io import savemat
import numpy
import matplotlib.pyplot as plt

start_scope()
#Network parameters

NE=800                    # Excitatory neurons, for fast simulations in the article we use 80
NI=200                     # Inhibitory neurons, for fast simulations in the article we use 20
# NE=80                    # Excitatory neurons, for fast simulations in the article we use 80
# NI=20                     # Inhibitory neurons, for fast simulations in the article we use 20


#Biophysical parameters
tauav = 2*ms            # tau AMPA decay on vACC, this parameter was used to simulate MDD. Mild MDD (2.5%) = 2.05; Moderate MDD (5%) = 2.1; Severe MDD (7.5%) = 2.15
taun = 100*ms          # tau NMDA decay
taux = 2*ms              # tau NMDA rise
taug = 10*ms            # tau GABA decay
Vt  ='V > -50*mV'           # spike threshold
Vr  ='V = -55*mV'           # reset value
Elv  =-70*mvolt         # resting potential ventral ACC, this parameter was used to simulate SSRI treatment.
El  =-70*mvolt           # resting potential dlPFC
Ven = 16.129*mV      
refE= 2*ms       # refractory periods piramidal cell
refI= 1*ms                 # refractory period inhibitory cell
cmE= 500*pF             #capacitance piramidal cel   
cmI= 200*pF              #capacitance interneuron  
tauE =20*ms              #tau piramidal cel
tauI =10*ms               #tau interneuron 
alpha =0.5*kHz



S=1	#Connectivity sparsensess; S=1, all-to-all connectivity was used in the article; use S<1 for sparse random connectivity
N=800/NE                  # Factor for rescaling the weights according to the number of neurons 
#Connection parameters
wgEEN = 0.001761*(1/S)*N     #weight  excitatory to excitatory through NMDA
wgEEA = 0.0009454*(1/S)*N    #weight  excitatory to excitatory through AMPA
wgEIN = 0.0012*(1/S)*N         #weight  excitatory to inhibitory through NMDA
wgEIA = 0.0004*(1/S)*N        #weight  excitatory to inhibitory through AMPA
wgIEG = 0.005*(1/S)*N          #weight inhibitory to excitatory through GABA
wgIIG = 0.004865*(1/S)*N     #weight inhibitory to  inhibitory through GABA



#equations excitatory cell vACC 
eqsE1 = '''
dV/dt = (-gea*V-gen*V/(1+exp(-V/Ven)/3.57)-gi*(V+70*mV)-(V-Elv))/(tauE)  : volt (unless refractory)
dgea/dt = -gea/(tauav)           : 1
dgi/dt = -gi/(taug)              : 1
dspre/dt = -spre/(taun)+alpha*xpre*(1-spre) : 1
dxpre/dt = -xpre/(taux)                     : 1
gen: 1

'''
#equations inhibitory cell vACC 
eqsI1 = '''
dV/dt = (-gea*V-gen*V/(1+exp(-V/Ven)/3.57)-gi*(V+70*mV)-(V-El))/(tauI)    : volt (unless refractory)
dgea/dt = -gea/(tauav)           : 1
dgi/dt = -gi/(taug)             : 1
dspre/dt = -spre/(taun)+alpha*xpre*(1-spre) : 1
dxpre/dt = -xpre/(taux)                     : 1
gen: 1


'''



Pev = NeuronGroup(NE, model= eqsE1, threshold=Vt, reset= Vr, refractory=refE)      #vACC excitatory neurons 
Pev.V = -55*mV
Piv = NeuronGroup(NI, model= eqsI1, threshold=Vt, reset= Vr, refractory=refI)          #vACC inhibitory neurons
Piv.V = -55*mV

selfnmda_v = Synapses(Pev, Pev,   on_pre='xpre += 1.0') 
selfnmda_v.connect(j='i')


Ceeav = Synapses(Pev, Pev, on_pre='gea += wgEEA') #AMPA connections, excitatory to excitatory neurons in vACC 
Ceiav = Synapses(Pev, Piv,  on_pre='gea += wgEIA') #AMPA connections, excitatory to inhibitory neurons in vACC
Ciev = Synapses(Piv, Pev, on_pre='gi += wgIEG') # GABA connections, inhibitory to excitatory neurons in vACC
Ciiv = Synapses(Piv, Piv,  on_pre='gi += wgIIG') # GABA connections, inhibitory to inhibitory neurons in vACC

Ceeav.connect(p=S)
Ceiav.connect(p=S)
Ciev.connect(p=S)
Ciiv.connect(p=S)



# NMDA synapses
E_nmda_v = asarray(Pev.spre)
E_gen_v = asarray(Pev.gen)
I_gen_v = asarray(Piv.gen)


# Calculate NMDA contributions

@network_operation(when='start')
def update_nmda():
    sE=sum(E_nmda_v)
    E_gen_v[:] = 0.0012 *wgEEN/wgEEA * sE
    I_gen_v[:] = 0.0005 *wgEIN/wgEIA * sE

    
#External noise:
extinput1E=PoissonGroup(NE,rates=1800*Hz)
extinput1I=PoissonGroup(NI,rates=1800*Hz)

input1_coE=Synapses(extinput1E,Pev,on_pre='gea += 0.082708')
input1_coE.connect(j='i')
input1_coI=Synapses(extinput1I,Piv,on_pre='gea += 0.081')
input1_coI.connect(j='i')



# #Sadness task,
exttaskinput1_on=1000*ms
exttaskinput1_off=1250*ms
exttaskinput1E=PoissonGroup(NE,rates= '(t>exttaskinput1_on)*(t<exttaskinput1_off)*200*Hz')

taskinput1_coE=Synapses(exttaskinput1E,Pev,  on_pre='gea += 0.0955')
taskinput1_coE.connect(j='i')

#  cognitive signal to vACC
exttaskinput2_on=2750*ms
exttaskinput2_off=3000*ms
exttaskinput1I=PoissonGroup(NI,rates= '(t>exttaskinput2_on)*(t<exttaskinput2_off)*200*Hz')

taskinput1_coI=Synapses(exttaskinput1I, Piv,  on_pre='gea += 0.2')
taskinput1_coI.connect(j='i')

#Save files
Miv = SpikeMonitor(Piv)
Mev = SpikeMonitor(Pev)

Mv=PopulationRateMonitor(Pev)
Mivp=PopulationRateMonitor(Piv)

Mv_V = StateMonitor(Pev, 'V', record=True)
Miv_V = StateMonitor(Piv, 'V', record=True)

Mev_gen = StateMonitor(Pev, 'gen', record=True)
Mev_gea = StateMonitor(Pev, 'gea', record=True)
Mev_gi = StateMonitor(Pev, 'gi', record=True)
#run
run(3.5*second) 

subplot(2,1,1)
plot(Mev.t/ms, Mev.i, '.k')

subplot(2,1,2)
plot(Miv.t/ms, Miv.i, '.k')

subplot(2,1,1)
plt.plot(Mv.smooth_rate(window='flat', width=100*ms)[1000:40000])
subplot(2,1,2)
plt.plot(Mivp.smooth_rate(window='flat', width=100*ms)[1000:40000])

plot(Mev_gen.t/ms, Mev_gen.gen[0])
xlabel('Time (ms)')
ylabel('gen')
plot(Mev_gen.t/ms, Mev_gen.gen[0]*Mv_V.V[0]/(1+exp(-Mv_V.V[0]/Ven)/3.57))
xlabel('Time (ms)')
ylabel('I_nmda')

plt.figure(figsize = (32,8))
plot(Mev_gen.t/ms, (-Mev_gea.gea[0]*Mv_V.V[0]-Mev_gen.gen[0]*Mv_V.V[0]/(1+exp(-Mv_V.V[0]/Ven)/3.57)-Mev_gi.gi[0]*(Mv_V.V[0]+70*mV)-(Mv_V.V[0]-Elv))/(tauE))
plot(Mev_gen.t/ms, (-Mev_gea.gea[0]*Mv_V.V[0]-Mev_gi.gi[0]*(Mv_V.V[0]+70*mV)-(Mv_V.V[0]-Elv))/(tauE))
# plot( M.spre[10])
xlabel('Time (ms)')
ylabel('gen')

Iam also looking at the code related to Brunel_Wang_2001 (Example: Brunel_Wang_2001 — Brian 2 2.5.1 documentation). The model is similar and the model has been coded without update_nmda for NMDA synapse.