# Rate-based model of coupled neuron populations to 
# perform discrimination using integral feedback control.
# Discrete attractor version, so r2 is determined by 
# the sum of 20 units, r2a to r2t.
par w1_2=0.065 w2_1=-0.7 tau=0.005 di=0.001 i2l=0
par tstim=1 t1=4.0 t2=10.0 tauapp=0.005 dtapp=0.001
par rsp=20 tauNMDA=0.03
par f1=20 f2=24 

init r1=0.0 r2=0.0 
init r2a=0.0 r2b=0.0 r2c=0.0 r2d=0.0 r2e=0.0 
init r2f=0.0 r2g=0.0 r2h=0.0 r2i=0.0 r2j=0.0 
init r2k=0.0 r2l=0.0 r2m=0.0 r2n=0.0 r2o=0.0 
init r2p=0.0 r2q=0.0 r2r=0.0 r2s=0.0 r2t=0.0 
# init sapp1=0.0 sapp2=0.0 
init s1=0.0 s2=0.0

sapp1=f1/((1+exp(-2*(t-t1)/dtapp))*(1+exp(-2*(t1+tstim-t)/dtapp)))
sapp2=f2/((1+exp(-2*(t-t2)/dtapp))*(1+exp(-2*(t2+tstim-t)/dtapp)))

ds1/dt=(sapp1-s1)/tauapp
ds2/dt=(sapp2-s2)/tauapp

i2=w1_2*r1*heav(r1)+r2a+r2b+r2c+r2d+r2e+r2f+r2g+r2h+r2i+r2j+r2k+r2l+r2m+r2n+r2o+r2p+r2q+r2r+r2s+r2t-i2l
dr2/dt=(-r2+i2*heav(i2))/tauNMDA
i1=rsp+w2_1*r2+s1+s2
dr1/dt=(-r1+i1)/tau
dr2a/dt=(-r2a+2*heav(r2-1))/tau
dr2b/dt=(-r2b+3*heav(r2-4))/tau
dr2c/dt=(-r2c+3*heav(r2-7))/tau
dr2d/dt=(-r2d+3*heav(r2-10))/tau
dr2e/dt=(-r2e+3*heav(r2-13))/tau
dr2f/dt=(-r2f+3*heav(r2-16))/tau
dr2g/dt=(-r2g+3*heav(r2-19))/tau
dr2h/dt=(-r2h+3*heav(r2-22))/tau
dr2i/dt=(-r2i+3*heav(r2-25))/tau
dr2j/dt=(-r2j+3*heav(r2-28))/tau
dr2k/dt=(-r2k+3*heav(r2-31))/tau
dr2l/dt=(-r2l+3*heav(r2-34))/tau
dr2m/dt=(-r2m+3*heav(r2-37))/tau
dr2n/dt=(-r2n+3*heav(r2-40))/tau
dr2o/dt=(-r2o+3*heav(r2-43))/tau
dr2p/dt=(-r2p+3*heav(r2-46))/tau
dr2q/dt=(-r2q+3*heav(r2-49))/tau
dr2r/dt=(-r2r+3*heav(r2-52))/tau
dr2s/dt=(-r2s+3*heav(r2-55))/tau
dr2t/dt=(-r2t+3*heav(r2-58))/tau

# dsapp1/dt=(s1-sapp1)/tauapp
# dsapp2/dt=(s2-sapp2)/tauapp
# ds1/dt=0.0
# ds2/dt=0.0
aux i2out=i2
aux di2=-r2+i2
# global 1 {t-t1} {s1=f1}
# global 2 {t-t1-tstim} {s1=0.0}
# global 3 {t-t2} {s2=f2}
# global 4 {t-t2-tstim} {s2=0}

@total=12,bound=100,dt=.001,dtmin=1e-5,dtmax=10,atoler=1e-4
@toler=1e-5,xhi=12,yhi=50,ylo=0 njmp=50,bell=0
@bell=off,nout=50                                                             
done