Stochastic models of ion channel gating/JAW SCM gillespie vclamp.m
From Scholarpedia
% gillespie_vclamp.m % runs the gillespie simultation for a two-state channel under voltage % clamp % function [output] = gillespie_vclamp NV = [1 10 100 1000]; % number of channels ts = 10; % time to switch V, ms tend = 2*ts; % time to end simulation, ms figure hold for i = 1:length(NV) [t, g] = gillespie_function(NV(i)); output{i,1} = [t' g']; stairs(t,g/NV(i)+(length(NV)+1-i)*1.2) % plot normalized, offset conductance % axis([0 tend*1.05 -0.1 1.1]) % normalize plot axes end % plot deterministic solution tdet = [0:0.05:20]; vi = -80; % initial voltage, mV vf = -20; % final voltage, mV [a0 b0 m0 ta0] = mparams(vi); % initial values of alpha, beta, minf and taum [a1 b1 m1 ta1] = mparams(vf); % final values of alpha, beta, minf, and taum gdet = m0 + (tdet>ts).*(m1-m0).*(1-exp(-(tdet-ts)/ta1)); output{length(NV)+1,1} = [tdet' gdet']; plot(tdet,gdet,'k') tscale = [15 20]; gscale = [-0.1 -0.1]; plot(tscale,gscale,'k') gdivide = [-1 6]; plot([ts ts],gdivide,':') end % main function function [tv,N_open] = gillespie_function(N) gamma = 10; % open-channel conductance, pS %N = 100; % number of channels vi = -80; % initial voltage, mV vf = -20; % final voltage, mV ts = 10; % time to switch V, ms tend = 2*ts; % time to end simulation, ms [a0 b0 m0 ta0] = mparams(vi); % initial values of alpha, beta, minf and taum [a1 b1 m1 ta1] = mparams(vf); % final values of alpha, beta, minf, and taum % seed random number generator rand('twister',sum(1000*clock)); No = round(m0*N); % initial number of open channels Nc = N-No; % initial number of closed channels t = 0; % initial time i = 1; % counter for vector of transition times % simulate first portion tv(1) = 0; N_open(1)=No; while (t<tend) i = i + 1; if (t<ts) lm(i) = Nc*a0 + No*b0; beta = b0; else lm(i) = Nc*a1 + No*b1; beta = b1; end r = rand(2,1); dt(i) = 1/lm(i)*log(1/r(1)); if (No*beta > r(2)*lm(i)) No = No - 1; Nc = Nc + 1; else No = No + 1; Nc = Nc - 1; end N_open(i)=No; N_closed(i)=Nc; t = t+dt(i); tv(i) = t; end % while loop end % gillespie_function % mparams(v) % returns taum and minf as functions of v % HH (6 deg C) parameters function [a,b,minf,taum] = mparams(v) a = alm(v); b = bem(v); taum = 1/(alm(v)+bem(v)); minf = alm(v)*taum; end function alm = alm(v) x = -(v+40); y = 10; if abs(x/y) < 1.e-6 alm = 0.1*y*(1 - x/y/2); % this "trap" catches divide by zero errors else alm = 0.1*x/(exp(x/y) - 1); end % if/then function end % function call function bem = bem(v) bem = 4*exp(-(v+65)/18); end