;; Cerebellar Purkinje Cell: resurgent Na current and high frequency
;; firing (Khaliq et al 2003).

nemo-model 

 Khaliq03

  input v (cai from ion-pools) (ica from ion-currents)

  const ena = 60
  const ek = -88
  const ca0 = 1e-4

  component (type ionic-current) (name CaBK)
   ;: BK-type Purkinje calcium-activated potassium current

              component (type gate)

			 ;; constants
			 const ztau = 1.0


			 ;; rate functions

			 CaBK_v = (v + 5)

			 minf = (let ((vh -28.9)
			              (k  6.2))
		                      (1.0 / (1.0 + exp (neg ((CaBK_v - vh) / k)))))

			 mtau = (let
	                          ((y0   0.000505)
		                   (vh1  -33.3)
		                   (k1   -10.0)
		                   (vh2  86.4)
		                   (k2   10.1))
	                          ((1e3) * (y0 + 1 / (exp ((CaBK_v + vh1) / k1) + 
	                                                    exp ((CaBK_v + vh2) / k2)))))

			 hinf = (let 
                                  ((y0  0.085)
		                   (vh  -32.0)
		                   (k   5.8))
	                          (y0 + (1 - y0) / (1 + exp ((CaBK_v - vh) / k))))


			 htau = (let
                                  ((y0   0.0019)
		                   (vh1  -54.2)
		                   (k1   -12.9)
		                   (vh2   48.5)
		                   (k2    5.2))
	                          ((1e3) * (y0 + 1 / (exp ((CaBK_v + vh1) / k1) + exp ((CaBK_v + vh2) / k2)))))


			 zinf = (let ((k 0.001))
	                          (1 / (1 + (k / cai))))

			 z_alpha = (zinf / ztau)
			 z_beta  = ((1 - zinf) / ztau)

			 reaction
			  z
                           transitions (<-> O C z_alpha z_beta)
                           conserve  (1 = (O + C))
                           initial   (let ((k 0.001)) 
                                        (1 / (1 + k / ca0)))
                           open O
                           power 2
			 
			 output z

			 
			 hh-ionic-gate 
			  CaBK  ;; ion name: exported variables will be of the form {ion}_{id}
	                   initial-m  (minf)
	                   initial-h  (hinf)
	                   m-power    3
	                   h-power    1
	                   m-inf      (minf)
	                   m-tau      (mtau)
	                   h-inf      (hinf)
	                   h-tau      (htau)
              
              component (type pore)
			 const  gbar_CaBK  = 0.007
			 output gbar_CaBK 
              
              component (type permeating-ion) (name k)
			 const e_CaBK = ek
			 output e_CaBK
              
              ;; end BK current

   
  component (type ionic-current) (name CaP)
              ;; HH P-type Calcium current
              
              component (type gate)


			 ;; rate functions
			 inf =  (let ((cv  -19)  (ck  5.5))
	                          (1.0 / (1.0 + exp (neg ((v - cv) / ck)))))

			 tau = ((1e3) *
	                        (if (v > -50) 
	                         then (0.000191 + (0.00376 * exp (neg (((v + 41.9) / 27.8) ^ 2))))
	                         else (0.00026367 + (0.1278 * exp (0.10327 * v)))))


			 hh-ionic-gate 
			  CaP  ;; ion name: exported variables will be of the form {ion}_{id}
	                   initial-m  (inf)
	                   m-power    1
	                   h-power    0
	                   m-inf      inf
	                   m-tau      tau

              component (type permeability)

			 fun ghk (v ci co)
	                   (let ((F  9.6485e4)
	                   	 (R  8.3145)
   	                         (T  (22 + 273.19))  
                                 (Z  2)
				 (E  ((1e-3) * v)))
                             (let ((k0 ((Z * (F * E)) / (R * T))))
                                (let ((k1 (exp (neg (k0))))
	                             (k2 (((Z ^ 2) * (E * (F ^ 2))) / (R * T))))
                                (1e-6) * (if (abs (1 - k1) < 1e-6) 
                                          then (Z * F * (ci - (co * k1)) * (1 - k0))
                                          else (k2 * (ci - (co * k1)) / (1 - k1))))))
			 
			 const pcabar  = 0.00005
			 const cao     = 2.4
			 pca           = (pcabar * ghk (v cai cao))

			 output pca
              
              component (type permeating-ion) (name ca)
              
              
              ;; end CaP current

  component (type ionic-current) (name K1)
   ;; HH TEA-sensitive Purkinje potassium current

              component (type gate)

			 ;; constants

			 ;; rate functions

			 K1_v = (v + 11) ;; account for junction potential

			 minf = (let ((mivh -24)
		                      (mik  15.4))
	                              (1 / (1 + exp (neg (K1_v - mivh) / mik))))


			 mtau = (let ((mty0   0.00012851)
	                              (mtvh1  100.7)
		                      (mtk1   12.9)
		                      (mtvh2  -56.0)
		                      (mtk2   -23.1))
	                             (1e3 * (if (K1_v < -35) 
				  	then (3.0 * (3.4225e-5 + 0.00498 * exp (neg (K1_v) / -28.29)))
					else (mty0 + 1.0 / (exp ((K1_v + mtvh1) / mtk1) + exp ((K1_v + mtvh2) / mtk2)))
					)))

			 hinf = (let ((hiy0  0.31)
		                      (hiA   0.78)
		                      (hivh  -5.802)
		                      (hik   11.2))
	                           (hiy0 + hiA / (1 + exp ((K1_v - hivh) / hik))))


			 htau =  (1e3 * (if ( K1_v > 0 ) 
			                   then (0.0012 + 0.0023 * exp (-0.141 * K1_v))
			                   else (1.2202e-05 + 0.012 * exp (neg (((K1_v - (-56.3)) / 49.6) ^ 2)))))

			 hh-ionic-gate 
			  K1  ;; ion name: exported variables will be of the form {ion}_{id}
	                   initial-m  (minf)
	                   initial-h  (hinf)
	                   m-power    3
	                   h-power    1
	                   m-inf      (minf)
	                   m-tau      (mtau)
	                   h-inf      (hinf)
	                   h-tau      (htau)
              
              component (type pore)
			 const  gbar  = 0.004
			 output gbar
              
              component (type permeating-ion) (name k)
			 const e = ek
			 output e
              
              ;; end K1 current

  component (type ionic-current) (name K2)
   ;; HH Low TEA-sensitive Purkinje potassium current

              component (type gate)

			 ;; constants

			 ;; rate functions

			 K2_v = (v + 11) ;; account for junction potential

			 minf = (let ((mivh -24)
		                      (mik  20.4))
	                            (1 / (1 + exp ((neg (K2_v - mivh)) / mik))))


			 mtau = ((1e3) * (if (K2_v < -20) 
					  then (0.000688 + 1 / (exp ((K2_v + 64.2) / 6.5) + exp ((K2_v - 141.5) / -34.8)))
					  else (0.00016 + 0.0008 * exp (-0.0267 * K2_v))))

			 
			 hh-ionic-gate 
			  K2  ;; ion name: exported variables will be of the form {ion}_{id}
	                   initial-m  (minf)
	                   m-power    4
	                   h-power    0
	                   m-inf      (minf)
	                   m-tau      (mtau)
              
              component (type pore)
			 const  gbar  = 0.002
			 output gbar
              
              component (type permeating-ion) (name k)
			 const e = ek
			 output e
              
              ;; end K2 current



	
  component (type ionic-current) (name K3)
   ;; HH slow TEA-insensitive Purkinje potassium current

              component (type gate)

			 ;; constants

			 ;; rate functions

			 K3_v = (v + 11) ;; account for junction potential

			 minf = (let ((mivh -16.5)
		                      (mik  18.4))
	                            (1 / (1 + exp ((neg (K3_v - mivh)) / mik))))


			 mtau = ((1e3) * (0.000796 + 1.0 / (exp ((K3_v + 73.2) / 11.7) + exp ((K3_v - 306.7) / -74.2))))
			 
			 hh-ionic-gate 
			  K3  ;; ion name: exported variables will be of the form {ion}_{id}
	                   initial-m  (minf)
	                   m-power    4
	                   h-power    0
	                   m-inf      (minf)
	                   m-tau      (mtau)

              
              component (type pore)
			 const  gbar  =  0.004
			 output gbar
              
              component (type permeating-ion) (name k)
			 const e = ek
			 output e
              ;; end K3 current

  component (type ionic-current) (name Narsg) 

              ;; constants
                     
              component (type gate)
				
                         const Con   = 0.005
                         const Coff  = 0.5
                         const Oon   = 0.75
                         const Ooff  = 0.005


                         const alfac = (pow ((Oon / Con) (1.0 / 4.0)))
                         const btfac = (pow ((Ooff / Coff) (1.0 / 4.0)))
                         
                         const alpha = 150
                         const beta  = 3
                         const gamma = 150
                         const delta = 40
                         const epsilon = 1.75
                         const zeta = 0.03
                         const x1 = 20
                         const x2 = -20
                         const x3 = 1e12
                         const x4 = -1e12
                         const x5 = 1e12
                         const x6 = -25
  
                         ;; rate functions
                         
                         f01 = (4.0 * alpha * exp (v / x1))
                         f02 = (3.0 * alpha * exp (v / x1))
                         f03 = (2.0 * alpha * exp (v / x1))
                         f04 = (alpha * exp (v / x1))
                         f0O = (gamma * exp (v / x3))
                         fip = (epsilon * exp (v / x5))
                         f11 = (4.0 * alpha * alfac * exp (v / x1))
                         f12 = (3.0 * alpha * alfac * exp (v / x1))
                         f13 = (2.0 * alpha * alfac * exp (v / x1))
                         f14 = (alpha * alfac * exp (v / x1))
                         f1n = (gamma * exp (v / x3))
                         
                         fi1 = (Con)
                         fi2 = (Con * alfac)
                         fi3 = (Con * alfac * alfac)
                         fi4 = (Con * alfac * alfac * alfac)
                         fi5 = (Con * alfac * alfac * alfac * alfac)
                         fin = (Oon)
		 
                         b01 = (beta * exp (v / x2))
                         b02 = (2.0 * beta * exp (v / x2))
                         b03 = (3.0 * beta * exp (v / x2))
                         b04 = (4.0 * beta * exp (v / x2))
                         b0O = (delta * exp (v / x4))
                         bip = (zeta * exp (v / x6))
                         
                         b11 = (beta * btfac * exp (v / x2))
                         b12 = (2.0 * beta * btfac * exp (v / x2))
                         b13 = (3.0 * beta * btfac * exp (v / x2))
                         b14 = (4.0 * beta * btfac * exp (v / x2))
                         b1n = (delta * exp (v / x4))
                         
                         bi1 = (Coff)
                         bi2 = (Coff * btfac)
                         bi3 = (Coff * btfac * btfac)
                         bi4 = (Coff * btfac * btfac * btfac)
                         bi5 = (Coff * btfac * btfac * btfac * btfac)
                         bin = (Ooff)

                         reaction
                          z
                           transitions
                            <-> C1 C2 f01 b01
                            <-> C2 C3 f02 b02
                            <-> C3 C4 f03 b03
                            <-> C4 C5 f04 b04
                            <-> C5 O  f0O b0O
                            <-> O  B  fip bip
                            <-> O  I6 fin bin
                            <-> C1 I1 fi1 bi1
                            <-> C2 I2 fi2 bi2
                            <-> C3 I3 fi3 bi3
                            <-> C4 I4 fi4 bi4
                            <-> C5 I5 fi5 bi5
                            <-> I1 I2 f11 b11
                            <-> I2 I3 f12 b12
                            <-> I3 I4 f13 b13
                            <-> I4 I5 f14 b14
                            <-> I5 I6 f1n b1n
                           conserve (1 = (I1 + I2 + I3 + I4 + I5 + I6 + C1 + C2 + C3 + C4 + C5 + O + B))
                           open O   
                           power 1
                   
                         output z


		component (type pore)
	                   const  gbar  = 0.015
	                   output gbar
		
		component (type permeating-ion) (name na)
	                   const e = ena
	                   output e

  component (type ionic-current) (name Ih)
              
              component (type gate)
                         
			 ;; rate functions
			 
			 inf = (1.0 / (1.0 + exp ((v + 90.1) / 9.9)))
			 
			 tau = ((1e3) * (0.19 + 0.72 * exp (neg (((v - (-81.5)) / 11.9) ^ 2))))

			 hh-ionic-gate 
			  Ih  ;; ion name: exported variables will be of the form {ion}_{id}
	                   initial-m (inf)
	                   m-power   1
	                   h-power   0
	                   m-inf     (inf)
	                   m-tau     (tau)

              component (type pore)
			 const  gbar  = 0.0001
			 output gbar
              
              component (type permeating-ion) (name non-specific)
			const e = -30
			output e
              
              ;; end Ih current

	
  component (type ionic-current) (name Leak)
              
              component (type pore)
			 const  gbar  = 5e-5
			 output gbar
              
              component (type permeating-ion) (name non-specific)
			 const e = -60
			 output e
              
              ;; end leak current


  component (type decaying-pool) (name ca)

              const  F = 96485.0
              const  ca_depth   = 0.1
              const  ca_beta    = 1.0
              
              d (ca) =  ((neg (ica) / (2 * ca0 * F * ca_depth)) -
                          ((if (ca < ca0) then ca0 else ca) * ca_beta))
                        ((initial ca0))

              cac = (if (ca < ca0) then ca0 else ca)

              output cac


  component (type membrane-capacitance)
           const C_m  = (1e-3)
           output C_m


  component (type voltage-clamp) (name K1)
           
           const vchold   = -71
           const vcbase   = -69
           const vcinc    = 10
           const vcsteps  = 8
           const vchdur   = 30
           const vcbdur   = 100
           
           output vchold vcbase vcinc vcsteps vchdur vcbdur

              
  component (type voltage-clamp) (name K2)
           const vchold   = -71
           const vcbase   = -69
           const vcinc    = 10
           const vcsteps  = 9
           const vchdur   = 30
           const vcbdur   = 100

           output vchold vcbase vcinc vcsteps vchdur vcbdur


  component (type voltage-clamp) (name K3)

           const vchold   = -71
           const vcbase   = -61
           const vcinc    = 10
           const vcsteps  = 8
           const vchdur   = 30
           const vcbdur   = 100

           output vchold vcbase vcinc vcsteps vchdur vcbdur


  component (type voltage-clamp) (name CaBK)

           const vchold   = -90
           const vcbase   = -40
           const vcinc    = 10
           const vcsteps  = 5
           const vchdur   = 5
           const vcbdur   = 20

           output vchold vcbase vcinc vcsteps vchdur vcbdur


  component (type voltage-clamp) (name CaP)

           const vchold   = -90
           const vcbase   = -90
           const vcinc    = 10
           const vcsteps  = 11
           const vchdur   = 5
           const vcbdur   = 10

           output vchold vcbase vcinc vcsteps vchdur vcbdur


  component (type voltage-clamp) (name Ih)

           const vchold   = -50
           const vcbase   = -60
           const vcinc    = -10
           const vcsteps  = 7
           const vchdur   = 300
           const vcbdur   = 1200

           output vchold vcbase vcinc vcsteps vchdur vcbdur



  component (type voltage-clamp) (name Leak)

           const vchold   = -71
           const vcbase   = -91
           const vcinc    = 10
           const vcsteps  = 5
           const vchdur   = 30
           const vcbdur   = 100

           output vchold vcbase vcinc vcsteps vchdur vcbdur



  component (type voltage-clamp) (name Narsg)

           const vchold   = -71
           const vcbase   = -60
           const vcinc    = 10
           const vcsteps  = 9
           const vchdur   = 30
           const vcbdur   = 100
           
           output vchold vcbase vcinc vcsteps vchdur vcbdur