# Brian1 to brian2 conversion; issue could be with update\_nmda

**URL:** <https://brian.discourse.group/t/brian1-to-brian2-conversion-issue-could-be-with-update-nmda/981>\
**Category:** Support\
**Created:** [28 April 2023 12:03 UTC](https://brian.discourse.group/t/brian1-to-brian2-conversion-issue-could-be-with-update-nmda/981 "2023-04-28T12:03:21Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![selvamaran](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/selvamaran/32/599_2.png) [@selvamaran](https://brian.discourse.group/u/selvamaran)\
**Post date:** [28 April 2023 12:03 UTC](https://brian.discourse.group/t/brian1-to-brian2-conversion-issue-could-be-with-update-nmda/981/1 "2023-04-28T12:03:21Z")

</div>

# 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

```python
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

```python
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

---

<div class="post-metadata">

**Author:** ![genevievefahey](https://avatars.discourse-cdn.com/v4/letter/g/f14d63/32.png) [@genevievefahey](https://brian.discourse.group/u/genevievefahey)\
**Post date:** [3 May 2023 00:31 UTC](https://brian.discourse.group/t/brian1-to-brian2-conversion-issue-could-be-with-update-nmda/981/2 "2023-05-03T00:31:02Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [3 May 2023 09:57 UTC](https://brian.discourse.group/t/brian1-to-brian2-conversion-issue-could-be-with-update-nmda/981/3 "2023-05-03T09:57:41Z")

</div>

> [@genevievefahey](#):
>
> I copied your code as it is but I noticed that there are a lot of syntax issues and missing multiplication signs.

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 😊 ), please post Python code like this:

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

````

i.e. enclosed in triple backticks 🙏.

---

<div class="post-metadata">

**Author:** ![selvamaran](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/selvamaran/32/599_2.png) [@selvamaran](https://brian.discourse.group/u/selvamaran)\
**Post date:** [8 May 2023 04:22 UTC](https://brian.discourse.group/t/brian1-to-brian2-conversion-issue-could-be-with-update-nmda/981/4 "2023-05-08T04:22:05Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![genevievefahey](https://avatars.discourse-cdn.com/v4/letter/g/f14d63/32.png) [@genevievefahey](https://brian.discourse.group/u/genevievefahey)\
**Post date:** [10 May 2023 02:49 UTC](https://brian.discourse.group/t/brian1-to-brian2-conversion-issue-could-be-with-update-nmda/981/5 "2023-05-10T02:49:44Z")

</div>

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:

```auto
 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)

---

<div class="post-metadata">

**Author:** ![selvamaran](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/selvamaran/32/599_2.png) [@selvamaran](https://brian.discourse.group/u/selvamaran)\
**Post date:** [11 May 2023 13:04 UTC](https://brian.discourse.group/t/brian1-to-brian2-conversion-issue-could-be-with-update-nmda/981/6 "2023-05-11T13:04:02Z")

</div>

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:

```auto
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](https://brian2.readthedocs.io/en/stable/examples/frompapers.Brunel_Wang_2001.html?highlight=nmda#sample-specific-persistent-activity)). The model is similar and the model has been coded without update\_nmda for NMDA synapse.