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