;; Reference: Theta-Frequency Bursting and Resonance in Cerebellar Granule Cells:Experimental
;; Evidence and Modeling of a Slow K+-Dependent Mechanism
;; Egidio D'Angelo,Thierry Nieus,Arianna Maffei,Simona Armano,Paola Rossi,Vanni Taglietti,
;; Andrea Fontana and Giovanni Naldi

(nemo-model Golgi

  (
   (input v celsius
    (cai from ion-pools)
    (ica from ion-currents)
    (ca2i from ion-pools)
    (ca2o from ion-pools)
    (ica2 from ion-currents))



   (component (type decaying-pool) (name ca)
   	      (const  F       = 96485.0)
   	      (const  d       = .2)
   	      (const  cao     = 2.)
   	      (const  cai0    = 50e-6)
   	      (const  beta    = 1.3)

   	      (d (ca) =  ((neg (ica) / (2 * F * d)) * 1e4 -
   			  (beta * (ca - cai0)))
   	      	      	 (initial cai0))
	      
   	      (output ca)
   	      )


   (component (type decaying-pool) (name ca2)
              (const valence  = 2)
   	      (const  F       = 96485.0)
   	      (const  d       = .2)
   	      (const  ca2o     = 2.)
   	      (const  ca2i0    = 50e-6)
   	      (const  beta    = 1.3)
	      
   	      (d (ca2) =  ((neg (ica2) / (2 * F * d)) * 1e4 -
   			  (beta * (ca2 - ca2i0)))
   	      	      	 (initial ca2i0))
	      
   	      (output ca2 valence)
   	      )

   
   (component (type ionic-current) (name CaHVA )
	      
	      (component (type gate)

			 ;; rate constants

			 (Q10 = (pow (3 ((celsius - 20) / 10))))

			 (const Aalpha_s  = 0.04944)
			 (const Kalpha_s  =  15.87301587302)
			 (const V0alpha_s = -29.06)
	
			 (const Abeta_s  = 0.08298)
			 (const Kbeta_s  =  -25.641)
			 (const V0beta_s = -18.66)

			 (const Aalpha_u  = 0.0013)
			 (const Kalpha_u  =  -18.183)
			 (const V0alpha_u = -48)
			 
			 (const Abeta_u = 0.0013)
			 (const Kbeta_u = 83.33)
			 (const V0beta_u = -48)
			 
			 ;; rate functions

			 (defun alpha_s (v Q10)
			   (Q10 * Aalpha_s * exp((v - V0alpha_s) / Kalpha_s)))

			 (defun beta_s (v Q10)
			   (Q10 * Abeta_s * exp((v - V0beta_s) / Kbeta_s)))

			 (defun alpha_u (v Q10)
			   (Q10 * Aalpha_u * exp((v - V0alpha_u) / Kalpha_u)))

			 (defun beta_u (v Q10)
			   (Q10 * Abeta_u * exp((v - V0beta_u) / Kbeta_u)))


			 (hh-ionic-gate 
			  (CaHVA  ;; ion name: exported variables will be of the form {ion}_{id}
			   (initial-m  ((alpha_s (v Q10))/(alpha_s (v Q10) + beta_s (v Q10)) ))
			   (initial-h  ((alpha_u (v Q10))/(alpha_u (v Q10) + beta_u (v Q10)) ))
			   (m-power    2)
			   (h-power    1)
			   (m-inf      (alpha_s (v Q10) / (alpha_s (v Q10) + beta_s (v Q10)) ))
			   (m-tau      (1 / (alpha_s (v Q10) + beta_s (v Q10)) ))
			   (h-inf      (alpha_u (v Q10) / (alpha_u (v Q10) + beta_u (v Q10)) ))
			   (h-tau      (1 / (alpha_u (v Q10) + beta_u (v Q10)) ))
			   ))
			 
			 )
	      
	      (component (type pore)
			 (const  gbar  = 460e-6)
			 (output gbar ))
	      
	      (component (type permeating-ion) (name ca)
			 (const e = 129.33)
			 (output e ))
	      
	      ) ;; end CaHVA current



   (component (type voltage-clamp) (name CaHVA)
	      
	      (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 ionic-current) (name CaLVA )
	      
	      (component (type gate)

                         (const shift   = 2) ; screening charge for Ca_o = 2 mM

                         (const v0_m_inf = -50)
                         (const v0_h_inf = -78)
                         (const k_m_inf  = -7.4)
                         (const k_h_inf  = 5.0)
	
                         (const C_tau_m   = 3)
                         (const A_tau_m   = 1.0)
                         (const v0_tau_m1 = -25)
                         (const v0_tau_m2 = -100)
                         (const k_tau_m1  = 10)
                         (const k_tau_m2 = -15)
	
                         (const C_tau_h   = 85)
                         (const A_tau_h   = 1.0)
                         (const v0_tau_h1 = -46)
                         (const v0_tau_h2 = -405)
                         (const k_tau_h1  = 4)
                         (const k_tau_h2  = -50)

			 
			 ;; rate functions

                         (phi_m = (pow (5.0 ((celsius - 24) / 10))))
                         (phi_h = (pow (3.0 ((celsius - 24) / 10))))
	
                         (m_inf = (1.0 / ( 1 + exp ((v + shift - v0_m_inf) / k_m_inf)) ))
                         (h_inf = (1.0 / ( 1 + exp ((v + shift - v0_h_inf) / k_h_inf)) ))
	
                         (tau_m = ( (C_tau_m + A_tau_m / ( exp ((v + shift - v0_tau_m1) / k_tau_m1) + 
                                                          exp ((v + shift - v0_tau_m2) / k_tau_m2) ) ) / phi_m) )

                         (tau_h = ( (C_tau_h + A_tau_h / ( exp ((v + shift - v0_tau_h1 ) / k_tau_h1) + 
                                                          exp ((v + shift - v0_tau_h2) / k_tau_h2) ) ) / phi_h) )


			 (hh-ionic-gate 
			  (CaLVA  ;; ion name: exported variables will be of the form {ion}_{id}
			   (initial-m  m_inf)
			   (initial-h  h_inf)
			   (m-power    2)
			   (h-power    1)
			   (m-inf      m_inf)
			   (m-tau      tau_m)
			   (h-inf      h_inf)
			   (h-tau      tau_h)
			   ))
			 )
	      

	      (component (type pore)
			 (const  gbar  = 2.5e-4)
			 (output gbar ))

	      
	      (component (type permeating-ion) (name ca2)
                         (const F = 96485.0)
                         (const R = 8.3145)
			 (e = ((1e3) * (R * (celsius + 273.15)) / (2 * F) * log10 (ca2o / ca2i)))
			 (output e ))

	      
	      ) ;; end CaLVA current



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


      
   (component (type ionic-current) (name HCN1)

   	      (const Ehalf = -72.49)
   	      (const c = 0.11305)
	      
   	      (const rA = 0.002096)
   	      (const rB = 0.97596)
	      
   	      (defun r (potential)
   		((rA * potential) + rB))
	      
   	      (defun tau (potential t1 t2 t3)
   		(exp (((t1 * potential) - t2) * t3)))
	      
   	      (defun o_inf (potential Ehalf c)
   		(1 / (1 + exp ((potential - Ehalf) * c))))
              

   	      (component (type gate)
                         

   			 ;; rate constants
	
                         (const tCs = 0.01451)
                         (const tDs = -4.056)
                         (const tEs = 2.302585092)

                         (o_slow_inf = ((1 - r (v)) * o_inf (v Ehalf c)))
			
                         (tau_s =  (tau (v tCs tDs tEs)) )

                         (d (o_slow) = ((o_slow_inf - o_slow) / tau_s) 
                            (initial o_slow_inf))
                         
                         (output o_slow)

   			 )
	      
   	      (component (type gate)
                         

   			 ;; rate constants
	
                         (const tCf = 0.01371)
                         (const tDf = -3.368)
                         (const tEf = 2.302585092)

                         (o_fast_inf = (r (v) * o_inf (v Ehalf c)))
			
                         (tau_f =  (tau (v tCf tDf tEf)) )

                         (d (o_fast) = ((o_fast_inf - o_fast) / tau_f) 
                            (initial o_fast_inf))
                         
                         (output o_fast)

   			 )
	      
   	      (component (type pore)
   			 (const  gbar  = 5e-5)
   			 (output gbar))
	      
   	      (component (type permeating-ion) (name non-specific)
   			 (const e = -20)
   			 (output e ))
	      
   	      ) ;; end HCN1 current


              
   (component (type voltage-clamp) (name HCN1)
	      
   	      (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 ionic-current) (name HCN2)
                         

   	      (const Ehalf = -81.95)
   	      (const c = 0.1661)
	      
   	      ;; rate constants

   	      (const rA = -0.0227)
   	      (const rB = -1.4694)

   	      (defun r (potential r1 r2)
   		(if (potential >= -64.70) 
   		    then  0.0
   		    else (if (potential <= -108.70) 
   			     then 1.0
   			     else (r1 * potential) + r2)))

   	      (defun o_inf (potential Ehalf c)
   		(1 / (1 + exp((potential - Ehalf) * c))))
              
   	      (component (type gate)
                         

   			 ;; rate constants

                         (const tCs = 0.0152)
                         (const tDs = -5.2944)
                         (const tEs = 2.3026)
                         
                         (defun tau_slow (potential t1 t2 t3)
                           (exp (t3 * ((t1 * potential) - t2))))

                         (o_slow_inf = ((1 - r (v rA rB)) * o_inf (v Ehalf c)))

                         (tau_s =  (tau_slow(v tCs tDs tEs)))

                         (d (o_slow) = ((o_slow_inf - o_slow) / tau_s) 
                            (initial o_slow_inf))

                         (output o_slow)

   			 )
	      
   	      (component (type gate)
                         

   			 ;; rate constants

                         (const tCf = 0.0269)
                         (const tDf = -5.6111)
                         (const tEf = 2.3026)
                         
                         (defun tau_fast (potential t1 t2 t3)
                           (exp (t3 * ((t1 * potential) - t2))))

                         (o_fast_inf = (r (v rA rB) * o_inf (v Ehalf c)))

                         (tau_f =  (tau_fast(v tCf tDf tEf)))

                         (d (o_fast) = ((o_fast_inf - o_fast) / tau_f) 
                            (initial o_fast_inf))

                         (output o_fast)

   			 )
	      
   	      (component (type pore)
   			 (const  gbar  = 8e-5)
   			 (output gbar))
	      
   	      (component (type permeating-ion) (name non-specific)
   			 (const e = -20)
   			 (output e ))
	      
   	      ) ;; end HCN2 current



   (component (type voltage-clamp) (name HCN2)
	      
   	      (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 ionic-current) (name KA )

	      (defun sigm (x y) (1 / (exp (x / y) + 1)))
	      
	      (component (type gate)
			 
			 ;; rate constants
			 
			 (Q10 = (pow (3 ((celsius - 25.5) / 10))))

			 (const Aalpha_a  = 0.8147)
			 (const Kalpha_a  = -23.32708)
			 (const V0alpha_a = -9.17203)
			 (const Abeta_a   = 0.1655)
			 (const Kbeta_a   = 19.47175)
			 (const V0beta_a  = -18.27914)

			 (const Aalpha_b  = 0.0368)
			 (const Kalpha_b  = 12.8433)
			 (const V0alpha_b = -111.33209)
			 (const Abeta_b   = 0.0345)
			 (const Kbeta_b   = -8.90123)
			 (const V0beta_b  = -49.9537)

			 (const V0_ainf = -38)
			 (const  K_ainf = -17)

			 (const V0_binf   = -78.8)
			 (const K_binf    = 8.4)

			 
			 ;; rate functions

			 (defun alpha_a (v Q10)
			   (Q10 * Aalpha_a * sigm((v - V0alpha_a) Kalpha_a)))

			 (defun beta_a (v Q10)
			   (Q10 * (Abeta_a / exp((v - V0beta_a) / Kbeta_a))))

			 (defun alpha_b (v Q10)
			   (Q10 * Aalpha_b * sigm((v - V0alpha_b) Kalpha_b)))

			 (defun beta_b (v Q10)
			   (Q10 * Abeta_b * sigm((v - V0beta_b) Kbeta_b)))



			 (a_inf = (1 / (1 + exp ((v - V0_ainf) / K_ainf)))) 
			 (tau_a = (1 / (alpha_a (v Q10) + beta_a (v Q10)) ))
			 (b_inf = (1 / (1 + exp ((v - V0_binf) / K_binf))))
			 (tau_b = (1 / (alpha_b (v Q10) + beta_b (v Q10)) ))


			 (hh-ionic-gate 
			  (KA  ;; ion name: exported variables will be of the form {ion}_{id}
			   (initial-m  (a_inf))
			   (initial-h  (b_inf))
			   (m-power    3)
			   (h-power    1)
			   (m-inf      a_inf)
			   (m-tau      tau_a)
			   (h-inf      b_inf)
			   (h-tau      tau_b)
			   ))

			 )

	      	      
	      (component (type pore)
			 (const  gbar  = 0.008)
			 (output gbar ))

	      
	      (component (type permeating-ion) (name k)
			 (const e = -84.69)
			 (output e ))
	      
	      ) ;; end KA current
	      


   (component (type voltage-clamp) (name KA)
	      
	      (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 ionic-current) (name KCa )
	      
	      (component (type gate)

			 ;; rate constants
			 (Q10 = (pow (3 ((celsius - 30) / 10))))
			 
			 (const Aalpha_c = 7)
			 (const Balpha_c = 1.5e-3)

			 (const Kalpha_c =  -11.765)

			 (const Abeta_c = 1.)
			 (const Bbeta_c = 0.15e-3)

			 (const Kbeta_c = -11.765)

			 ;; rate functions
			 (defun alpha_c (v cai Q10)
			   (Q10 * Aalpha_c / (1 + (Balpha_c * exp(v / Kalpha_c) / cai)))) 

			 (defun beta_c (v cai Q10)
			   (Q10 * Abeta_c / (1 + (cai / (Bbeta_c * exp (v / Kbeta_c))) )))


			 (hh-ionic-gate 
			  (KCa  ;; ion name: exported variables will be of the form {ion}_{id}
			   (initial-m  ((alpha_c (v cai Q10)) / (alpha_c (v cai Q10) + beta_c (v cai Q10)) ))
			   (m-power    1)
			   (h-power    0)
			   (m-alpha    (alpha_c (v cai Q10)))
			   (m-beta     (beta_c (v cai Q10)))
			   ))
			 
			 )

	      (component (type pore)
			 (const  gbar  = 0.003)
			 (output gbar ))

	      
	      (component (type permeating-ion) (name k)
			 (const e = -84.69)
			 (output e ))

	      (component (type modulating-ion) (name ca) )
	      
	      ) ;; end KCa current

	      (component (type voltage-clamp) (name KCa)
           
			 (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 ionic-current) (name KM )
	      
	      (component (type gate)

			 ;; rate constants
			 (const Aalpha_n = 0.0033)

			 (const Kalpha_n  = 40)
			 (const V0alpha_n = -30)
			 (const Abeta_n   = 0.0033)

			 (const Kbeta_n  = -20)
			 (const V0beta_n = -30)
			 (const V0_ninf  = -35)
			 (const   B_ninf = 6)
			 
			 (Q10 = (pow (3 ((celsius - 22) / 10))))
			 
			 ;; rate equations
			 (defun alpha_n (v Q10)
			   (Q10 * Aalpha_n * exp((v - V0alpha_n) / Kalpha_n) ))

			 (defun beta_n (v Q10)
			   (Q10 * Abeta_n * exp((v - V0beta_n) / Kbeta_n) ))

			 (hh-ionic-gate 
			  (KM  ;; ion name: exported variables will be of the form {ion}_{id}
			   (initial-m  (1 / (1 + exp((neg (v - V0_ninf)) / B_ninf))))
			   (m-power    1)
			   (h-power    0)
			   (m-tau      (1 / (alpha_n(v Q10) + beta_n (v Q10)) ))
			   (m-inf      (1 / (1 + exp((neg (v - V0_ninf)) / B_ninf))))
			   ))
			 )
	      
	      (component (type pore)
			 (const  gbar  = 0.001)
			 (output gbar ))
	      
	      (component (type permeating-ion) (name k)
			 (const e = -84.69)
			 (output e ))
	      
	      ) ;; end KM current



   (component (type voltage-clamp) (name KM)
	      
	      (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 ionic-current) (name KV )
	      
	      (defun linoid (x y) 
		(if ((abs (x / y)) < 1e-6) 
		    then (y * (1 - ((x / y) / 2)))
		    else (x / (exp (x / y) - 1))
		    ))


	      (component (type gate)

			 ;; rate constants

			 (const Aalpha_n = -0.01)
			 (const Kalpha_n = -10)
			 (const V0alpha_n = -26)
			 (const Abeta_n = 0.125)
	
			 (const Kbeta_n = -80)
			 (const V0beta_n = -36)

			 ;; rate functions
			 (Q10 = (pow (3 ((celsius - 6.3) / 10))))

			 (defun alpha_n (v Q10) 
			   (Q10 * Aalpha_n * linoid ((v - V0alpha_n) Kalpha_n)))

			 (defun beta_n (v Q10) 
			   (Q10 * Abeta_n * exp((v - V0beta_n) / Kbeta_n) ))

			 (hh-ionic-gate 
			  (KV  ;; ion name: exported variables will be of the form {ion}_{id}
			   (initial-m  ((alpha_n (v Q10)) / (alpha_n (v Q10) + beta_n (v Q10)) ))
			   (m-power    4)
			   (h-power    0)
			   (m-alpha    (alpha_n(v Q10)))
			   (m-beta     (beta_n(v Q10)))
			   ))
			 )

	      (component (type pore)
			 (const  gbar  = 0.032)
			 (output gbar ))
	      
	      (component (type permeating-ion) (name k)
			 (const e = -84.69)
			 (output e ))
	      
	      ) ;; end KV current
	      

   (component (type voltage-clamp) (name KV)
	      
	      (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 ionic-current) (name SK2 )
	      
   	      (component (type gate)

                         (Q10 = (pow (3 ((celsius - 23) / 10))))

                         (const diff = 3) ;; Diffusion factor

                         ;; rates ca-independent
                         (const invc1 = 80e-3)
                         (const invc2 = 80e-3)
                         (const invc3 = 200e-3)

                         (const invo1 = 1)
                         (const invo2 = 100e-3)
                         (const diro1 = 160e-3)
                         (const diro2 = 1.2)

                         ;; rates ca-dependent
                         (const dirc2 = 200)
                         (const dirc3 = 160)
                         (const dirc4 = 80)

                         (invc1_t = (invc1 * Q10  ))
                         (invc2_t = (invc2 * Q10 ))
                         (invc3_t = (invc3 * Q10 ))
                         (invo1_t = (invo1 * Q10 ))
                         (invo2_t = (invo2 * Q10 ))
                         (diro1_t = (diro1 * Q10 ))
                         (diro2_t = (diro2 * Q10 ))
                         (dirc2_t = (dirc2 * Q10 ))
                         (dirc3_t = (dirc3 * Q10 ))
                         (dirc4_t = (dirc4 * Q10 ))


                         (dirc2_t_ca = (dirc2_t * (cai / diff) ))
                         (dirc3_t_ca = (dirc3_t * (cai / diff) ))
                         (dirc4_t_ca = (dirc4_t * (cai / diff) ))


                         (reaction
   			  (SK2_z
   			   (transitions 
                            (<-> c1 c2 dirc2_t_ca invc1_t)
                            (<-> c2 c3 dirc3_t_ca invc2_t) 
                            (<-> c3 c4 dirc4_t_ca invc3_t) 
                            (<-> c3 o1 diro1_t invo1_t) 
                            (<-> c4 o2 diro2_t invo2_t) 
                            )
                           (conserve  (1 = (c1 + c2 + c3 + c4 + o2 + o1)))
   			   (open o1 o2) 
                           (power 1)))

                         (output SK2_z)
 			 
   			 )

   	      (component (type pore)
   			 (const  gbar  = 0.038)
   			 (output gbar ))

	      
   	      (component (type permeating-ion) (name k)
   			 (const e = -84.69)
   			 (output e ))

   	      (component (type modulating-ion) (name ca) )
	      
   	      ) ;; end SK2 current



   (component (type voltage-clamp) (name SK2)
           
   	      (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 ionic-current) (name Na )


   (defun linoid (x y) 
    	 (if ((abs (x / y)) < 1e-6) 
             then (y * (1 - ((x / y) / 2)))
	     else (x / (1 - exp (x / y) ))
	      ))

	      
   	      (component (type gate)

   			 ;; rate constants
   			 (Q10 = (pow (3 ((celsius - 20) / 10))))

   			 (const Aalpha_u  = 0.3)
   			 (const Kalpha_u  = -10)
   			 (const V0alpha_u = -25)
	
   			 (const Abeta_u  = 12)
   			 (const Kbeta_u  = -18.182)
   			 (const V0beta_u = -50)

   			 (const Aalpha_v  = 0.21)
   			 (const Kalpha_v  = -3.333)
   			 (const V0alpha_v = -50)
 
   			 (const Abeta_v  = 3)
   			 (const Kbeta_v  = -5)
   			 (const V0beta_v = -17)

   			 ;; rate functions
   			 (defun alpha_u (v Q10)	
   			   (Q10 * Aalpha_u * linoid((v - V0alpha_u) Kalpha_u) ))

   			 (defun beta_u (v Q10)
   			   (Q10 * Abeta_u * exp((v - V0beta_u) / Kbeta_u) ))

   			 (defun alpha_v (v Q10)
   			   (Q10 * Aalpha_v * exp((v - V0alpha_v) / Kalpha_v) ))	
				
   			 (defun beta_v (v Q10)
   			   (Q10 * Abeta_v / (1 + exp((v - V0beta_v) / Kbeta_v) )))

			 
   			 (hh-ionic-gate 
   			  (Na  ;; ion name: exported variables will be of the form {ion}_{id}
   			   (initial-m  ((alpha_u (v Q10)) / (alpha_u (v Q10) + beta_u (v Q10)) ))
   			   (initial-h  ((alpha_v (v Q10)) / (alpha_v (v Q10) + beta_v (v Q10)) ))
   			   (m-power    3)
   			   (h-power    1)
   			   (m-alpha    (alpha_u (v Q10)))
   			   (m-beta     (beta_u (v Q10)))
   			   (h-alpha    (alpha_v (v Q10)))
   			   (h-beta     (beta_v (v Q10)))
   			   ))
			 
   			 )
	      
   	      (component (type pore)
   			 (const  gbar  = 0.048)
   			 (output gbar ))

	      
   	      (component (type permeating-ion) (name na)
   			 (const e = 87.39)
   			 (output e ))
	      
   	      ) ;; end Na current


   (component (type voltage-clamp) (name Na)

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

   
   (component (type ionic-current) (name NaP )
	      
   	      (defun linoid (x y) 
   		(if ((abs (x / y)) < 1e-6) 
   		    then (y * (1 - ((x / y) / 2)))
   		    else (x / (exp (x / y) - 1))
   		    ))


   	      (component (type gate)

   			 ;; rate constants
   			 ( Q10 = (pow (3 ((celsius - 30) / 10))))

                         (const Aalpha_m  = -0.91)
                         (const Kalpha_m  = -5)
                         (const V0alpha_m = -40)
                         (const Abeta_m	  = 0.62)
                         (const Kbeta_m   = 5)
                         (const V0beta_m  = -40)
                         (const V0_minf   = -43)
                         (const B_minf    = 5)

   			 ;; rate functions

   			 (defun alpha_m (v Q10) 
   			   (Q10 * Aalpha_m * linoid( ( v - V0alpha_m ) Kalpha_m) ))

   			 (defun beta_m (v Q10) 
   			   (Q10 * Abeta_m * linoid( ( v - V0beta_m )  Kbeta_m ) ))

   			 (hh-ionic-gate 
   			  (NaP  ;; ion name: exported variables will be of the form {ion}_{id}
   			   (initial-m  (1 / (1 + exp ((neg (v - V0_minf)) / B_minf))))
   			   (m-power    1)
   			   (h-power    0)
   			   (m-inf     (1 / (1 + exp ((neg (v - V0_minf)) / B_minf))))
   			   (m-tau     (5 / (alpha_m (v Q10) + beta_m (v Q10)) ))
   			   ))
   			 )

   	      (component (type pore)
   			 (const  gbar  = 0.00019)
   			 (output gbar ))

   	      (component (type permeating-ion) (name na)
   			 (const e = 87.39)
   			 (output e ))
	      
   	      ) ;; end NaP current


   (component (type voltage-clamp) (name NaP)

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

   
   (component (type ionic-current) (name NaR )
	      
   	      (component (type gate)

   			 ;; rate constants
   			 (Q10 = (pow (3 ((celsius - 20) / 10))))

   			 (const Aalpha_s     = -0.00493)
   			 (const V0alpha_s    = -4.48754)
   			 (const Kalpha_s     = -6.81881)
   			 (const Shiftalpha_s = 0.00008)

   			 (const Abeta_s     = 0.01558)
   			 (const V0beta_s    = 43.97494)
   			 (const Kbeta_s     =  0.10818)
   			 (const Shiftbeta_s = 0.04752)

   			 (const Aalpha_f  = 0.31836)
   			 (const V0alpha_f = -80)
   			 (const Kalpha_f  = -62.52621)

   			 (const Abeta_f  = 0.01014)
   			 (const V0beta_f = -83.3332)
   			 (const Kbeta_f  = 16.05379)

   			 ;; rate functions
   			 (defun alpha_s (v Q10) 
   			   (Q10 * (Shiftalpha_s + ((Aalpha_s * ((v + V0alpha_s) )) / (exp ((v + V0alpha_s) / Kalpha_s) - 1)))))

   			 (defun beta_s (v Q10) 
                           (let ((x1 ((v + V0beta_s) / Kbeta_s)))
   			     (Q10 * (Shiftbeta_s + Abeta_s * ((v + V0beta_s) / (exp((if (x1 > 200) then 200 else x1)) - 1))))))

   			 (defun alpha_f (v Q10) 
   			   (Q10 * Aalpha_f * exp( ( v - V0alpha_f ) / Kalpha_f)))

   			 (defun beta_f (v Q10) 
   			   (Q10 * Abeta_f * exp( ( v - V0beta_f ) / Kbeta_f )  ))

   			 (hh-ionic-gate 
   			  (NaR  ;; ion name: exported variables will be of the form {ion}_{id}
   			   (initial-m  ((alpha_s (v Q10)) / (alpha_s (v Q10) + beta_s (v Q10)) ))
   			   (initial-h  ((alpha_f (v Q10)) / (alpha_f (v Q10) + beta_f (v Q10)) ))
   			   (m-power    1)
   			   (h-power    1)
   			   (m-alpha    (alpha_s (v Q10)))
   			   (m-beta     (beta_s (v Q10)))
   			   (h-alpha    (alpha_f (v Q10)))
   			   (h-beta     (beta_f (v Q10)))
   			   ))
   			 )

   	      (component (type pore)
   			 (const  gbar  = 0.0017)
   			 (output gbar ))

   	      (component (type permeating-ion) (name na)
   			 (const e = 87.39)
   			 (output e ))
	      
   	      ) ;; end Nar current


   (component (type voltage-clamp) (name NaR)

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

      
   (component (type ionic-current) (name Lkg)
	      
	      (component (type pore)
			 (const  gbar  = (21e-6))
			 (output gbar))
	      
	      (component (type permeating-ion) (name non-specific)
			 (const e = -55)
			 (output e ))
	      
	      ) ;; end leak current

   ))