// Potassium channel from original HH model
// Voltage clamp simulations with non-stationary noise analysis
// Coupled activation particles (5-state channel), Diffusion approximation algorithm
// Steady-state approximation of variables in stochastic terms.

stacksize('max');
nsim=200;  //number of sweeps to be simulated
Tstop=6; dt=0.001; //Total time and dt in ms
points = round(Tstop/dt) //number of points per sweep
NK=300; //number of potassium channels

Vhold=-90; //voltage for t=0
Vtest=70;
rand('normal');

p=1;
Norec = zeros(points,nsim);

v = Vhold*ones(1,nsim); //calculus of equilibrium state at t=0
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);

//Now we change voltage for the rest of the simulation
v = Vtest*ones(1,nsim);
an=0.01*(v+55)./(1-exp(-(v+55)/10));
bn=0.125*exp(-(v+65)/80);

tic();

tint=1; //period for reporting simulation time (see lines 34 and 62)
for tt=dt:tint:Tstop //Nested FORs are only for the purpose of reporting the time (see line 62)
    for t = tt:dt:tt+tint-dt
        
        Norec(p,:) = n(5,:)*NK;
        p=p+1;

        trans_n=[-4*an.*n(1,:)+bn.*n(2,:); //Deterministic part of state transitions
        -(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,:)];
        
        // Stochastic terms R1-R4 do not include variable (n) values but only kinetic constants
        R=rand(4,nsim).*sqrt((dt/NK)*[8*an.*bn.^(4);
        24*an.^(2).*bn.^3;
        24*an.^(3).*bn.^2;
        8*an.^(4).*bn]./(ones(4,1)*(an+bn).^4));

        Wtn=[R(1,:);
        -R(1,:)+R(2,:);
        -R(2,:)+R(3,:);
        -R(3,:)+R(4,:);
        -R(4,:);]
        
        n=n+dt*trans_n+Wtn;
        n(1,:) = ones(1,nsim) - sum(n(2:5,:),1)  //adjusting the sum of states equal to 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);