% Current through this channel is modelled as a MICHAELIS-MENTEN reaction

F = 96485;
Cm = 153400e-9; % Unit microF

% Define volumes (unit pL)
W_i = 38;
W_e = 5.182;

% stoichiometric matrices
N_f = [1	0;
0	0;
1	0;
0	1;];

N_r = [0	0;
0	1;
0	1;
1	0;];

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

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

M = [I N_fT; I N_rT];

fkc = 1e6;
smr = 1e-3;
Km = 0.5e-3; % [=] mM
Vmax = 1.15; %[=] uA/uF
kcat = Vmax*Cm*1e9/(F*W_i);
E0 = Vmax/kcat;
km = fkc; % guess
kp = (km + kcat)/Km;
kf = [kp; fkc];
kr = [km;smr];
k = [kf;kr];

% Treat the channel and complex as interior (W_i)
W = [ones(size(N,2),1);
    W_i;W_e;W_i;W_i];

lambdaW = exp(pinv(M)*log(k));
lambda = lambdaW./W;
kappa = lambda(1:size(N,2));   % [=] mol/s if corresponding v [=] mol/s
K = lambda(size(N,2)+1:end);   % [=] 1/mol if corresponding q [=] mol

save(['output\K_kappa_pCa.mat'],'kappa','K');

%% Checks
N_rref = rref(N);
R_zdvat = null(N,'r');

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

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

kdiff = [k,k_est];

if diff > 1
    warning('lol DIFF is greater than 1, nublet.');
end
