function [nNaT,nKT,MnNa,MnK,C0Na,C0K,QNa,QK,mINa,mIK,sINa,sIK] = simulations_CovN_full_MC()

% coded by Roberto Fernández Galán, rfgalan@case.edu, December 2011

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   gammaNa  :   single Na channel conductance  (20 pS)
%
%   gammaK   :   single Na channel conductance  (20 pS)
%
%   memS     :   membrane patch size (example: 100 um^2)
%
%   v        :   holding potential in voltage-clamp (example value: -10 mV)
%
%   T        :   time window for integration (example value: 1000 ms)
%
%   trans    :   initial interval with nonstationary dynamics to be removed
%              (should not be smaller than 50 ms and must be larger than T)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% parameters definition

dt = 0.01;
T = 2500;
trans = 500;

v = -38;
memS = 100;

N = round(T/dt);
Ntrans = round(trans/dt);
time = (0:1:N-1)*dt;

ENa = 50;
EK = -77;

gammaNa = 20e-9; % in mS per channel
gammaK = 20e-9;  % in mS per channel

% gammaNa = gammaNa*10^6; %so that the current is in nA
% gammaK = gammaK*10^6; %so that the current is in nA

% input parameters check

if trans >  T
    
    disp('T must be larger than the nonstationary transient\n')
    disp('Setting transient to 10%T\n')

    trans = round(T/trans);
    
end

% initializing variables

NNa = 250*memS; 
NK = 50*memS; 

% homogeneous occupancy at the beginning of the simulations

    temp = NNa/8*ones(8,1); 
    nNa = round(temp);
    [~,k] = max(nNa);
    nNa(k) = nNa(k) - (sum(nNa) - NNa);

    temp = NK/5*ones(5,1);
    nK = round(temp);
    [~,k] = max(nK);
    nK(k) = nK(k) - (sum(nK) - NK);

nNaT = zeros(8,T);
nKT = zeros(5,T);

nNaT(:,1) = nNa;
nKT(:,1) = nK;

% time integration

for t = 1:N-1

    ma = malfa(v);
    mb = mbeta(v);
    
    ha = halfa(v);
    hb = hbeta(v);
    
    na = nalfa(v);
    nb = nbeta(v);
    
    % stochastic gating of sodium channels
    
    check = 0;
    
    while check == 0
          
      PNaij = [0     sum(rand(1,nNa(1)) < 3*ma*dt)    0       0      sum(rand(1,nNa(1)) < ha*dt)      0       0      0;...
               sum(rand(1,nNa(2)) < mb*dt)            0       sum(rand(1,nNa(2)) < 2*ma*dt)    0      0      sum(rand(1,nNa(2)) < ha*dt)      0      0;...
               0     sum(rand(1,nNa(3)) < 2*mb*dt)    0       sum(rand(1,nNa(3)) < ma*dt)     0       0      sum(rand(1,nNa(3)) < ha*dt)      0;...
               0     0       sum(rand(1,nNa(4)) < 3*mb*dt)    0      0       0       0      sum(rand(1,nNa(4)) < ha*dt);...
               sum(rand(1,nNa(5)) < hb*dt)    0       0       0      0       sum(rand(1,nNa(5)) < 3*ma*dt)    0      0;...
               0     sum(rand(1,nNa(6)) < hb*dt)      0       0      sum(rand(1,nNa(6)) < mb*dt)      0       sum(rand(1,nNa(6)) < 2*ma*dt)   0;...
               0     0       sum(rand(1,nNa(7)) < hb*dt)      0      0       sum(rand(1,nNa(7)) < 2*mb*dt)    0      sum(rand(1,nNa(7)) < ma*dt);...
               0     0       0       sum(rand(1,nNa(8)) < hb*dt)     0       0       sum(rand(1,nNa(8)) < 3*mb*dt)   0];  
       
      dnNa = (PNaij' - PNaij)*ones(8,1);

      if sum(nNa + dnNa < 0) > 0
       
           check = 0;

       else

           check = 1;
           nNa = nNa + dnNa;

      end
       
      nNaT(:,t+1) = nNa;
      
    end
    
   % stochastic gating of potassium channels
    
   check = 0;
    
   while check == 0
          
      PKij = [0      sum(rand(1,nK(1)) < 4*na*dt)    0      0      0;...
              sum(rand(1,nK(2)) < nb*dt)             0      sum(rand(1,nK(2)) < 3*na*dt)    0      0;...
              0      sum(rand(1,nK(3)) < 2*nb*dt)    0      sum(rand(1,nK(3)) < 2*na*dt)    0;...
              0      0      sum(rand(1,nK(4)) < 3*nb*dt)    0      sum(rand(1,nK(4)) < na*dt);...
              0      0       0      sum(rand(1,nK(5)) < 4*nb*dt)   0];
 
       
      dnK = (PKij'-PKij)*ones(5,1);

      if sum(nK + dnK < 0) > 0
       
           check = 0;

       else

           check = 1;
           nK = nK + dnK;

      end
       
      nKT(:,t+1) = nK; 
      
   end
             
end

%% plotting traces of the occupation numbers

time = time - time(Ntrans+1);

figure
plot(time(Ntrans+1:end),nNaT(:,Ntrans+1:end))

figure
plot(time(Ntrans+1:end),nKT(:,Ntrans+1:end))

%% mean and covariances of the occupation numbers and their change

MnNa = mean(nNaT(:,Ntrans+1:N),2);
MnK = mean(nKT(:,Ntrans+1:N),2);

C0Na = cov(nNaT(:,Ntrans+1:N)');
C0K = cov(nKT(:,Ntrans+1:N)');

QNa = cov(diff(nNaT(:,Ntrans+1:N),1,2)');
QK = cov(diff(nKT(:,Ntrans+1:N),1,2)');

%% mean and standard deviation of the currents (in uA)

mINa = gammaNa*MnNa(8)*(v-ENa);
mIK = gammaK*NK*MnK(5)*(v-EK);

sINa = gammaNa*sqrt(C0Na(8,8))*(v-ENa);
sIK = gammaK*sqrt(C0K(5,5))*(v-EK);

%% saving results

save simulations_CovN_full_MC.mat

%% plotting means and covariances

figure
subplot(1,2,1)
bar(MnNa)
title('mean state occupancy for Na')
subplot(1,2,2)
bar(MnK)
title('mean state occupancy for K')

figure
subplot(2,2,1)
imagesc(QNa)
colorbar
title('Q_{Na}')

subplot(2,2,2)
imagesc(QK)
colorbar
title('Q_K')

subplot(2,2,3)
imagesc(C0Na)
colorbar
title('C0_{Na}')

subplot(2,2,4)
imagesc(C0K)
colorbar
title('CO_{K}')

end

%% functions for gating variables (transition rates)

function ma = malfa(V)
    
    if abs(V+40) < 0.1
		ma = 1;
	else
		ma = 0.1*(V+40)/(1-exp(-(V+40)/10));
    end
    
end

function mb = mbeta(V)

	mb = 4*exp(-(V+65)/18);

end
    
function ha = halfa(V)
	
    ha = 0.07*exp(-(V+65)/20);

end

function hb = hbeta(V)

	hb = 1/(1+exp(-(V+35)/10));

end    

function na = nalfa(V)

    if abs(V+55) < 0.1

		na = 0.1;

	else

		na = 0.01*(V+55)/(1-exp(-0.1*(V+55)));

    end
    
end
    
function nb = nbeta(V)

	nb = 0.125*exp(-(V+65)/80);
    
end