2018020/KATP_stochastic.mod at main · ModelDBRepository/2018020 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
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
		}
	}
}