TITLE Propagated Hodgkin-Huxley Action Potential with EC Potentials
:  Simulation Resources, Inc.  1/31/93, rev 3/22/2001
:  Reference - (Membrane model) A.L. Hodgkin and A.F. Huxley,
:    J. Physiol. (Lond.) 117:500 (1952).
:    (Cable solution) J.W. Moore, F. Ramon, R.W. Joyner,
:    Biophys. J. 15:11 (1975).
:  Default output - Extracellular potentials at four locations
:    from the propagated action potential, plus potential from
:    bipolar electrode at one location. See HH3.MOD for a
:    description of the action potential model.

DEFINE segments  100

PARAMETER
{
  radius = .025		(cm)
  Ri = .05		(kOhm.cm)
  Ro = .01		(kOhm.cm)
  Cm = 1.		(uF/cm2)
  Temp = 279.3		(K)
  q10g = 3.0			: q10 for gating
  q10c = 1.3			: q10 for ion current
  GbarNa = 120. (mS/cm^2)
  GbarK = 36.	(mS/cm^2)
  Gl = 0.3	(mS/cm^2)
  Nai = 44.5	(mM)
  Nao = 440	(mM)
  Ki = 200.8	(mM)
  Ko = 10	(mM)
  El = -49.387	(mV)
  Istim = 5.		(uA)
  Istim2 = 25	(uA/cm^2)
  Tdelay = 0.0	(msec)
  Tdelay2 = 15	(msec)
  Tdur = 0.5	(msec)
  delta_time = .01	(msec)
  delta_x = .05		(cm)
  electrode_pos = 4.0	(cm)
  electrode_spacing = 0.05 (cm)
}

ASSIGNED
{
  mesh_ratio
  x[segments] FROM 0 TO 2.5	(cm)
  f[segments]  g[segments]
  Vec[segments] FROM -1 TO 4	(mV)
  Vec_measured FROM -1 TO 1	(mV)
  electrode_indx_1
  electrode_indx_2
  g_mult c_mult
  am ah an bm bh bn
  ENa FROM 0 TO 60	(mV)
  EK FROM -100 TO -50	(mV)
  istim FROM 0 TO 200	(uA/cm^2)
}

INDEPENDENT
{
  time FROM 0 TO 6 WITH 300  (msec)
}

STATE
{
  Vrest FROM -100 TO 60 START -60 (mV)
  V[segments] FROM -80 TO 60	  (mV)
  m[segments] FROM 0 TO 1
  n[segments] FROM 0 TO 1
  h[segments] FROM 0 TO 1
}

PLOT Vec[10],Vec[20],Vec[40],Vec[80],Vec_measured VS time

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

  loc_v = loc_v + 60.0
  s = 2.5 - .1 * loc_v
  IF (fabs(s) >= .001)
  {
      am = s / (exp(s) - 1.) * g_mult
  }
  ELSE
  {
      am = 1. * g_mult
  }
  bm = 4. * exp(-loc_v / 18.) * g_mult
  ah = .07 * exp(-loc_v / 20.) * g_mult
  bh = 1. / (exp(3. - .1 * loc_v) + 1.) * g_mult
  s = 1. - .1 * loc_v
  IF (fabs(s) >= .001)
  {
      an = .1 * s / (exp(s) - 1.) * g_mult
  }
  ELSE
  {
      an = .1 * g_mult
  }
  bn = .125 * exp(-loc_v / 80.) * g_mult
}

NONLINEAR restingV
{
  LOCAL fm,fm3,fh,fn,GNa,GK

  hhco(Vrest)
  fm = am / (am + bm)
  fm3 = fm * fm * fm
  fh = ah / (ah + bh)
  fn = an / (an + bn)
  GNa = GbarNa * fm3 * fh
  GK = GbarK * fn * fn * fn * fn

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

INITIAL
{
  mesh_ratio = delta_time * radius / (2. * Ri * Cm * delta_x * delta_x)
  g_mult = q10g^(0.1 * (Temp - 279.3))	: Temp effect on gating process
  c_mult = q10c^(0.1 * (Temp - 279.3))	: Temp effect on ionic current
  ENa = .086114 * Temp * log(Nao / Nai)
  EK = .086114 * Temp * log(Ko / Ki)
  electrode_indx_1 = electrode_pos / delta_x + delta_x * 0.5
  electrode_indx_2 = (electrode_pos + electrode_spacing) / delta_x + delta_x * 0.5

  SOLVE restingV
  VERBATIM
    switch (error)
    {
      case 1:
	puts("Error: Jacobian matrix is singular");
	break;
      case 2:
	puts("Error: Newton method did not converge");
    }
  ENDVERBATIM

  hhco(Vrest)
  FROM indx = 0 TO segments - 1
  {
    m[indx] = am / (am + bm)
    h[indx] = ah / (ah + bh)
    n[indx] = an / (an + bn)
    V[indx] = Vrest
    Vec[indx+1] = 0.
    x[indx] = indx * delta_x
    f[indx] = radius / (2. * Ri * Cm)
  }
}

DERIVATIVE channels
{
    FROM indx = 0 TO segments - 1
    {
    hhco(V[indx])
    m'[indx] = am * (1. - m[indx]) - bm * m[indx]
    h'[indx] = ah * (1. - h[indx]) - bh * h[indx]
    n'[indx] = an * (1. - n[indx]) - bn * n[indx]
    }
}

PARTIAL cable
{
  LOCAL GNa,GK,INa,IK,m3

  SOLVE channels METHOD heun
  istim = pulse(time, Tdelay, Istim, Tdur)
	   + pulse(time, Tdelay2, Istim2, Tdur)
  FROM indx = 0 TO segments - 1
  {
    m3 = m[indx] * m[indx] * m[indx]
    GNa = GbarNa * m3 * h[indx]
    INa = GNa * (V[indx] - ENa)
    GK = GbarK * n[indx] * n[indx] * n[indx] * n[indx]
    IK = GK * (V[indx] - EK)
    g[indx] = - (INa + IK + Gl * (V[indx] - El)) * c_mult / Cm
  }
~ V' = f * DEL2(V) + g
~ DEL V[FIRST] = - Ri * istim / (3.14159 * radius * radius)
~ DEL V[LAST]  = 0.
}

BREAKPOINT
{
  LOCAL indx
  SOLVE cable
  FROM indx = 1 TO segments
  {
    Vec[indx] = (V[indx-1] - V[indx]) * Ro / Ri
  }
  Vec_measured = (Vec[electrode_indx_1] - Vec[electrode_indx_2])
}