(nemo-model CA3

  (
   (input v celsius)

   (const  F = 96485.3)
   (const  R = 8.3134)
   (const  KTOMV = .0853 (mV/degC)


   (component (type defaults) 
      (const V_t = -35)
      (const celsius = 23)
      (output celsius V_t)
      )


   (component (type membrane-capacitance)
  	 (const c = 1)
	 (output c  ))


   (component (type decaying-pool) (name ca)

      (input (ica from ion-currents))

   	      (const  d       = .1)
   	      (const  cao     = 2.)
   	      (const  cai0    = 50e-6)
   	      (const  tau     = 100.)

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

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

      (input
       (cai from ion-pools))
	      
      (component (type gate)


                 (const d1 = .84)
                 (const d2 = 1.)
                 (const k1 = .48e-3)
                 (const k2 = .13e-6)
                 (const abar = .28)
                 (const bbar = .48)

		 ;; rate functions
                 (defun exp1 (k d v)
                   (k * exp ((-2. * d * F * v) / (R * (273.15 + celsius))))
                   )

		 (defun alpha_c (v c)
                   ((c * abar) / (c + exp1 (k1 d1 v))))
		 
		 (defun beta_c (v c)
                   (bbar / (1 + c / exp1 (k2 d2 v))))

                 (mtau = (1 / (alpha_c (v c) + beta_c (v c))))
                 (minf = ((alpha_c (v c) * tau)))
             
		 (hh-ionic-gate 
		  (CaGK  ;; ion name: exported variables will be of the form {ion}_{id}
		   (initial-m  (minf))
		   (m-power    1)
		   (h-power    0)
                   (m-inf      (minf))
                   (m-tau      (mtau))
		   ))
		 
		 )
      
      (component (type pore)
		 (const  gbar  = 0.01)
		 (output gbar ))
      
      
      (component (type permeating-ion) (name k)
		 (e = ek)
		 (output e ))
      
      (component (type modulating-ion) (name ca) )
      
      ) ;; end KCa current


   
   (component (type voltage-clamp) (name CaGK)
	      
	      (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 CaL )

      (input
       (cai from ion-pools))
	      
      (component (type gate)

                 (const ki     = .001)
                 (const cai    = 50.e-6)
                 (const cao    = 2)
                 (const q10    = 5)
                 (const mmin   = 0.2)
                 (const tfa    = 1)
                 (const a0m    = 0.1)
                 (const zetam  = 2)
                 (const vhalfm = 4)
                 (const gmm    = 0.1)

                 (defun h2 (cai) (ki / (ki + cai)))
                 
                 (defun ghk  (v ci co celsius)
                   (let ((f (KTF (celsius) / 2))
                         (nu (v / f)))
                     (neg (f) * (1. - (ci / co) * exp (nu)) * efun (nu))
                     ))

                 (defun KTF (celsius)
                   ((25. / 293.15) * (celsius + 273.15))
                   )
                 
                 (defun efun (z)
                   (if (abs (z) < 1e-4) 
                       then (1 - z/2)
                       else (z / (exp (z) - 1)))
                   )

                 (defun alpha (v)
                   (15.69 * (neg (v) + 81.5) / (exp ((neg (v) + 81.5) / 10.0) - 1.0)))

                 (defun beta (v)
                   (0.29 * exp (neg (v) / 10.86)))

                 (defun alphamt (v)
                   (exp (0.0378 * zetam * (v - vhalfm)))
                   )

                 (defun betamt (v) 
                   (exp (0.0378 * zetam * gmm * (v - vhalfm))))

                 (qt = (pow (q10 ((celsius - 25) / 10))))
                 
                 (a = alpha (v))
                 (b = (1 / ((a + beta (v)))))
                 (minf = (a * b))

                 (tau0  = (betamt (v) / (qt * a0m * (1 + alphamt (v)))))
                 (mtau  = (if (tau0 < mmin / qt) then (mmin / qt) else tau0))

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

                 )


      (component (type permeability)

                 (p   = (gcalbar * ghk (v cai cao celsius) * h2 (cai)))

                 (output p)
                 )

	      
      (component (type permeating-ion) (name ca) 
                 )

      )

      
      ;; Borg-Graham type generic K-AHP channel

NEURON {
	SUFFIX KahpM95
	USEION k READ ek WRITE ik
        USEION ca READ cai
        RANGE  gbar,ik, gkahp
        GLOBAL inf,tau
}

UNITS {
	(mA) = (milliamp)
	(mV) = (millivolt)
}

PARAMETER {
	celsius = 6.3	(degC)
	gbar	= .003 	(mho/cm2)
        n	= 4
        cai	= 50.e-6 (mM)
        a0	= 1e8 (/ms-mM-mM-mM-mM)		:b0/(20e-4^4)
        b0	= .5e-2  (/ms)			:0.5/(0.100e3)
        v       	 (mV)
        ek      	 (mV)
	q10=3
}

STATE {	w }

ASSIGNED {
	ik 		(mA/cm2)
        gkahp  		(mho/cm2)
        inf
        tau
}

BREAKPOINT {
	SOLVE state METHOD cnexp
	gkahp = gbar*w
	ik = gkahp*(v-ek)
}

INITIAL {
	rate(cai)
	w=inf
}

FUNCTION alp(cai (mM)) {
  alp = a0*cai^n
}

DERIVATIVE state {     : exact when v held constant; integrates over dt step
        rate(cai)
        w' = (inf - w)/tau
}

PROCEDURE rate(cai (mM)) { :callable from hoc
        LOCAL a,qt
        qt=q10^((celsius-24)/10)
        a = alp(cai)
        tau = 1/(qt*(a + b0))
        inf = a*tau*qt
}


)