# RANDOM, SPARSELY CONNECTED NETWORK 4/10/02
# Produces "episodic" rhythm with somewhat irregular
# firing within each episode.
# Used to generate Figure 3 in J. Neuro. paper.
# There are 10 STN and GPe cells.
# Each STN cell receives input from two GPe cells chosen at random.
# Each GPe cell receives input from one STN cell.
# There is all to all coupling among the GPe cells.


########  STN Parameters  #######

p vl=-60,vna=55.,vk=-80.,thetam=30,sigmam=15
p gl=2.25,gna=37.5,gk=45,tn=1.,th=0.05
p gahp=9.,gca=.5,vca=140.,k1=15.,eps=5e-05
p kca=22.5,thetas=39.,sigmas=8.,xp=1.,i=0.
p thetah=-39,sigmah=3.1,thetan=-32.,sigman=-8.
p taun0=1,taun1=100.,thn=80.,sigmant=26.
p tauh0=1,tauh1=500,thh=57.,sigmaht=3.,phi=.75
p thetaa=-63.,sigmaa=-7.8,gt=.5,phir=.5
p thetar=-67,sigmar=2.,taur0=7.1,taur1=17.5,thr=-68,sigmart=2.2
##############################
#these parameters are not needed
p eps1=.01
p root1=-60,root2=-40,root3=-10,root4=-35,add=0.6
p scale=10000.,gd=1.0
###############################
p alpha=5,beta=1.,ab=-30.
p gGtoS=2.5,vGtoS=-100
p thetab=.25,sigmab=-.07


#######  GPe Parameters  #######

p gnag=120.,gkg=30.,gahpg=30.,gtg=.5,gcag=.15,glg=.1
p vnag=55.,vkg=-80.,vcag=120.,vlg=-55.
p thag=-57.,sigag=2.,thsg=-35.,sigsg=2.
p thrg=-70.,sigrg=-2.,taurg=30.
p thmg=-37.,sigmg=10.
p thng=-50.,signg=14.
p taun0g=.05,taun1g=.27,thngt=-40,sng=-12
p thhg=-58,sighg=-12
p tauh0g=.05,tauh1g=.27,thhgt=-40,shg=-12
p k1g=30.,kcag=20.,epsg=0.0001
p phig=1.,phing=.05,phihg=.05
p iapp=-1.2,gGtoG=0.0,vGtoG=-100.
p gStoG=0.016,vStoG=0,alphag=2,betag=.08,abg=-20


#######  STN Functions  #######

sinf(v)=1./(1.+exp(-(v+thetas)/sigmas))
minf(v)=1./(1.+exp(-(v+thetam)/sigmam))
hinf(v)=1./(1.+exp((v-thetah)/sigmah))
ninf(v)=1./(1.+exp((v-thetan)/sigman))
taun(v)=taun0+taun1/(1+exp((v+thn)/sigmant))
tauh(v)=tauh0+tauh1/(1+exp((v+thh)/sigmaht))
rinf(v)=1/(1+exp((v-thetar)/sigmar))
taur(v)=taur0+taur1/(1+exp((v+thr)/sigmart))
ainf(v)=1/(1+exp((v-thetaa)/sigmaa))
binf(r)=1/(1+exp((r-thetab)/sigmab))-1/(1+exp(-thetab/sigmab))


####### GPe Functions  #######

ainfg(v)=1/(1+exp(-(v-thag)/sigag))
sinfg(v)=1/(1+exp(-(v-thsg)/sigsg))
rinfg(v)=1/(1+exp(-(v-thrg)/sigrg))
minfg(v)=1./(1.+exp(-(v-thmg)/sigmg))
ninfg(v)=1./(1.+exp(-(v-thng)/signg))
taung(v)=taun0g+taun1g/(1+exp(-(v-thngt)/sng))
hinfg(v)=1./(1.+exp(-(v-thhg)/sighg))
tauhg(v)=tauh0g+tauh1g/(1+exp(-(v-thhgt)/shg))


####### STN Currents  #######

il(v)=gl*(v-vl)
ina(v,h)=gna*(minf(v))^3*h*(v-vna)
ik(v,n)=gk*n^4*(v-vk)
iahp(v,ca)=gahp*(v-vk)*ca/(ca+k1)
ica(v)=gca*((sinf(v))^2)*(v-vca)
it(v,r)=gt*(ainf(v)**3)*(binf(r)^2)*(v-vca)
#############################
#The choice of T-current is a bit unusual.
#Here we use binf(r) instead of just r.
#This gives a more acurrate rebound.
#############################
isyn1=gGtoS*(sg10+sg2)*(v1-vGtoS)
isyn2=gGtoS*(sg7+sg9)*(v2-vGtoS)
isyn3=gGtoS*(sg4+sg6)*(v3-vGtoS)
isyn4=gGtoS*(sg5+sg3)*(v4-vGtoS)
isyn5=gGtoS*(sg1+sg8)*(v5-vGtoS)
isyn6=gGtoS*(sg2+sg9)*(v6-vGtoS)
isyn7=gGtoS*(sg5+sg6)*(v7-vGtoS)
isyn8=gGtoS*(sg1+sg7)*(v8-vGtoS)
isyn9=gGtoS*(sg4+sg10)*(v9-vGtoS)
isyn10=gGtoS*(sg3+sg8)*(v10-vGtoS)


#######  GPe Currents  #######

itg(vg,rg)=gtg*(ainfg(vg)^3)*rg*(vg-vcag)
inag(vg,hg)=gnag*(minfg(vg)^3)*hg*(vg-vnag)
ikg(vg,ng)=gkg*(ng^4)*(vg-vkg)
iahpg(vg,cag)=gahpg*(vg-vkg)*cag/(cag+k1g)
icag(vg)=gcag*((sinfg(vg))^2)*(vg-vcag)
ilg(vg)=glg*(vg-vlg)

isyng1=gStoG*s1*(vg1-vStoG)
isyng2=gStoG*s2*(vg2-vStoG)
isyng3=gStoG*s3*(vg3-vStoG)
isyng4=gStoG*s4*(vg4-vStoG)
isyng5=gStoG*s5*(vg5-vStoG)
isyng6=gStoG*s6*(vg6-vStoG)
isyng7=gStoG*s7*(vg7-vStoG)
isyng8=gStoG*s8*(vg8-vStoG)
isyng9=gStoG*s9*(vg9-vStoG)
isyng10=gStoG*s10*(vg10-vStoG)
stot=sg1+sg2+sg3+sg4+sg5+sg6+sg7+sg8+sg9+sg10
isyngg(vg,sg)=gGtoG*stot*(vg-vGtoG)


#######  STN Equations  #######

v[1..10]'=-(il(v[j])+ina(v[j],h[j])+ik(v[j],n[j])+iahp(v[j],ca[j])+ica(v[j])+it(v[j],r[j]))-isyn[j]
h[1..10]'=phi*( hinf(v[j])-h[j] )/tauh(v[j])
n[1..10]'=phi*( ninf(v[j])-n[j] )/taun(v[j])
r[1..10]'=phir*(rinf(v[j])-r[j])/taur(v[j])
ca[1..10]'=phi*eps*(-ica(v[j])-it(v[j],r[j]) - kca*ca[j])
s[1..10]'=alpha*(1-s[j])*sinf(v[j]+ab)-beta*s[j]


#######  GPe Equations  #######

vg[1..10]'= -(itg(vg[j],rg[j])+inag(vg[j],hg[j])+ikg(vg[j],ng[j])+iahpg(vg[j],cag[j])+icag(vg[j])+ilg(vg[j])) +iapp -isyngg(vg[j],stot)-isyng[j]
ng[1..10]'= phing*(ninfg(vg[j])-ng[j])/taung(vg[j])
hg[1..10]'= phihg*(hinfg(vg[j])-hg[j])/tauhg(vg[j])
rg[1..10]'= phig*(rinfg(vg[j])-rg[j])/taurg
cag[1..10]'=epsg*(-icag(vg[j])-itg(vg[j],rg[j]) - kcag*cag[j])
sg[1..10]'=alphag*(1-sg[j])*sinfg(vg[j]+abg)-betag*sg[j]

@ dt=.2,total=999,meth=qualrk,xp=t,yp=v5,xlo=0,xhi=1000,ylo=-80,yhi=20.,bound=5000

done