Modeling synaptic depression(STD) in cortical network model

Description of problem
want to model synaptic depression in biophysical detailed synapse so to observe switching between two excitatory sub-populations of neurons. But switiching happens just once and it release probability Prel saturates. also confused where the Prel should be multiplied in case of NMDA synapses. it is brunel wang 2001 paper with added synaptic depression mechanism within selective groups .

Minimal code to reproduce problem

What you have aready tried

start_scope()

defaultclock.dt = 0.05*ms

N= 1000

Ne = 0.8*N

Ni =0.2*N

# parameters for both cells

VL=-70*mV

VThr= -50*mV

Vreset=-55*mV

# parameters for pyramidal cells

Cm_Pc =0.5*nF

gm_Pc= 25*nS

taurp_Pc =2*ms# refractory period

gAMPA_ext_Pc=2.08*nS

gAMPA_rec_Pc=0.104*(800.0/Ne)*nS

gNMDA_rec_Pc=0.328*(800.0/Ne)*nS

gGABA_Pc=1.25*(200.0/Ni)*nS

# parameters for the interneurons

Cm_In =0.2*nF

gm_In= 20*nS

taurp_In =1*ms

gAMPA_ext_In=1.62*nS

gAMPA_rec_In=0.081*(800.0/Ne)*nS

gNMDA_rec_In=0.258*(800.0/Ne)*nS

gGABA_In=0.973*(200.0/Ni)*nS

#parameters for the synapses of a fully connected network

Ce =Ne

Ci =Ni

Ve=0*mV #excitatory

Vi=-70*mV # inhibitory

Mg =1

tau_ampa=2*ms

tau_nmda_decay=100*ms

tau_gaba=10*ms

alpha=0.5/ms

tau_nmda_rise=2*ms

# writing equations for Pyramidal cell

eqs_Pc ='''

label : integer (constant)

dv/dt=(-gm_Pc*(v-VL)-Isyn)/Cm_Pc : volt (unless refractory)

Isyn= IAmpa_ext + IAmpa_rec + INmda_rec+ IGaba: ampere

IAmpa_ext=gAMPA_ext_Pc*(v-Ve)*s_ampa_ext: ampere

ds_ampa_ext/dt= -s_ampa_ext/tau_ampa :1

IAmpa_rec=gAMPA_rec_Pc*(v-Ve)*s_ampa_rec: ampere

ds_ampa_rec/dt= -s_ampa_rec/tau_ampa :1

INmda_rec=(gNMDA_rec_Pc*(v-Ve))/(1+Mg*exp(-0.062*v/mV)/3.57)*s_nmdaT: ampere

s_nmdaT:1

IGaba=gGABA_Pc*(v-Vi)*s_gaba: ampere

ds_gaba/dt= -s_gaba/tau_gaba :1

'''

# writing equation for Interneuron

eqs_In ='''

label : integer (constant)

dv/dt=(-gm_In*(v-VL)-Isyn)/Cm_In : volt (unless refractory)

Isyn= IAmpa_ext + IAmpa_rec + INmda_rec + IGaba: ampere

IAmpa_ext=gAMPA_ext_In* (v-Ve)*s_ampa_ext: ampere

ds_ampa_ext/dt= -s_ampa_ext/tau_ampa :1

IAmpa_rec=gAMPA_rec_In* (v-Ve)*s_ampa_rec: ampere

ds_ampa_rec/dt= -s_ampa_rec/tau_ampa :1

INmda_rec=(gNMDA_rec_In*(v-Ve))/(1+Mg*exp(-0.062*v/mV)/3.57)*s_nmdaT:ampere

s_nmdaT:1

IGaba=gGABA_In*(v-Vi)*s_gaba: ampere

ds_gaba/dt= -s_gaba/tau_gaba :1

'''
G_Pc = NeuronGroup(Ne,model=eqs_Pc,threshold='v>-50*mV',reset='v=Vreset', refractory =2*ms ,method='rk2')

G_Pc.v=VL

G_In = NeuronGroup(Ni,model=eqs_In,threshold='v>-50*mV',reset='v=Vreset', refractory =1*ms ,method='rk2')

G_In.v=VL
#create subgroups for excitatory and inhibitory neurons

#2 subpopulations

f= 0.08

number =int(N*f)# no of neurons in each subpopulation

popS1 = G_Pc[:number]

popS2 = G_Pc[number:2*number]

pop_nS=G_Pc[2*number:]

popS1.label = 0

popS2.label = 1

pop_nS.label = 2

G_In.label=3

w_plus = 2.1

w_minus = (0.8 - f*w_plus) / (0.8 - f)
# terms for depression

fD = 0.988# Depression factor

P0 = 1.0 # Baseline release probability

tauP = 1000*ms # Recovery time constant

#write synaptic equations

eqs_syn='''

s_nmdaT_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)

dPrel/dt = (P0 - Prel) / tauP *int(label_pre ==label_post)*int(label_pre < 2)+int(label_pre == label_post)*int(label_pre==2)*0*Hz+int(label_pre != label_post)*0*Hz:1(clock-driven)

w:1

'''

eqs_pre_glutamate = '''

s_ampa_rec += w*Prel*int(label_pre == label_post)*int(label_pre < 2) + int(label_pre == label_post)*int(label_pre==2)*w+ int(label_pre != label_post)*w

x += 1*Prel*int(label_pre == label_post)*int(label_pre < 2) + int(label_pre == label_post)*int(label_pre==2)*1+int(label_pre != label_post)*1

Prel=clip(Prel* fD, 0, 1)*int(label_pre == label_post)*int(label_pre < 2) + int(label_pre == label_post)*int(label_pre==2)*1+int(label_pre != label_post)*1

'''

eqs_pre_gaba = '''
s_gaba += 1.02
'''
#excitatory connections

# excitattory to inhibitory

# E to I

S3 = Synapses(G_Pc, G_In, model=eqs_syn, on_pre=eqs_pre_glutamate,delay=0.5*ms, method='rk2')

S3.connect()

S3.w[:] = 1

#excitatory to excitatory

# E to E

S4 = Synapses(G_Pc,G_Pc, model=eqs_syn, on_pre=eqs_pre_glutamate,delay=0.5*ms,method='rk2')

S4.connect('i != j')

S4.w[:] = 1

S4.w["label_pre == label_post and label_pre < 2"] = 1*w_plus

S4.w["label_pre != label_post and label_post < 2"] = 1*w_minus

S4.Prel = 1
Cext =800

#bacground spontaneous activity

rate =3.0*Hz

rate_input1 = 3.1*Hz

noise=TimedArray(np.r_[np.ones(240)], dt=25*ms)

noise_nS=PoissonInput(pop_nS, 's_ampa_ext', Cext, rate,'noise(t)')

noise_In=PoissonInput(G_In, 's_ampa_ext', Cext, rate,'noise(t)')

stimuli_1 = TimedArray(np.r_[np.ones(20),np.zeros(220)], dt=25*ms)

input_S1 = PoissonInput(popS1, 's_ampa_ext', Cext, rate, 'stimuli_1(t)')

input_S2 = PoissonInput(popS2, 's_ampa_ext', Cext, rate, 'stimuli_1(t)')

stimuli_0 = TimedArray(np.r_[np.zeros(20),np.ones(220)], dt=25*ms)

input_1 = PoissonInput(popS1, 's_ampa_ext', Cext, rate_input1, 'stimuli_0(t)')

input_2 = PoissonInput(popS2, 's_ampa_ext', Cext, rate_input1, 'stimuli_0(t)')

rate_S1=PopulationRateMonitor(popS1)

rate_S2=PopulationRateMonitor(popS2)
run_time=6*second

run(run_time, report = "text")

plt.plot(rate_S1.t / ms, rate_S1.smooth_rate(width=25*ms) / Hz, label='pop1',color="red",linewidth=2)

plt.plot(rate_S2.t / ms, rate_S2.smooth_rate(width=25*ms) / Hz, label='pop2',color="blue",linewidth=2)

show()

the other way i had tried modelling depression

eqs_syn='''

s_nmdaT_post = w*s_nmda*Prel*int(label_pre == label_post)*int(label_pre < 2) + int(label_pre == label_post)*int(label_pre==2)*w*s_nmda +int(label_pre != label_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)

dPrel/dt = (P0 - Prel) / tauP *int(label_pre ==label_post)*int(label_pre < 2)+int(label_pre == label_post)*int(label_pre==2)*0*Hz+int(label_pre != label_post)*0*Hz:1(clock-driven)

w:1

'''

eqs_pre_glutamate = '''

w=w*Prel*int(label_pre == label_post)*int(label_pre < 2) + int(label_pre == label_post)*int(label_pre==2)*w+ int(label_pre != label_post)*w

s_ampa_rec += w

x += 1

Prel=clip(Prel* fD, 0, 1)*int(label_pre == label_post)*int(label_pre < 2) + int(label_pre == label_post)*int(label_pre==2)*1+int(label_pre != label_post)*1

'''

eqs_pre_gaba = '''

s_gaba += 1.02

'''

Expected output (if relevant)

Actual output (if relevant)

Full traceback of error (if relevant)