TITLE SDH Neuron - Single Firing with AdEx 

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

Implemented by L Medlock 
ENDCOMMENT


NEURON {
  SUFFIX Single_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)  : same as delay
  el = -70 (mV)

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

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

  : Iadapt (AHP-like current) from Ratte et al. (2015)
  gadaptbar = 0.04 (S/cm2)        : same as tonic
  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_A (inactivating (A-type) potassium)
  gkabar = 0 (S/cm2)  : Changes per neuron type
  eka = -100 (mV) 
  phia = 1.0
  phib = 1.0

  : 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
  : Adapting (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)))
}