;; 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
		
		;; end Narsg component





   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