TITLE DiFrancesco-Noble Model for Purkinje Cell
: National Biomedical Simulation Resource, 8-2-90
: Reference - D. DiFrancesco and D. Noble, Phil. Trans. R. Soc. Lond.
:    B 307:353-398, 1985.
: Parameter values are for Purkinje tissue.
: The system of units used is as follows - nA uS sec uF mV.
: Default output - action potential with depolarization at time=0
:    (see Fig. 6 of DiF-N paper).  Note that there are some differences
:    between the output of this program and the figures in the paper,
:    probably due to differences in the integration method.
:  $Id: difrnobl.mod,v 1.2 90/11/01 15:13:40 mdf Exp $

LOCAL alpha, beta

PARAMETER
{
  R = 8.3143                    (J/K/mol)
  T = 307.73                    (K)
  F = 96.487                    (Coul/mmol)
  Cm = 12                       (uF/cm2)
  radius = 0.005                (cm)
  length = 0.2                  (cm)
  fVecs = .1
  Istim =  250.                 (nA)    :Parameters for current stimulus
  Idelay = 0.                   (sec)
  Idur = .005                   (sec)
  vclamp_mode = 0.                      :Set to 1 for voltage clamp mode
  Vdelay = .05                  (sec)
  Vclamp = 90.                  (mV)    :Size of step relative to Em0
  Nao = 140.                    (mM)
  Cao = 2.                      (mM)
  Kb = 4.                       (mM)
  Kmf = 45.                     (mM)
  gfNa = 3.                     (uS)
  gfK = 3.                      (uS)
  iKmax = 180.                  (nA)
  gK1 = 920.                    (uS)
  Km1 = 210.                    (mM)
  ito_scale = 1.                        :Permits turning off of ito (Fig.5)
  Kmto = 10.                    (mM)
  Kact4 = .0005                 (mM)
  gbNa = .18                    (uS)
  gbCa = .02                    (uS)
  ipmax = 125.                  (nA)
  KmK = 1.                      (mM)
  KmNa = 40.                    (mM)
  KNaCa = .02                   (mM)
  nNaCa = 3.
  dNaCa = .001
  gamma = .5
  gNa = 750.                    (uS)
  Psi = 15.
  Kmf2 = .001                   (mM)
  Caup_bar = 5.                 (mM)
  KmCa = .001                   (mM)
  rr = 2.
  tau_rel = .05                 (sec)
  tau_rep = 2.                  (sec)
  tau_up = .025                 (sec)
  P = .7                        (1/sec)

}

ASSIGNED
{
  RTF
  FRT
  Vup                           (cm3)
  Vrel                          (cm3)
  Vi                            (cm3)
  Ve                            (cm3)
  surface_area                  (cm2)
  itot FROM -100 TO 200         (nA)
  C                             (uF)
  i_f FROM -50 TO 50            (nA)
  EK
  ENa
  ECa
  ENaCa
  Emh
  iK FROM -50 TO 50             (nA)
  iK1 FROM -50 TO 50            (nA)
  ito FROM -50 TO 50            (nA)
  ibNa FROM -50 TO 50           (nA)
  ibCa FROM -50 TO 50           (nA)
  ip FROM -50 TO 50             (nA)
  iNaCa FROM -50 TO 50          (nA)
  iNa
  isiCa FROM -50 TO 50          (nA)
  isiK FROM -50 TO 50           (nA)
  isiNa FROM -50 TO 50          (nA)
  isi FROM -50 TO 50            (nA)
  iCa FROM -50 TO 50            (nA)
  irel FROM 0 TO 1000           (nA)
  itr FROM 0 TO 1000            (nA)
  iup FROM 0 TO 1000            (nA)
  imK
  ipulse
}

STATE
{
  Em FROM -100 TO 50 START -75  (mV)
  y  FROM 0 TO 1
  x  FROM 0 TO 1
  r  FROM 0 TO 1
  m  FROM 0 TO 1
  h  FROM 0 TO 1
  d  FROM 0 TO 1
  f  FROM 0 TO 1
  f2 FROM 0 TO 1
  pp FROM 0 TO 1 START 1
  Nai FROM 0 TO 25 START 8      (mM)
  Cai FROM 0 TO .01 START .00005 (mM)
  Caup FROM 0 TO 10 START 2     (mM)
  Carel FROM 0 TO 10 START 1    (mM)
  Kc FROM 0 TO 10 START 4       (mM)
  Ki FROM 0 TO 200 START 140    (mM)
}

INDEPENDENT
{
  time FROM 0 TO 2.5 WITH 250   (sec)
}

INITIAL
{
  RTF = R * T / F
  FRT = 1. / RTF
  surface_area = 2. * 3.14159 * radius * length
  Vi = (1. - fVecs) * 3.14159 * radius * radius * length
  Ve = Vi * fVecs / (1. - fVecs)
  Vup = .05 * Vi
  Vrel = .02 * Vi
  C = Cm * surface_area 
:  ycoef(Em)                            :The constant initial values shown
:  y = alpha / (alpha + beta)           :are those used by DiF-N.  The
  y = 0.2                               :other initialization method (in
:  xcoef(Em)                            :comments) would permit correct
:  x = alpha / (alpha + beta)           :initialization at other voltages.
  x = .01
  rcoef(Em)
  r = alpha / (alpha + beta)
:  mcoef(Em)
:  m = alpha / (alpha + beta)
  m = .01
:  hcoef(Em)
:  h = alpha / (alpha + beta)
  h = 0.8
:  dcoef(Em)
:  d = alpha / (alpha + beta)
  d = .005
:  fcoef(Em)
:  f = alpha / (alpha + beta)
  f = 1.
:  f2 = 1. - Cai / (Cai + Kmf2)
  f2 = 1.
}

DERIVATIVE rates
{
  LOCAL temp

  EK = RTF * log(Kc / Ki)
  ENa = RTF * log(Nao / Nai)
  ECa = RTF * log(Cao / Cai) / 2.
  ENaCa = (nNaCa * ENa - 2. * ECa) / (nNaCa - 2.)
  Emh = RTF * log((Nao + .12 * Kc) / (Nai + .12 * Ki))

  ycoef(Em)             : Hyperpolarizing-activated current 
  y' = alpha * ( 1. - y ) - beta * y
  i_f = y * (Kc / (Kc + Kmf)) * (gfK * (Em - EK) + gfNa * (Em - ENa))

  xcoef(Em)             : Delayed rectifier current 
  x' = alpha * ( 1. - x ) - beta * x
  iK = x * iKmax * (Ki - Kc * exp(- Em * FRT)) / 140.

			: Background K current 
  iK1 = gK1 * (Kc / (Kc + Km1)) * ((Em - EK) / (1. + exp(2. * FRT *
		(Em - EK + 10.))))

  rcoef(Em)             : Transient outward current 
  r' = alpha * (1. - r) - beta * r
  IF (fabs(Em + 10.) < .001)
	{temp = 5. }
  ELSE
	{temp = (Em + 10.) / (1. - exp(-.2 * (Em + 10.))) }
  ito = ito_scale * r * .28 * ((.2 + Kc) / (Kmto + Kc)) * (Cai / (Kact4 + Cai))
	* temp * (Ki * exp(.02 * Em) - Kc * exp(-.02 * Em))

			: Background sodium current 
  ibNa = gbNa * (Em - ENa)

			: Resting Ca leak 
  ibCa = gbCa * (Em - ECa)

			: Na,K pump current 
  ip = ipmax * (Kc / (KmK + Kc)) * (Nai / (KmNa + Nai))

			: Na,Ca exchange current 
  iNaCa = KNaCa * (exp(gamma * (nNaCa - 2.) * Em * FRT) * (Nai^nNaCa) * Cao
	- exp(- (1. - gamma) * (nNaCa - 2.) * Em * FRT) * (Nao^nNaCa) * Cai) /
	(1. + dNaCa * (Cai * (Nao^nNaCa) + Cao * (Nai^nNaCa)))

			: Fast sodium current 
  mcoef(Em)
  m' = alpha * (1. - m) - beta * m
  hcoef(Em)
  h' = alpha * (1. - h) - beta * h
  iNa = gNa * m * m * m * h * (Em - Emh)

			: Slow inward currents 
  dcoef(Em)
  d' = alpha * (1. - d) - beta * d
  fcoef(Em)
  f' = alpha * (1. - f) - beta * f
  alpha = alpha * 0.1
  beta = beta * 0.1
  pp' = alpha * (1. - pp) - beta * pp
  f2coef(Em)
  f2' = alpha * (1. - f2) - beta * Cai * f2
  temp = Em - 50.
  if (fabs(temp) < .001) {
	isiCa = 2. * Psi
	isiK  = .01 * Psi
	isiNa = .01 * Psi
  }
  else {
	isiCa = 4. * Psi * d * f * f2 * temp * FRT / (1. - exp(- temp * 2. * FRT))
	isiK = .01 * Psi * d * f * f2 * temp * FRT / (1. - exp(- temp * FRT))
	isiNa = .01 * Psi * d * f * f2 * temp * FRT / (1. - exp(- temp * FRT))
  }
  isiCa = isiCa * (Cai * exp(100. * FRT) - Cao * exp(- 2. * temp * FRT))
  isiK = isiK * (Ki * exp(50. * FRT) - Kc * exp(- temp * FRT))
  isiNa = isiNa * (Nai * exp(50. * FRT) - Nao * exp(- temp * FRT))
  iCa = (isiCa + isiK + isiNa)
  isi = iCa + iNaCa

			: Internal sodium concentration 
  Nai' = - (iNa + ibNa + i_f/2. + 3 * ip + (nNaCa/(nNaCa - 2)) * iNaCa)
		* 1.e-6 / (Vi * F)

			: Internal calcium concentration 
  Caupcoef(Em)
  iup = alpha * Cai * (Caup_bar - Caup) - beta * Caup
  Catrcoef(Em)
  itr = alpha * pp * (Caup - Carel)
  Carelcoef(Em)
  irel = alpha * Carel * (Cai ^ rr) / (Cai^rr + KmCa^rr)
  Caup' = (iup - itr) / (4.e6 * Vup * F)
  Carel' = (itr - irel) / (4.e6 * Vrel * F)
  Cai' = - (isiCa + ibCa - 2 * iNaCa / (nNaCa - 2.) + iup - irel) /
		(2.e6 * Vi * F)

			: Cleft potassium concentration 
			: Uniform Ki; equation 63       
  imK = iK1 + iK + i_f/2. - 2. * ip
  Kc' = - P * (Kc - Kb) + imK / (1.e6 * Ve * F)

			: Internal potassium concentration 
  Ki' = - imK / (1.e6 * Ve * F)

			: Membrane current and voltage
  ipulse = pulse(time, Idelay , Istim, Idur)
  itot = i_f + iK + iK1 + ito + ibNa + ibCa + ip + iNaCa + iNa + iCa 
  IF (vclamp_mode == 1)
	{Em' = 0
	Em = Em0 + step(time, Vdelay, Vclamp) }
  ELSE
	{Em' = (ipulse - itot) / C }
}

BREAKPOINT
{
  SOLVE rates METHOD clsoda
}

PROCEDURE ycoef(V)
{
  alpha = .05 * exp(- .067 * (V + 52))
  IF ( fabs(V+52.) < .001 )
	{beta = 5. }
  ELSE
	{beta = (V + 52) / (1. - exp(-.2 * (V + 52))) }
}

PROCEDURE xcoef(V)
{
  alpha = .5 * exp(.0826 * (V + 50.)) / (1. + exp(.057 * (V + 50.)))
  beta = 1.3 * exp(-.06 * (V + 20.)) / (1. + exp(-.04 * (V + 20.)))
}

PROCEDURE rcoef(V)
{
  alpha = .033 * exp(-V / 17.)
  beta = 33. / (1. + exp( - (V + 10.) / 8.))
}

PROCEDURE mcoef(V)
{
  IF (fabs(V + 41.) < .001)
	{alpha = 2000. }
  ELSE
	{alpha = 200. * (V + 41.) / (1. - exp(- .1 * (V + 41.))) }
  beta = 8000. * exp( - .056 * (V + 66.))
}

PROCEDURE hcoef(V)
{
  alpha = 20. * exp(-.125 * (V + 75.))
  beta = 2000. / (1. + 320. * exp(-.1 * (V + 75.)))
}

PROCEDURE dcoef(V)
{
  LOCAL tmp

  tmp = V + 24.
  IF (fabs(tmp) < .001) {
	alpha = 120.
	beta = 120.
  }
  ELSE {
	alpha = 30. * tmp / (1. - exp(- tmp / 4.))
	beta = 12. * tmp / (exp(tmp / 10.) - 1.)
  }
}

PROCEDURE fcoef(V)
{
  LOCAL tmp

  tmp = V + 34.
  IF (fabs(tmp) < .001)
	{alpha = 25. }
  ELSE
	{alpha = 6.25 * tmp / (exp( tmp / 4.) - 1.) }
  beta = 50. / (1. + exp(- tmp / 4.))
}

PROCEDURE f2coef(V)
{
  alpha = .01
  beta = alpha / Kmf2
}

PROCEDURE Caupcoef(V)
{
  alpha = 2.e6 * F * Vi / (tau_up * Caup_bar)
  beta = 0.
}

PROCEDURE Catrcoef(V)
{
  alpha = 2.e6 * F * Vrel / tau_rep
  beta = 0.
}

PROCEDURE Carelcoef(V)
{
  alpha = 2.e6 * F * Vrel / tau_rel
}