;; 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 ion-channel) (name CaBK)
   ;: BK-type Purkinje calcium-activated potassium current

	      (component (type gate)

			 ;; constants
			 (const CaBK_ztau = 1.0)


			 ;; rate functions

			 (CaBK_v = (v + 5))

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

			 (CaBK_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))))))

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


			 (CaBK_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))))))


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

			 (CaBK_z_alpha = (CaBK_zinf / CaBK_ztau))
			 (CaBK_z_beta  = ((1 - CaBK_zinf) / CaBK_ztau))

			 (reaction
			  (CaBK_z
			   (transitions (<-> O C CaBK_z_alpha CaBK_z_beta))
                           (conserve  (1 = (O + C)))					     
			   (initial   (let ((k 0.001))
					(1 / (1 + k / ca0))))
			   (open O)   (power 2)))
			 
			 (output CaBK_z )  

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

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

			 (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)))))))
			 
			 ;; rate functions
			 (CaP_inf =  
			   (let ((cv  -19)  (ck  5.5))
			     (1.0 / (1.0 + exp (neg ((v - cv) / ck))))))

			 (CaP_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  (CaP_inf))
			   (m-power    1)
			   (h-power    0)
			   (m-inf      CaP_inf)
			   (m-tau      CaP_tau)))
			 
			 )

	      (component (type pore)
			 (const gmax_CaP  = 0.00005)
			 (const cao       = 2.4)
			 (gbar_CaP        = (gmax_CaP * ghk (v cai cao)))
			 (output gbar_CaP ))
	      
	      (component (type accumulating-substance) (name ca) )
	      
	      
	      ) ;; end CaP current

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

	      (component (type gate)

			 ;; constants

			 ;; rate functions

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

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


			 (K1_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)))
					))))

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


			 (K1_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  (K1_minf))
			   (initial-h  (K1_hinf))
			   (m-power    3)
			   (h-power    1)
			   (m-inf      (K1_minf))
			   (m-tau      (K1_mtau))
			   (h-inf      (K1_hinf))
			   (h-tau      (K1_htau))))
			 
			 )
	      
	      (component (type pore)
			 (const  gbar_K1  = 0.004)
			 (output gbar_K1 ))
	      
	      (component (type permeating-substance) (name k)
			 (const e_K1 = ek)
			 (output e_K1 ))
	      
	      ) ;; end K1 current

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

	      (component (type gate)

			 ;; constants

			 ;; rate functions

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

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


			 (K2_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  (K2_minf))
			   (m-power    4)
			   (h-power    0)
			   (m-inf      (K2_minf))
			   (m-tau      (K2_mtau))))
			 
			 )
	      
	      (component (type pore)
			 (const  gbar_K2  = 0.002)
			 (output gbar_K2 ))
	      
	      (component (type permeating-substance) (name k)
			 (const e_K2 = ek)
			 (output e_K2 ))
	      
	      ) ;; end K2 current



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

	      (component (type gate)

			 ;; constants

			 ;; rate functions

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

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


			 (K3_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  (K3_minf))
			   (m-power    4)
			   (h-power    0)
			   (m-inf      (K3_minf))
			   (m-tau      (K3_mtau))))
			 
			 )
	      
	      (component (type pore)
			 (const  gbar_K3  =  0.004)
			 (output gbar_K3 ))
	      
	      (component (type permeating-substance) (name k)
			 (const e_K3 = ek)
			 (output e_K3 ))
	      
	      ) ;; end K3 current

   (component (type ion-channel) (name Narsg) 

	      ;; constants
		     
	      (component (type gate)
				
			 (const Na_Con   = 0.005)
			 (const Na_Coff  = 0.5)
			 (const Na_Oon   = 0.75)
			 (const Na_Ooff  = 0.005)


			 (const Na_alfac = (pow ((Na_Oon / Na_Con) (1.0 / 4.0))))
			 (const Na_btfac = (pow ((Na_Ooff / Na_Coff) (1.0 / 4.0))))
			 
			 (const Na_alpha = 150)
			 (const Na_beta  = 3)
			 (const Na_gamma = 150)
			 (const Na_delta = 40)
			 (const Na_epsilon = 1.75)
			 (const Na_zeta = 0.03)
			 (const Na_x1 = 20)
			 (const Na_x2 = -20)
			 (const Na_x3 = 1e12)
			 (const Na_x4 = -1e12)
			 (const Na_x5 = 1e12)
			 (const Na_x6 = -25)
   
			 ;; rate functions
			 
			 (f01 = (4.0 * Na_alpha * exp (v / Na_x1)))
			 (f02 = (3.0 * Na_alpha * exp (v / Na_x1)))
			 (f03 = (2.0 * Na_alpha * exp (v / Na_x1)))
			 (f04 = (Na_alpha * exp (v / Na_x1)))
			 (f0O = (Na_gamma * exp (v / Na_x3)))
			 (fip = (Na_epsilon * exp (v / Na_x5)))
			 (f11 = (4.0 * Na_alpha * Na_alfac * exp (v / Na_x1)))
			 (f12 = (3.0 * Na_alpha * Na_alfac * exp (v / Na_x1)))
			 (f13 = (2.0 * Na_alpha * Na_alfac * exp (v / Na_x1)))
			 (f14 = (Na_alpha * Na_alfac * exp (v / Na_x1)))
			 (f1n = (Na_gamma * exp (v / Na_x3)))
			 
			 (fi1 = (Na_Con))
			 (fi2 = (Na_Con * Na_alfac))
			 (fi3 = (Na_Con * Na_alfac * Na_alfac))
			 (fi4 = (Na_Con * Na_alfac * Na_alfac * Na_alfac))
			 (fi5 = (Na_Con * Na_alfac * Na_alfac * Na_alfac * Na_alfac))
			 (fin = (Na_Oon))
			 
			 (b01 = (Na_beta * exp (v / Na_x2)))
			 (b02 = (2.0 * Na_beta * exp (v / Na_x2)))
			 (b03 = (3.0 * Na_beta * exp (v / Na_x2)))
			 (b04 = (4.0 * Na_beta * exp (v / Na_x2)))
			 (b0O = (Na_delta * exp (v / Na_x4)))
			 (bip = (Na_zeta * exp (v / Na_x6)))
			 
			 (b11 = (Na_beta * Na_btfac * exp (v / Na_x2)))
			 (b12 = (2.0 * Na_beta * Na_btfac * exp (v / Na_x2)))
			 (b13 = (3.0 * Na_beta * Na_btfac * exp (v / Na_x2)))
			 (b14 = (4.0 * Na_beta * Na_btfac * exp (v / Na_x2)))
			 (b1n = (Na_delta * exp (v / Na_x4)))
			 
			 (bi1 = (Na_Coff))
			 (bi2 = (Na_Coff * Na_btfac))
			 (bi3 = (Na_Coff * Na_btfac * Na_btfac))
			 (bi4 = (Na_Coff * Na_btfac * Na_btfac * Na_btfac))
			 (bi5 = (Na_Coff * Na_btfac * Na_btfac * Na_btfac * Na_btfac))
			 (bin = (Na_Ooff))
			 
		(reaction
		    (Na_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 Na_z )  
		   
		   )
			
		(component (type pore)
			   (const  gbar  = 0.015)
			   (output gbar ))
		
		(component (type permeating-substance) (name na)
			   (const e = ena)
			   (output e ))
		
		) ;; end Narsg component





   (component (type ion-channel) (name Ih)
	      
	      (component (type gate)
			 
			 ;; rate functions
			 
			 (Ih_inf = (1.0 /(1.0 + exp ((v + 90.1) / 9.9))))
			 
			 (Ih_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 (Ih_inf))
			   (m-power   1)
			   (h-power   0)
			   (m-inf     (Ih_inf))
			   (m-tau     (Ih_tau))
			   ))
			 
			 )

	      (component (type pore)
			 (const  gbar_Ih  = 0.0001)
			 (output gbar_Ih ))
	      
	      (component (type permeating-substance) (name non-specific)
			 (const e_Ih = -30)
			 (output e_Ih ))
	      
	      ) ;; end Ih current

	
   (component (type ion-channel) (name Leak)
	      
	      (component (type pore)
			 (const  gbar_Leak  = 5e-5)
			 (output gbar_Leak ))
	      
	      (component (type permeating-substance) (name non-specific)
			 (const e_Leak = -60)
			 (output e_Leak ))
	      
	      ) ;; 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 l = 20)
           (const r = 10)
           (const C_m = (2 * PI * l * r * 1e-6))
	   (output C_m))

))