- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2022-02-24 10:32:47+13:00
- Desc:
- Init
- Permanent Source URI:
- https://models.cellml.org/workspace/82d/rawfile/fbfa69a8b36b675965eeab93973da6276c42496a/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_Ks_channel = 5369/N_A*1e15; % channel density: unit fmol (From Pan K channel)
%% Load stoichiometric matrices
N_f = [1 0 0 0 0; % includes GHK
0 0 0 0 0;
0 1 0 1 0;
0 0 0 0 1;
0 0 1 0 0;
0 0 0 0 0;
];
N_r = [0 0 0 0 0;
1 0 0 0 0;
0 0 0 0 0;
0 1 0 0 0;
0 0 0 1 0;
0 0 1 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 'ks_G_GHK.mat']);
P_k_ks = G_GHK/F * 1e12; % Unit pL/s
%% Set up the vectors
load([storage_dir 'ks_xs1_parameters.mat']);
params_xs1 = params_vec;
load([storage_dir 'ks_xs2_parameters.mat']);
params_xs2 = params_vec;
% unit s^-1
alpha_xs1 = params_xs1(1); % unit s^-1
beta_xs1 = params_xs1(3); % unit s^-1
alpha_xs2 = params_xs2(1); % unit s^-1
beta_xs2 = params_xs2(3); % unit s^-1
kf_Ks = [P_k_ks/x_Ks_channel; ... % R_GHK
alpha_xs1; ... % Rx1_0
alpha_xs1; ... % Rx1_1
alpha_xs2; ... % Rx2_0
alpha_xs2]; % Rx2_1
kr_Ks = [P_k_ks/x_Ks_channel; ... % R_GHK
beta_xs1; ... % Rx1_0
beta_xs1; ... % Rx1_1
beta_xs2; ... % Rx2_0
beta_xs2]; % Rx2_1
% Overall
kf = kf_Ks;
kr = kr_Ks;
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
save([output_dir 'K_kappa_1.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
