www.neuron.yale.edu

Code: Select all

: an ion channel with gating that depends on an external variable other than voltage
: this could be membrane stretch or light or something else. Here it is called "displacement"
: the external variable has to be played in via a vector.play in the hoc
: stochastic gating is implemented via Gillespie algorithm, 
: performance should be checked when Ntotal - number of channels - increases
: for this two state implementation, 10 to 1000 channels ran as fast as the mean field implementation
: 05/2019 andreas neef
 

NEURON {
	SUFFIX 	stomet
	GLOBAL 	displacement
	NONSPECIFIC_CURRENT icat
	RANGE	gbar, Ntotal, N_open, updateCount
}

UNITS	{
   (mV) = (millivolt)
   (cm2)= (centimeter2)
   (mA) = (milliamp)
   (S) 	= (siemens)
}

PARAMETER {
	dt   (ms)
	v    (mV)
	rate (1/ms)
	time (ms)
	updateCount
	Ntotal
	N_open
	closingRate (1/ms) 
	openingRate (1/ms)
	rateCO	    (1/ms)
	rateTOTAL   (1/ms)
}	
CONSTANT {	
	: parameters estimated from figures in Beurg et al. 2014 for outer hair cell 
	: mechanotransduction channels
	maxrate	= 1.69 (1/ms)
	s_half_o= 204  (nm)
	s_half_c= 15   (nm)
	k_o		= 50.6 (1/nm)
	k_c		= 51.4 (1/nm)
}		

ASSIGNED {
	displacement (nm)
	icat (mA/cm2)
	gbar (S/cm2)
	}
STATE {C O}

BREAKPOINT {
	closingRate = maxrate/(exp( (displacement - s_half_c)/k_c) + 1) 
	openingRate = maxrate/(exp(-(displacement - s_half_o)/k_o) + 1)
	rateCO    = C*openingRate
	rateTOTAL = rateCO +  O*closingRate
	SOLVE kinetics METHOD sparse
	icat = v * O * gbar / Ntotal
}

INITIAL	{
	C = Ntotal - N_open
	O = N_open
	time = 0
	closingRate = maxrate/(exp( (displacement - s_half_c)/k_c) + 1) 
	openingRate = maxrate/(exp(-(displacement - s_half_o)/k_o) + 1)
	rateCO=         C*openingRate
	: the rate with which a transition from C to O occurs:
	: occupancy of C times opening rate
	rateTOTAL = rateCO +  O*closingRate
	: the rate at which ANY transition occurs
}

KINETIC kinetics {
	~ C <-> O ( 0 , 0) 	: just here in case the call "SOLVE"
						: requires a recipe like this
						
	time = updateTime()
	: presumably, most updates end here, because the update time is beyond dt/2
	: only for Markov models this is not a problem, because the process has no memory
	: one dt/2 later the transition probabilities are unchanged 
	WHILE ( time < dt/2) : using dt/2 instead of dt, because the block is executed twice per dt update
	{
	updateStates()
	updateRates()
	time = time + updateTime()
	}
	updateCount = updateCount + 1
	: checks how often the entire update 
	: has been run
	: just there to debug
	: it turns out, the kinetics is call TWICE per dt (as described in the NMODL instructions)
	: therefore, the limit in the WHILE clause is now set to dt/2
	: in principle the updateCount could now be removed, but it is reassuring to be able to check whether 
	: a simulation of T/dt time steps produced indeed 2*T/dt update steps 
	N_open = O
}
PROCEDURE updateStates(){
: this determines (stochastically)
: which transition (out of all possible ones) actually occured
: and updates the states
	LOCAL ran
	ran=scop_random()
	ran=ran*rateTOTAL
	IF (ran < rateCO) : an opening event has occured
	: in the general case, here not only an upper bound has to be tested, but multiple comparisons with several brakets 
	: have to be made. Ideally, sort the brakets first, so that the most likely transition is tested first
	: the sorting overhead might also kill performance, so test whether its worth it
	{
	C = C - 1
	O = O + 1 
	}
	: this is the place to add more states
	ELSE : something has occured, so realize the last option
	{
	C = C + 1
	O = O - 1 
	}	
}
	
PROCEDURE updateRates() {
	rateCO    = C*openingRate
	rateTOTAL = rateCO +  O*closingRate
}

FUNCTION updateTime(){
: this draws the time intervall until the next event
: and returns it
	LOCAL ran
	ran=scop_random()
	updateTime = -log(ran)/rateTOTAL
}
	: it should be possible to use the random number from exponential distribution instead
	: but I have zero documentation
	: best guess is this line:
	: updateTime = exprand()/rateTOTAL