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)