TITLE Ebihara-Johnson Model for Sodium Current
: National Biomedical Simulation Resource   11-89
:  Reference - L. Ebihara and E.A. Johnson
:    Biophys. J. 32:779-790 (1980). This program generates an action
:    potential depolarization with the switch variable VclampOn=0 and
:    with VclampOn=1, V is fixed to simulate a voltage clamp experiment.
:    A potassium current with fixed conductance has been added to give
:    a stable rest potential and partial repolarization.
:  Default output - Depolarization of action potential from chick
:    embryonic heart cells.  V starts at about -72 mV and
:    rises linearly to about -42 mV during the stimulus. Then an s-shaped
:    depolarization begins with V rising to a peak of +14.5 mV at 1.60 msec.
:    Finally, there is a slow repolarization, reaching +4.38 mV at 5 msec.

PARAMETER
{
  Temp = 310	(K)
  q10g = 3.0			: q10 for gating
  q10c = 1.3			: q10 for ion current
  Cm = 1.3	(uF/cm^2)

  GbarNa = 23.0 (mS/cm^2)
  GK = .05	(mS/cm^2)
  ENa = 29.0	(mV)
  EK = -72.0	(mV)

  Istim = 80	(uA/cm^2)
  Tdelay = 0.0	(msec)
  Tdelay2 = 15	(msec)
  Tdur = 0.5	(msec)

  VclampOn = 0
  Vclamp1 = -20	(mV)
  Vclamp2 = -50	(mV)
}

INDEPENDENT
{
  time FROM 0 TO 5 WITH 100 (msec)
}

PLOT V VS time

STATE
{
  Vrest FROM -100 TO 60 START -72 (mV)
  V FROM -80 TO 60	(mv)
  m FROM 0 TO 1
  h FROM 0 TO 1
}

ASSIGNED
{
  g_mult c_mult
  am bm ah bh m3

  GNa FROM 0 TO 40	(mS/cm^2)
  INa FROM -500 TO 100	(uA/cm^2)
  Iion FROM -500 TO 100 (uA/cm^2)
  istim FROM 0 TO 200	(uA/cm^2)
}

PROCEDURE ejco(loc_v)
: Subroutine for alphas & betas
{
  LOCAL s

  s = loc_v + 47.13
  IF (fabs(s) >= .001)
  {
      am = .32 * s / (1. - exp(-s)) * g_mult
  }
  ELSE
  {
      am = .32 * g_mult
  }
  bm = .08 * exp(-loc_v / 11.) * g_mult
  IF ( loc_v >= -40 )
  {
      ah = 0.
      bh = 1. / (0.13 * (exp(- (loc_v + 10.66)/11.1) + 1.)) * g_mult
  }
  ELSE
  {
      ah = 0.135 * exp(-(80 + loc_v)/6.8) * g_mult
      bh = (3.56 * exp(0.079 * loc_v) + 3.1e5 * exp(0.35 * loc_v) ) * g_mult
  }
}

PROCEDURE currents (Vloc)
{
  m3 = m * m * m
  GNa = GbarNa * m3 * h
  INa = GNa * (Vloc - ENa) * c_mult
  Iion = INa + GK * (Vloc - EK) * c_mult
}

NONLINEAR restingV
{

  ejco(Vrest)
  m = am / (am + bm)
  h = ah / (ah + bh)
  currents(Vrest)

~ 0 = GNa * (Vrest - ENa) + GK * (Vrest - EK)
}

INITIAL
{
  g_mult = q10g^(0.1 * (Temp - 310))	: Temp effect on gating process
  c_mult = q10c^(0.1 * (Temp - 310))	: Temp effect on ionic current

  : Find resting potential
  SOLVE restingV
  : error = 0 if OK, = 1 if singular or ill-conditioned matrix, = 2 if the
  : newton method did not converge
  VERBATIM
    switch (error)
    {
      case 1:
	puts("Error: Jacobian matrix is singular");
	break;
      case 2:
	puts("Error: Newton method did not converge");
    }
  ENDVERBATIM

  V = Vrest
  ejco(V)
  m = am / (am + bm)
  h = ah / (ah + bh)
  currents(V)

}

DERIVATIVE voltageclamp
{
    V = Vrest
	+ threshold(time, Tdelay, "min")*(Vclamp1 - Vrest)
	+ threshold(time, Tdelay2, "min")*(Vclamp2 - Vclamp1)
    ejco(V)
    m' = am * (1. - m) - bm * m
    h' = ah * (1. - h) - bh * h
    currents(V)
}

DERIVATIVE currentclamp
{
    istim = pulse(time, Tdelay, Istim, Tdur)
    currents(V)
    V' = (istim - Iion) / Cm
    ejco(V)
    m' = am * (1. - m) - bm * m
    h' = ah * (1. - h) - bh * h
}

BREAKPOINT
{
    IF (VclampOn != 0.0)
    {
	SOLVE voltageclamp
    }
    ELSE 
    {
	SOLVE currentclamp
    }
}