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 vector.play from the .hoc file:
svec.play(&extPara_mec,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?
Thanks!
Andreas
Stochastic gating as a function of external variable (light, membrane tension)
Re: Stochastic gating as a function of external variable (light, membrane tension)
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.
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 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
Re: Stochastic gating as a function of external variable (light, membrane tension)
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.
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)
Oh, I never thought of that! Thanks for this idea!and have another mechanism that modifies the ligand concentration according to your desire function
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?
Thanks!
Andreas