Page 1 of 1

Stochastic gating as a function of external variable (light, membrane tension)

Posted: Tue May 14, 2019 6:30 am
by aNNe
I would like to implement a mechanism that reflects ion channels gated by an external parameter, such as light intensity or membrane tension. The rate's dependency on the external parameter is not simply linear (quadratic, cubic), as in the case of gating triggered by binding one (two, three) ligand molecules. Instead, the external parameter influences the energy of the open and closed states similar to the impact of voltage on voltage-gated ion channels.
In the mean-field, deterministic case, there is no problem. I can write a .mod file where rates depend on the external parameter:
~ C <-> O ( maxrate/(exp(-(extPara- s_half_o)/k_o) + 1), maxrate/(exp((extPara- s_half_c)/k_c) + 1) )
and then provide a time varying extPara via from the .hoc file:,tvec,1)

But there are only about 20 externally gated channels in a given cell. So actual channel counts and stochasticity really matter. I am also interested in parameters, such as, "time to first opening" and cannot simply add a noise term to an average current.

I was hoping to use the Channel Builder to create a stochastic (point process) version of the mechanism, but as far as I can see, the GUI restricts the impact of external parameters to the case of ligands.

Does the Channel Builder provide a way to implement a channel using gating transitions with exponential dependence on external variables?

If the functional forms of the transition rates are fixed, could I create a table of rates and pick the relevant entry based on the value of an external parameter?

If the Channel Builder cannot provide a solution, is there a tried and true way to implement stochastic gating with an NMODL mechanism? I am aware of some implementations, that have been recommended earlier in the forum, e.g.the stochastic Ih for Kole et al. 2006. However, I have the feeling that the published implementation is close to the intended statistics only if the number of transitions per time step is much much smaller the total number of channels. That is likely with very small time steps and very many channels. In my particular case, the channel count is small and I would like to implement a mechanism that generalizes well for arbitrary channel numbers and also arbitrary gating schemes, not just C<->O. Therefore I was looking for a general implementation of stochastic gating but could not find one.

Any ideas on either Channel Builder with external variables or NMODL mechanisms for stochastic gating of a general gating scheme with dependence on external variables?

Re: Stochastic gating as a function of external variable (light, membrane tension)

Posted: Wed May 22, 2019 5:09 am
by aNNe
Because the timing of individual openings and the validity for small copy numbers mattered in my application, I ended up implementing a Gillespie algorithm in a NMODL mechanism. It reproduces the kinetics of the mean field variant in the average of many reproductions and also for many channels. Surprisingly, it does not limit the simulation time. Probably, because even with 100 channels, no transition occurs in most of the time steps. Therefore, after the first random number is drawn, the update is over and in which case the mechanism is faster than any implementation of a differential equation solver or any mechanism that always updates all channels separately.

Because it runs orders of magnitude faster than some of the brute-force implementations around and because I am not aware of a Gillespie algorithm implemented in NMODL, I thought it might could serve as a useful reference, even though the code is specific for two-state models and needs some extension for a general case.

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

	SUFFIX 	stomet
	GLOBAL 	displacement
	RANGE	gbar, Ntotal, N_open, updateCount

   (mV) = (millivolt)
   (cm2)= (centimeter2)
   (mA) = (milliamp)
   (S) 	= (siemens)

	dt   (ms)
	v    (mV)
	rate (1/ms)
	time (ms)
	closingRate (1/ms) 
	openingRate (1/ms)
	rateCO	    (1/ms)
	rateTOTAL   (1/ms)
	: 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)

	displacement (nm)
	icat (mA/cm2)
	gbar (S/cm2)

	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

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

Re: Stochastic gating as a function of external variable (light, membrane tension)

Posted: Wed May 22, 2019 10:02 am
by hines
I'm happy your nmodl mechanism resolves the issue. So this response is a bit stale. The channel builder with ligand would also suffice, Just
invent a ligand name with 0 charge and have another mechanism that modifies the ligand concentration according to your desire function.

Re: Stochastic gating as a function of external variable (light, membrane tension)

Posted: Wed May 22, 2019 8:41 pm
by aNNe
and have another mechanism that modifies the ligand concentration according to your desire function
Oh, I never thought of that! Thanks for this idea!

Also, could you please tell me on which basis, the KSChan method creates gating events. Is it also using the Gillespie algorithm, mapped onto discrete time steps?