Location: BG_Kr @ bd8f1a808e21 / BG_fit_parameters_CC / K_kappa.m

Author:

Shelley Fong <s.fong@auckland.ac.nz>

Date:

2022-02-24 16:09:46+13:00

Desc:

Add

Permanent Source URI:

https://staging.physiomeproject.org/workspace/82c/rawfile/bd8f1a808e21ed308da7eedca222fd1ae09cd992/BG_fit_parameters_CC/K_kappa.m

% 18 April: make it fit 6 state CC model

clear;

%% Set directories
main_dir = pwd;
data_dir = [main_dir '\data' filesep];
output_dir = [main_dir '\output' filesep];
storage_dir = [main_dir '\storage' filesep];

fixedz = false;

%% Load stoichiometric matrices  

N_f = [1	0	0	0	0	0;
0	0	0	0	0	0;
0	1	0	0	0	0;
0	0	1	0	0	0;
0	0	0	1	1	0;
0	0	0	0	0	1;
0	0	0	0	0	0;
];

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

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

N = N_r - N_f;
% N_T = N_rT - N_fT;

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

M = [I N_fT; I N_rT];
M_rref = rref(M);


%% Set up the vectors
% GHK permeability for K in Ks channel (see PSO_GHK_fitting_curve.m)
F = 96485;
load([storage_dir 'kr_G_GHK.mat']);
P_K = G_GHK/F * 1e12; % Unit pL/s
N_A = 6.022e23;
x_kr_channel = 5369/N_A*1e15; % channel density: unit fmol. Same as Pan K channel density

% load alpha_beta rates
if fixedz
    load([storage_dir 'kr_alphabeta_fixedz.mat']);
else
    load([storage_dir 'kr_alphabetaz_someIC.mat']);
%     load([storage_dir 'kr_alphabetaz_zeroIC.mat']);
end

% unit s^-1
% [rc3 rc2 rc1 rc1IF rOIF]
% kf kr [=] 1/s 
kf = [P_K/x_kr_channel % R_GHK
alpha];

kr = [P_K/x_kr_channel % R_GHK
beta];


k = [kf;kr];
% don't need W: all ones. only use Wi and We for GHK permeabilities : but
% use the P value from Pan.
% W = [ones(size(N,2),1);W_i;W_e;ones(10,1);W_i;W_e;ones(16,1);W_i;W_e;ones(12,1);...
%     ones(15,1);W_i;W_i;W_i;W_i;ones(6,1);W_i;W_i;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

if fixedz
    save([output_dir 'K_kappa_Kr_fixedz.mat'],'kappa','K');
else
    save([output_dir 'K_kappa_Kr_variablez.mat'],'kappa','K');
end

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

kdiff = [k,k_est];

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

%% display text for cellml
rname = {'Kr','rC3', 'rC2', 'rC1', 'rC1IF', 'rOIF'};
Kname = {'Ki','Ke','C3','C2','C1','O','I'};
for ik = 1:length(kappa)
    fprintf('var kappa_%s: fmol_per_sec {init: %g, pub: out};\n', rname{ik}, kappa(ik));
end
for ik = 1:length(K)
    fprintf('var K_%s: per_fmol {init: %g, pub: out};\n', Kname{ik}, K(ik));
end
for ik = 1:length(zf)
    fprintf('var z_f%s: dimensionless {init: %g, pub: out};\n', Kname{ik}, zf(ik));
end
for ik = 1:length(zr)
    fprintf('var z_r%s: dimensionless {init: %g, pub: out};\n', Kname{ik}, zr(ik));
end