# hh200x50-LEV01.ode
# Ref: Huang et al., 2009, J Physiol Pharmacol (in press)
# 200 e and 50 i HH equations
# random applied current, random conductances
# to get it started, just set the excitatory synapses
# to some random values between 0 and 1
# you will get persistent activity.
# here are the HH functions
# slowly inactivating K (Kv3.1-encoded) current was included

am(v)=0.1*(v+40)/(1-exp(-(v+40)/10))
bm(v)=4*exp(-(v+65)/18)
ah(v)=.07*exp(-(v+65)/20)
bh(v)=1/(1+exp(-(v+35)/10))
an(v)=0.01*(v+55)/(1-exp(-(v+55)/10))
bn(v)=0.25*exp(-(v+65)/80)
an3(v)=(1/Kn)*(-0.021)*(v+8.3+eta)/(exp(-(v+8.3+eta)/9.8)-1)
bn3(v)=(1/Kn)*0.0002*exp(-(V+23.6+eta)/20.7)

# Stimulus protocol
ff(t)=heav(t-t_on)*heav(t_off-t)
par t_on=20, t_off=200
par Kn=1
% gk3: 7-->4 or 1
par gk3=7
% eta: 0-->10 mV
par eta=0

# this is the current for each cell
ihh(v,m,h,n,n3,hk)=gna*h*(v-vna)*m^3+gk*(v-vk)*n^4+gk3*(v-vk)*n3^4*hk+gl*(v-vl)
# synaptic onset parameters
# s' = a(vpre)(1-s)-s/tau
ae(x)=ae0/(1+exp(-x/5))
ai(x)=ai0/(1+exp(-x/5))
par ae0=4, ai0=1

# dont recompute the random tables every time a parameter is changed
@ autoeval=0

# random synapses - 20 % connectivity
table wee % 40000 0 39999 ran(1)<.02
table wei % 10000 0 9999 ran(1)<.02
table wie % 10000 0 9999 ran(1)<.02
table wii % 2500 0 2499  ran(1)<.02
# multiply synapses by weights
special see=mmult(200,200,wee,se0)
special sei=mmult(200,50,wei,se0)
special sie=mmult(50,200,wie,si0)
special sii=mmult(50,50,wii,si0)
# random currents applied to each cell
table r_e % 200 0 199  ran(1)-.5
table r_i % 50 0 49 ran(1)-.5

# parameters
par taue=4, taui=10
par vna=50,  vk=-80,  vl=-49,  gna=40,  gk=5, gl=0.03
par ie0=6.5, ie1=0
par ii0=0, ii1=0
par gee=0.1, gie=0.1, gii=0.1, gei=0.1
par eex=0, ein=-75

# finally the ODEs
ve[0..199]'=ie0*ff(t)+ie1*r_e([j])-ihh(ve[j],me[j],he[j],ne[j],n3e[j],hke[j])-gee*see([j])*(ve[j]-eex)-gie*sie([j])*(ve[j]-ein)
vi[0..49]'=ii0+ii1*r_i([j])-ihh(vi[j],mi[j],hi[j],ni[j],n3i[j],hki[j])-gei*sei([j])*(vi[j]-eex)-gii*sii([j])*(vi[j]-ein)
# synapses...
se[0..199]'=-se[j]/taue+ae(ve[j])*(1-se[j])
si[0..49]'=-si[j]/taui+ai(vi[j])*(1-si[j])
# gating variables
me[0..199]'=am(ve[j])*(1-me[j])-bm(ve[j])*me[j]
he[0..199]'=ah(ve[j])*(1-he[j])-bh(ve[j])*he[j]
ne[0..199]'=an(ve[j])*(1-ne[j])-bn(ve[j])*ne[j]
hke[0..199]'=0.0005*(1-hke[j])-0.0014*n3e[j]^4*hke[j]
n3e[0..199]'=an3(ve[j])*(1-n3e[j])-bn3(ve[j])*n3e[j]
mi[0..49]'=am(vi[j])*(1-mi[j])-bm(vi[j])*mi[j]
hi[0..49]'=ah(vi[j])*(1-hi[j])-bh(vi[j])*hi[j]
ni[0..49]'=an(vi[j])*(1-ni[j])-bn(vi[j])*ni[j]
hki[0..49]'=0.0005*(1-hki[j])-0.0014*n3i[j]^4*hki[j]
n3i[0..49]'=an3(vi[j])*(1-n3i[j])-bn3(vi[j])*n3i[j]

# initial data
init ve[0..199]=-75
init vi[0..49]=-75
# some numerical settings
@ total=200, meth=euler, nout=10, dt=.01
@ xlo=0, xhi=200, yhi=40, ylo=-90
done