;; 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 gate-complex) (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 gate-complex) (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)

			 (defun 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 gate-complex) (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 gate-complex) (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 gate-complex) (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 gate-complex) (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 gate-complex) (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 gate-complex) (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 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)
	   )

))