# XPP code for ML-type model with HH-style subthreshold conductances
# from Ratte et al. eLife 2014; 3: e02370

# DIFFERENTIAL EQUATIONS

dv/dt = (Istim-gna*minf(V)*(V-Vna)-gk*w*(V-VK)-gl*(V-Vl)-gsubNa*m*(V-Vna)-gsubK*m*(V-Vk)-gadapt*z*(v-VK))/cap
dw/dt = phi*(winf(V)-w)/tauw(V)
dm/dt = ka*(((V-Va)/sa)/(exp((V-Va)/sa)-1))*(1-m)-kb*exp((V-Vb)/sb)*m
dz/dt = (1/(1+exp((beta_z-v)/gamma_z))-z)/tauz

#dm/dt = a*(1-m)-b*m
#a = ka*(((V-Va)/sa)/(exp((V-Va)/sa)-1))
#b = kb*exp((V-Vb)/sb)

# PARAMETERS AND FUNCTIONS

minf(v)=.5*(1+tanh((v-v1)/v2))
winf(v)=.5*(1+tanh((v-v3)/v4))
tauw(v)=1/cosh((v-v3)/(2*v4))

# Parameters for ML model
param Istim=0
param gsubNa=0,gsubK=0
# vary gsubNa and gsubK to reproduce spiking patterns show in figure 1.
param vna=50
param vk=-100,vl=-70
param gk=20,gl=2,gna=20

# v1 and v2 correspond to beta_m and gamma_m respectively
# v3 and v4 ''-----------------w-----------w-----------''   

param v1=-1.2,v2=14
param v3=-10,v4=10,phi=.15
param gadapt=0.5,tauz=300,beta_z=0,gamma_z=4
param cap=2

# Parameters for HH style conductances
param Va=-24,Vb=-24
param sa=-17,sb=-17
param ka=1,kb=1


# INITIAL CONDITIONS
V(0)=-70
w(0)=0.000025
m(0)=0

@ total=10000,dt=.05,xlo=-100,xhi=60,ylo=-.125,yhi=.6,xp=v,yp=w
@ meth=euler
@ MAXSTOR=1000000,bounds=10000

done