Simple and accurate Diffusion Approximation algor. for stochastic ion channels (Orio & Soudry 2012)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:141272
" ... We derived the (Stochastic Differential Equations) SDE explicitly for any given ion channel kinetic scheme. The resulting generic equations were surprisingly simple and interpretable – allowing an easy, transparent and efficient (Diffusion Approximation) DA implementation, avoiding unnecessary approximations. The algorithm was tested in a voltage clamp simulation and in two different current clamp simulations, yielding the same results as (Markov Chains) MC modeling. Also, the simulation efficiency of this DA method demonstrated considerable superiority over MC methods, except when short time steps or low channel numbers were used."
Reference:
1 . Orio P, Soudry D (2012) Simple, fast and accurate implementation of the diffusion approximation algorithm for stochastic ion channels with multiple states. PLoS One 7:e36670 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s):
Channel(s): I Na,t; I K;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; MATLAB; SciLab;
Model Concept(s): Action Potentials; Markov-type model; Stochastic simulation;
Implementer(s): Orio, Patricio [patricio.orio at uv.cl];
Search NeuronDB for information about:  I Na,t; I K;
% 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

nsim=200;
Tstop=6; dt=0.005; %ms
points = round(Tstop/dt);
NK=300;

Vhold=-90;
Vtest=70;

tic();

p=1;
Norec = zeros(points,nsim);
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
    fprintf('time %g ms\n',t)
end
time=toc()
figure(1);
clf
plot(dt:dt:Tstop,Norec)

figure(2);
clf
plot(dt:dt:Tstop,[mean(Norec,2),var(Norec,1,2)])

figure(3);
clf
plot(mean(Norec,2),var(Norec,1,2))


fprintf('time = %g\n',time);


Loading data, please wait...