- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2022-02-28 15:14:59+13:00
- Desc:
- Fix equation finding K Kappa and update parameters
- Permanent Source URI:
- https://models.cellml.org/workspace/82e/rawfile/ffe0f7085fd622917c43914d1045c8f96dd254df/parameter_fitting/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];
%% Define constants and volumes ( [=] pL)
R = 8.314; % unit J/mol/K
T = 310;
F = 96485;
W_i = 38;
W_e = 5.182;
N_A = 6.022e23;
x_to_channel = 5369/N_A*1e15; % channel density: unit fmol. From Pan K channel density
%% Load stoichiometric matrices
N_f = [1 0 0 0 0 0 0 0 0 0 0; % includes GHK for Ki and Ko
0 0 0 0 0 0 0 0 0 0 0;
0 1 0 0 0 0 0 1 0 0 0;
0 0 0 1 0 0 0 0 1 0 0;
0 0 0 0 0 1 0 0 0 1 0;
0 0 0 0 0 0 0 0 0 0 1;
0 0 1 0 0 0 0 0 0 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 0 0 0 0 0 0 0 0 0
];
N_r = [0 0 0 0 0 0 0 0 0 0 0;
1 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 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 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 1 0 0 0 0 0 1 0 0;
0 0 0 0 1 0 0 0 0 1 0;
0 0 0 0 0 0 1 0 0 0 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);
% GHK permeability for K in Ks channel (see PSO_GHK_fitting_curve.m)
load([storage_dir 'to_G_GHK.mat']);
P_K = G_GHK/F * 1e12; % Unit pL/s
%% Set up the vectors
% params_vec: [kf, zf, kr, zr];
load([storage_dir 'to_zdv_parameters.mat']);
params_zdv = params_vec;
load([storage_dir 'to_ydv_parameters.mat']);
params_ydv = params_vec;
% unit s^-1
alpha_zdv_bg = params_zdv(1); % unit s^-1
beta_zdv_bg = params_zdv(3); % unit s^-1
alpha_ydv_bg = params_ydv(1); % unit s^-1
beta_ydv_bg = params_ydv(3); % unit s^-1
z_zf = params_zdv(2);
z_zr = params_zdv(4);
z_yf = params_ydv(2);
z_yr = params_ydv(4);
zf = [z_zf, z_yf];
zr = [z_zr, z_yr];
kf_Na = [P_K/x_to_channel; ... % R_GHK
3*alpha_zdv_bg; ... % Rzdv_00
3*alpha_zdv_bg; ... % Rzdv_01
2*alpha_zdv_bg; ... % Rzdv_10
2*alpha_zdv_bg; ... % Rzdv_11
alpha_zdv_bg; ... % Rzdv_20
alpha_zdv_bg; ... % Rzdv_21
alpha_ydv_bg; ... % Rydv_00
alpha_ydv_bg; ... % Rydv_10
alpha_ydv_bg; ... % Rydv_20
alpha_ydv_bg]; % Rydv_30
kr_Na = [P_K/x_to_channel; ... % R_GHK
beta_zdv_bg; ... % Rzdv_00
beta_zdv_bg; ... % Rzdv_01
2*beta_zdv_bg; ... % Rzdv_10
2*beta_zdv_bg; ... % Rzdv_11
3*beta_zdv_bg; ... % Rzdv_20
3*beta_zdv_bg; ... % Rzdv_21
beta_ydv_bg; ... % Rydv_00
beta_ydv_bg; ... % Rydv_10
beta_ydv_bg; ... % Rydv_20
beta_ydv_bg]; % Rydv_30
% Overall
kf = kf_Na;
kr = kr_Na;
k = [kf;kr];
W = [ones(size(N,2),1);W_i;W_e;ones(8,1)];
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_dir 'K_kappa_to.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
%% display text for cellml
rname = {'to','z00','z01','z10','z11','z20','z21','y00','y10','y20','y30'};
Kname = {'Ki','Ke','S00','S10','S20','S30','S01','S11','S21','S31'};
zname = {'z', 'y'};
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_%sf: dimensionless {init: %g, pub: out};\n', zname{ik}, zf(ik));
end
for ik = 1:length(zr)
fprintf('var z_%sr: dimensionless {init: %g, pub: out};\n', zname{ik}, zr(ik));
end
if true
% component import parameters
for ik = 1:length(kappa)
fprintf('var kappa_%s: fmol_per_sec {pub: in};\n', rname{ik});
end
for ik = 1:length(K)
fprintf('var K_%s: per_fmol {pub: in};\n', Kname{ik});
end
for ik = 1:length(zf)
fprintf('var z_%sf: dimensionless {pub: in};\n', zname{ik});
end
for ik = 1:length(zr)
fprintf('var z_%sr: dimensionless {pub: in};\n', zname{ik});
end
% mapping
for ik = 1:length(kappa)
fprintf('vars kappa_%s and kappa_%s;\n', rname{ik}, rname{ik});
end
for ik = 1:length(K)
fprintf('vars K_%s and K_%s;\n', Kname{ik}, Kname{ik});
end
for ik = 1:length(zf)
fprintf('vars z_%sf and z_%sf;\n', zname{ik}, zname{ik});
end
for ik = 1:length(zr)
fprintf('vars z_%sr and z_%sr;\n', zname{ik}, zname{ik});
end
end
