// Potassium channel from original HH model
// Voltage clamp simulations with non-stationary noise analysis
// Coupled activation particles (5-state channel), Diffusion approximation algorithm

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 33 and 61)
for tt=dt:tint:Tstop //Nested FORs are only for the purpose of reporting the time (see line 61)
    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,:)];

        na=abs(n); //For the stochastic term we use absolute values of states
        // Stochastic terms R1-R4
        R=rand(4,nsim).*sqrt((dt/NK)*[4*an.*na(1,:)+bn.*na(2,:);
        3*an.*na(2,:)+2*bn.*na(3,:);
        2*an.*na(3,:)+3*bn.*na(4,:);
        an.*na(4,:)+4*bn.*na(5,:)]);

        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);