% finding K, kappa for simple mass-action reaction of Ca leaking from SR
% into cytosol.
% Not gated, and not voltage dependent.


W_i = 25.8;
W_e = 3.5182;
z = 1; % for Ca, z is 2  
N_A = 6.022e23;
x_Kp_channel = 725/N_A*1e15; % unit fmol
x_Nab_channel = x_Kp_channel; % number of channels in myocyte [=] fmol
F = 96485;


%% stoichiometric matrices
% include kf and kr in matrices
N_f = [1; 0]; % [Ca_SR, Cai]
N_r = [0; 1];

N_fT = transpose(N_f);
N_rT = transpose(N_r);
N = N_r - N_f;

num_cols = size(N,2); % number of reactions
I = eye(num_cols); 

M = [I N_fT; I N_rT];

% solve lambda
kf = [  channel_k]; 
kr = [  channel_k];


k = [kf;kr];
W = [1;W_i;W_e];

lambdaW = exp(pinv(M)*log(k));
lambda = lambdaW./W;
kappa = lambda(1:size(N,2));
K = lambda(size(N,2)+1:end);


save(['K_kappa.mat'],'kappa','K');

% Checks for kappa and K
N_rref = rref(N);
R_mat = null(N,'r');

K_eq = kf./kr;
zero_est = R_mat'*K_eq;

k_est = exp(M*log(lambdaW));
diff = sum(abs((k-k_est)./k));