// Potassium channel from original HH model // Voltage clamp simulations with non-stationary noise analysis // UNcoupled activation particles (2-state independent particles), // Goldwyn et al. (Phys Rev E 83:041908 (2011)) implementation of the // Diffusion Approximation. Coupled activation particles WITHOUT // steady state approximation of variables in the stochastic terms // See "StochHH_K2 F1 Vclamp noise.sci" for more comments stacksize('max'); nsim=200; Tstop=6; dt=0.001; //ms points = round(Tstop/dt) NK=300; //Iamp=10; Vhold=-90; Vtest=70; rand('normal'); tic(); p=1; Norec = zeros(points,nsim); // mrec = zeros(points,nsim); // nrec=mrec; v = Vhold*ones(1,nsim); an=0.01*(v+55)./(1-exp(-(v+55)/10)); bn=0.125*exp(-(v+65)/80); N=an./bn; Kstatesum=(1+N)^4; n=[ones(1,nsim);4*N;6*N.^2;4*N.^3;N.^4]./(ones(5,1)*Kstatesum); v = Vtest; an=0.01*(v+55)./(1-exp(-(v+55)/10)); bn=0.125*exp(-(v+65)/80); tint=1; for tt=dt:tint:Tstop for t = tt:dt:tt+tint-dt Norec(p,:) = n(5,:)*NK; p=p+1; trans_n=[-(bn+3*an).*n(2,:)+4*an.*n(1,:)+2*bn.*n(3,:); -(2*bn+2*an).*n(3,:)+3*an.*n(2,:)+3*bn.*n(4,:); -(3*bn+an).*n(4,:)+2*an.*n(3,:)+4*bn.*n(5,:); -4*bn.*n(5,:)+an.*n(4,:)]; trans_n=[-(bn+3*an).*n(2,:)+4*an.*n(1,:)+2*bn.*n(3,:); -(2*bn+2*an).*n(3,:)+3*an.*n(2,:)+3*bn.*n(4,:); -(3*bn+an).*n(4,:)+2*an.*n(3,:)+4*bn.*n(5,:); -4*bn.*n(5,:)+an.*n(4,:)]; na=abs(n); for ii = 1:nsim //because this algorithm needs a matrix operation (line 90) it has to be done iterating over independent simulations Dmtx = [4*an*na(1,ii)+(3*an+bn)*na(2,ii)+2*bn*na(3,ii),-(3*an*na(2,ii)+2*bn*na(3,ii)),0,0; -(3*an*na(2,ii)+2*bn*na(3,ii)),3*an*na(2,ii)+2*(an+bn)*na(3,ii)+3*bn*na(4,ii),-(2*an*na(3,ii)+3*bn*na(4,ii)),0; 0,-(2*an*na(3,ii)+3*bn*na(4,ii)),2*an*na(3,ii)+(an+3*bn)*na(4,ii)+4*bn*na(5,ii),-(an*na(4,ii)+4*bn*na(5,ii)); 0,0,-(an*na(4,ii)+4*bn*na(5,ii)),an*na(4,ii)+4*bn*na(5,ii)]/NK; Rvec_n(:,ii)=sqrtm(dt*Dmtx)*rand(4,1); end n(2:5,:)=n(2:5,:)+dt*trans_n+Rvec_n; n(1,:)=ones(1,nsim)-sum(n(2:5,:),1); end printf("time %g ms\n",t) end time=toc() scf(0); clf plot(dt:dt:Tstop,Norec) scf(1); clf plot(dt:dt:Tstop,[mean(Norec,2),variance(Norec,2)]) scf(2); clf plot(mean(Norec,2),variance(Norec,2)) printf("time = %g\n",time);