TITLE SDH Neuron - Delay Firing with AdEx 

COMMENT
AdEx + IL, Iadapt, Iadapt2, Isub, IK_A, IK_lt

Implemented by L Medlock 
ENDCOMMENT


NEURON {
  SUFFIX Delay_AdEx
  USEION k READ ek WRITE ik
  USEION na READ ena WRITE ina
  NONSPECIFIC_CURRENT il, iklt, ika, iadapt, isub, iadapt2
  RANGE gnabar, gkbar, gleak, gkltbar, gkabar, gadaptbar, gadaptslowbar, gsubbar, el, eklt, eadapt, eadaptslow, eka, esub, ina, ik, il, iklt, ika, iadapt, iadapt2, isub
  RANGE betam, gammam, phiw, betaw, gammaw, phiz, betaz, gammaz, phia, phib, betap, gammap, taup, betaq, gammaq, tauq, betapslow, gammapslow, taupslow
  GLOBAL minf, winf, tauw, zinf, tauz, ainf, taua, binf, taub, pinf, qinf, pslowinf, taub_1
  THREADSAFE
}

UNITS {
	(mV) = (millivolt)
	(mA) = (milliamp)
	(S)  = (siemens)
}

PARAMETER {
  v (mV)
  dt (ms)
  
  :Ileak
  gleak  = 0.00005 (S/cm2)
  el = -70 (mV)

  : IK_A (inactivating (A-type) potassium)
  gkabar = 0.002 (S/cm2) 
  eka = -100 (mV) 
  phia = 1.0
  phib = 1.0

  : Iadapt (AHP-like current) from Ratte et al. (2015)
  gadaptbar = 0.1 (S/cm2)        
  eadapt = -100 (mV) 
  betap = 0 (mV)
  gammap = 5 (mV)
  taup = 5 (ms)

  : INa
  gnabar = 0 (S/cm2)   
  ena = 50 (mV) 
  betam  = -1.2 (mV)
  gammam = 18 (mV)

  : IK_dr (delayed-rectifier K+ current)
  gkbar  = 0 (S/cm2)
  ek = -100 (mV)
  phiw = 0.15
  betaw = -10 (mV)
  gammaw = 10 (mV)

 : IK_lt  (low-threshold non-inactivating (Kv1-type) potassium)
  gkltbar = 0 (S/cm2) 
  eklt = -100 (mV)  
  phiz = 0.15
  betaz = -21 (mV)
  gammaz = 15 (mV)

  : Iadapt2 (slower AHP) from Ratte et al. (2015)
  gadaptslowbar = 0 (S/cm2)  
  eadaptslow = -100 (mV)       
  betapslow = 0 (mV)
  gammapslow = 5 (mV)
  taupslow = 40 (ms)

  : Isub (subthreshold inward current) from Ratte et al. (2015)
  gsubbar = 0 (S/cm2)       
  esub = 50 (mV) 
  betaq = -40 (mV)
  gammaq = 10 (mV)
  tauq = 2 (ms)

}

ASSIGNED {
  ina (mA/cm2)
	ik (mA/cm2)
  il (mA/cm2)
  iklt (mA/cm2)
  ika (mA/cm2)
  iadapt (mA/cm2)
  isub (mA/cm2)
  iadapt2 (mA/cm2)
  gk (S/cm2)  
  gna (S/cm2) 
  gklt (S/cm2)
  gka (S/cm2)
  gadapt (S/cm2)
  gsub (S/cm2)
  gadaptslow (S/cm2)
	minf
	winf
	tauw (ms)
  zinf
  tauz (ms)
  ainf
  taua (ms)
  binf
  taub (ms)
  taub_1(ms)
  pinf
  qinf
  pslowinf
}

STATE {
	m   FROM 0 TO 1
  w   FROM 0 TO 1
  z   FROM 0 TO 1
  a   FROM 0 TO 1
  b   FROM 0 TO 1
  p   FROM 0 TO 1
  q   FROM 0 TO 1
  pslow FROM 0 TO 1
}

INITIAL {
  rates(v)
  m = minf
  w = winf
  z = zinf
  a = ainf
  b = binf
  p = pinf
  q = qinf
  pslow = pslowinf
}

BREAKPOINT {
	SOLVE states METHOD cnexp
  gna = gnabar * m
  gk  = gkbar * w  
  gklt = gkltbar * z
  gka = gkabar * a^4 * b
  gadapt = gadaptbar * p
  gadaptslow = gadaptslowbar * pslow
  gsub = gsubbar * q

	ik  = gk * (v - ek)
  ina = gna * (v - ena)
  il  = gleak * (v - el)
  iklt = gklt * (v - eklt)
  ika = gka * (v - eka)
  iadapt = gadapt * (v - eadapt)
  iadapt2 = gadaptslow * (v - eadaptslow)
  isub = gsub * (v - esub)
}

DERIVATIVE states {
	rates(v)
  : Sodium Current (INa)
  m = minf
  : Delayed Rectifier K+ Current (IK_dr)
  w' = phiw*(winf-w)/tauw
  : Low-Threshold K+ Current (IK_lt)
  z' = phiz*(zinf-z)/tauz
  : A-type K+ Current (IK_A)
  a' = phia*(ainf-a)/taua
  b' = phib*(binf-b)/taub
  : Fast AHP Current (Iadapt)
  p' = (pinf-p)/taup
  : Slow AHP Current (Iadapt2)
  pslow' = (pslowinf-pslow)/taupslow
  : Inward Subthreshold Current (Isub)
  q' = (qinf-q)/tauq
}


PROCEDURE rates( v (mV) ) {
  TABLE minf, winf, tauw, zinf, tauz, ainf, taua, binf, taub_1, pinf, taub, pslowinf, qinf FROM -100 TO 100 WITH 200
  : Fast Activating Inward Current (INa):
  minf = 0.5*(1 + tanh((v-betam)/gammam))
  : Delayed Rectifier K+ Current (IK_dr)
  winf = 0.5*(1 + tanh((v-betaw)/gammaw))
  tauw = 1/(cosh((v-betaw)/(2*gammaw)))
  : Low-Threshold K+ Current (IK_lt)
  zinf = 0.5*(1 + tanh((v-betaz)/gammaz))
  tauz = 1/(cosh((v-betaz)/(2*gammaz)))
  : A-type K+ Current (IK_A)
  ainf = (1 / (1 + exp(-((v+65)/8.5))))
  taua = 1/(exp((v+35.82)/19.69)+exp(-(v+79.69)/12.7)) + 0.37
  binf = (1 / (1 + exp((v+78)/6)))
  taub_1 = 1/((exp((v+46.05)/5)+exp(-(v+238.4)/37.45)))
  if (v<-63) {
		taub = taub_1
		}
	else {
		taub = 19
		}
  : Fast AHP Current (Iadapt)
  pinf = (1 / (1 + exp((betap-v)/gammap)))
  : Slow AHP Current (Iadapt2)
  pslowinf = (1 / (1 + exp((betapslow-v)/gammapslow)))
  : Subthreshold Current (Isub)
  qinf = (1 / (1 + exp((betaq-v)/gammaq)))
 
}