% check fit of krkf output (alpha_fit) for kf zf kr zr by comparing against
% kinetic (alpha)

clear;



V = (-120:60)'; % unit mV

gate_state_name = {'C3','c2','c1','O','IF','IS'};

R = 8.314;
T = 310;
RT = R*T;
F =  96500;
vol = 2.58E-17; % L
Ko = 4.5;

 P_O_Na =  1.78e-6;
 P_C1 =  1.119e-4;
 P_C2 =  0.0071021;
 P_C3 =  0.8294071;
 P_IF =  0.10199797;
 P_IS =  0.0613826;

alpha = 55.5*exp(0.05547153*(V-12));
beta = 2.357*exp(-0.036588*V);
alpha_in = 2.172*ones(size(V));
beta_in = 1.077*ones(size(V));
alpha_alpha = 65.5*exp(0.05547153*(V-36));
beta_beta = 2.9357*exp(-0.02158*V);
beta_i = 439*exp(-0.02352*(V+25))*4.5/Ko;
alpha_i = 656*exp(0.000942*V).*((4.5/Ko).^ 0.3);
mu = alpha_i.*beta_beta.*alpha_alpha./(alpha_alpha.*beta_i);        

alpha = [alpha alpha_in alpha_alpha alpha_i alpha_alpha]; 
beta = [beta beta_in beta_beta beta_i mu];


% some BG parameters: from fit_alpha_beta_kappa_constant_z.m

krkf = [68.438114	2.1524652
2.172	1.077
21.333536	0.55120363
638.49363	55.998552
21.333536	9414.4161
];

krkf = [krkf(:,1) [1;0;1;0;1] krkf(:,2) [-1;0;-1;-1;4]];

%% test kf, kr, zf, zr for fits against alpha and beta
if true
    
%     % BG
%     zf_exponent = exp(bsxfun(@times, krkf(:,2)', 1e-3*F*V/RT));
%     zr_exponent = exp(bsxfun(@times, krkf(:,4)', 1e-3*F*V/RT));
%     alpha_fit = bsxfun(@times, krkf(:,1)', zf_exponent);
%     beta_fit = bsxfun(@times, krkf(:,3)', zr_exponent);

    for j = 1:5
        zf_exponent = exp(krkf(j,2) * 1e-3*F*V/RT);
        zr_exponent = exp(krkf(j,4) * 1e-3*F*V/RT);
        alpha_fit(:,j) = krkf(j,1) * zf_exponent;
        beta_fit(:,j) = krkf(j,3) * zr_exponent;
    end
    
%     figure, subplot(1,2,1)
%     plot(V, alpha, '.');
%     title('alpha_{OG}'); legend('1','2','3','4','5','6');
%     subplot(1,2,2)
%     plot(V, alpha_fit, '.'); ylim([0 7e4]);
%     title('alpha_{fit}'); legend('1','2','3','4','5','6');
%     
%     figure, subplot(1,2,1)
%     plot(V, beta, '.'); legend('1','2','3','4','5','6');
%     title('beta_{OG}');
%     subplot(1,2,2)
%     plot(V, beta_fit, '.'); ylim([0 14e4]); legend('1','2','3','4','5','6');
%     title('beta_{fit}');
        
    figure,
    for ik = 1:5
        subplot(2,3,ik)
        plot(V, alpha(:,ik), V, alpha_fit(:,ik));
        legend('OG','fit');
        xlabel('V'); ylabel('rate'); title(num2str(ik));
    end
    suptitle('alpha');
    
        figure,
    for ik = 1:5
        subplot(2,3,ik)
        plot(V, beta(:,ik), V, beta_fit(:,ik));
        legend('OG','fit');
        xlabel('V'); ylabel('rate'); title(num2str(ik));
    end
    suptitle('beta');
    
%     % plot Open state only
%     figure,
%     subplot(1,2,1)
%     semilogy(V, alpha(:,4), '.', V, alpha_fit(:,4),'.');
%     title('alpha P_{Open}'); legend('OG','fit'); xlabel('V');
%     
%     subplot(1,2,2)
%     plot(V, beta(:,4), '.', V, beta_fit(:,4),'.');
%     title('beta P_{Open}'); legend('OG','fit'); xlabel('V');
end

%% test v (dPdt) using kappa, K, zf, zr.

% BEWARE: not comparing like and like (OG kinetic: dim'less but BG: mol/s)

if false
        alpha = alpha'; beta = beta';
    
    % declare K_c, q_Sc
    xgg = 2e-19; % mol
%     q_sC3 = 1e-9;
%     q_sC2 = 1e-11;
%     q_sC1 = 1e-11;
%     q_sO = 1e-11;
%     q_sIF = 1e-11;
%     q_sIS = 1e-11;
    q_sC3 = xgg*P_C3;
    q_sC2 = xgg*P_C2;
    q_sC1 = xgg*P_C1;
    q_sO = xgg*P_O_Na;
    q_sIF = xgg*P_IF;
    q_sIS = xgg*P_IS;
    q0 = [q_sC3 q_sC2 q_sC1 q_sO q_sIF q_sIS];
    
    % one K per gate variable
    % K [=] per_(f)mol
    K_C3 =	535.0649335;
    K_C2 =	9.771415327;
    K_C1 =	0.268890314;
    K_O	= 579.0861639;
    K_IF = 116.818253;
    K_IS =	4.82E-07;
    
    % one kappa per reaction
    % kappa [=] (f)mol_per_sec
    kappa_rC3 =	19.61852542;
    kappa_rC2 =	951.5336816;
    kappa_rC1 =	139.9921881;
    kappa_rC1IF = 3.17E-04;
    kappa_rO =	0.07040671;
    kappa_rIF =	0.785664583;
    
    zf = krkf(:,2);
    zr = krkf(:,4);
    
    % state variable potentials
    
    mu_sC3 = RT*log(K_C3*q_sC3);
    mu_sC2 = RT*log(K_C2*q_sC2);
    mu_sC1 = RT*log(K_C1*q_sC1);
    mu_sO = RT*log(K_O*q_sO);
    mu_sIF = RT*log(K_IF*q_sIF);
    mu_sIS = RT*log(K_IS*q_sIS);
    
    Volt = V*1e-3;
    Af_rC3 = mu_sC3+zf(1)*F*Volt;
    Ar_rC3 = mu_sC2+zr(1)*F*Volt;
    Af_rC2 = mu_sC2+zf(2)*F*Volt;
    Ar_rC2 = mu_sC1+zr(2)*F*Volt;
    Af_rC1 = mu_sC1+zf(3)*F*Volt;
    Ar_rC1 = mu_sO+zr(3)*F*Volt;
    Ar_rC1IF = mu_sIF+zr(4)*F*Volt;
    Af_rC1IF = mu_sC1+zf(4)*F*Volt;
    Af_rO = mu_sO+zf(5)*F*Volt;
    Ar_rO = mu_sIF+zr(5)*F*Volt;
    Af_rIF = mu_sIF+zf(6)*F*Volt;
    Ar_rIF = mu_sIS+zr(6)*F*Volt;
    
    v_rC3 = kappa_rC3*(exp(Af_rC3/RT)-exp(Ar_rC3/RT));
    v_rC2 = kappa_rC2*(exp(Af_rC2/RT)-exp(Ar_rC2/RT));
    v_rC1 = kappa_rC1*(exp(Af_rC1/RT)-exp(Ar_rC1/RT));
    v_rC1IF = kappa_rC1IF*(exp(Af_rC1IF/RT)-exp(Ar_rC1IF/RT));
    v_rO = kappa_rO*(exp(Af_rO/RT)-exp(Ar_rO/RT));
    v_rIF = kappa_rIF*(exp(Af_rIF/RT)-exp(Ar_rIF/RT));
    
    % state variable velocities [=] mol/s
    
    %v = [ sc3 sc2 sc1 o if is]
    
    v_bg(:,1) = -v_rC3;
    v_bg(:,2) = v_rC3-v_rC2;
    v_bg(:,3) = v_rC2-v_rC1-v_rC1IF;
    v_bg(:,4)= v_rC1-v_rO;
    v_bg(:,5) = v_rC1IF+v_rO-v_rIF;
    v_bg(:,6) = v_rIF;
    
%     v_bg = v_bg * 1e3 / vol;  % unit: mM/s
    v_bg = 1e-3*v_bg ./ q0;          % unit: 1/s

   % plot rates
   figure, 
   for i = 1:6
       subplot(2,3,i)
       plot(V, v_CC(i,:),V, v_bg(:,i));
       xlabel('V'); ylabel('rate mM/s'); title('CC');
       legend('CC','BG');
       title(gate_state_name{i});
   end

   
end

















