TITLE multi-state reduced descrete Markov K-ATP channel

: ADP, ATP transitions require 4 partially indpendent sites - reduced to single transition
: fast 'flicker' state increases noise
: C0-O-Cf
: rates 0.079/0.025  :4.5 / 0.61


NEURON {
	SUFFIX katpStoch
	USEION adp READ adpi VALENCE 0
	USEION k WRITE ik
	THREADSAFE :but we make it so (this SHOULD not be a problem so long as each neuron or compartment is handled by one thread)
				: suppress errors associated with pointer
	RANGE   gkatpbar, nSites, unitCond, km, nOpen, nFast,gk, ik,alphm,n,ik

}

INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}

UNITS {
	(molar) = (1/liter)
	(mM) =  (millimolar)
	(um) =  (micron)
	(mA) =  (milliamp)
	FARADAY = (faraday) (coulomb) 
	PI = (pi) (1)
}

STATE {
   l : dummy variable method adapted from IH, NaV channels from Kole et al 2006  
}

PARAMETER {
		 v (mV)
         :dt (ms)
         adpi (mM)
         km = 0.0077 (mM)
         n = 1.3 (1)
         gkatpbar = 20e-6 (S/cm2)
         alphm = 0.079 (1/ms)
         betam = 0.025 (1/ms)
         cfast = 0: 4.5 (1/ms) 0 turns off flicker state
         ofast = 0.61 (1/ms) : might be flipped
         ek  =  -90.0 (mV)
         unitCond = 25e-12 (S) :25 pS per Tanner et al 2011
         seed
         }

ASSIGNED {
  ik  (mA/cm2)
  gk  
  dt (ms)
  nOpen
  nFast
  area (um2)
  nSites
}

INITIAL {
  nSites = ceil(fabs((1e-8*gkatpbar*area)/unitCond)) : yes the number of sites is a float, but the loop in noise uses a while loop.  This will lead to the conductance being slightly larger if not 
  l=0
  nOpen = nSites*sinf(adpi)*ofast/(ofast+cfast)
  set_seed(seed)
}

FUNCTION sinf(adpi){
   if(adpi > 1e-6){
     sinf = 1.0/(1.0 +  pow(km/adpi,n))
   } else {
     sinf = 0
   }
}



BREAKPOINT {
   SOLVE states METHOD cnexp
   if (nSites >0){
	gk = nOpen/nSites*gkatpbar
	} else {
	gk = 0
	}
   :ik = nOpen*(v-ek)*unitCond
   ik = gk*(v-ek)
}

DERIVATIVE states {
	l'=l
	noise(adpi)
}

PROCEDURE noise(adp) {
	LOCAL nClose,nOpenLoop,nFastLoop, nCloseLoop, pFO, pOF, pOC, pCO, aa, bb, rand,ii: duplicates to avoid looping over updating number of states
    nClose = nSites - nOpen - nFast
    nCloseLoop = nClose
    nOpenLoop = nOpen
    nFastLoop = nFast
    
    aa = alphm*sinf(adp)
    bb = alphm/2.0 : open/close rates balance when adpi=km
	pFO=1.0-exp(-dt*ofast)
	pOF=1.0-exp(-dt*cfast)
	pOC=1.0-exp(-dt*bb)
	pCO=1.0-exp(-dt*aa)
	FROM ii=1 TO nOpenLoop {
		if (scop_random()<=pOC) { 	: note this will produce unexpected results if pOC+pOF >= 1 (the solution to which is to decrease dt until sum of rates out << 1
			:nClose = nClose+1
			nOpen = nOpen-1	
		} else if (scop_random()<= pOF)	{			:scop_random uniform between 0 and 1
			nFast = nFast+1
			nOpen = nOpen-1
		}
	}
	
	FROM ii=1 TO nFastLoop {
		if (scop_random()<= pFO)	{			:scop_random uniform between 0 and 1
			nFast = nFast-1
			nOpen = nOpen+1
		} 
	}
	
	FROM ii=1 TO nCloseLoop {
	    if (scop_random()<=pCO) { 	: note this will produce unexpected results if pOC+pOF >= 1 (the solution to which is to decrease dt until sum of rates out << 1
			:nClose = nClose-1 : this is updated by normalization step
			nOpen = nOpen+1	
		}
	}
}