Location: BG_funny @ 60ceabbf4b0d / Kernik_matlab_parameter_fitting / K_kappa.m

Author:

Shelley Fong <sfon036@UoA.auckland.ac.nz>

Date:

2024-11-06 10:57:15+13:00

Desc:

Fixing text

Permanent Source URI:

https://models.physiomeproject.org/workspace/838/rawfile/60ceabbf4b0df7dd6b375b7d6924f6acad6f6caa/Kernik_matlab_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_f_channel = 5369/N_A*1e15; % channel density: unit fmol (From Pan K channel)

%% Load stoichiometric matrices  

Kname = {'Nai',	'Nae',	'Ki',	'Ke',	'X0',	'X1'};
rname = {'Na', 'K', 'R_xf'};
zname = {'xf'};
N_f = [1	0	0;
0	0	0;
0	1	0;
0	0	0;
0	0	1;
0	0	0;
];

N_r = [0	0	0;
1	0	0;
0	0	0;
0	1	0;
0	0	0;
0	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];
M_rref = rref(M);

% GHK permeability for Na and K in f channel (see PSO_GHK_fitting_curve.m)
load([storage_dir 'fNa_G_GHK.mat']);
P_f_Na = G_GHK/F * 1e12; % Unit pL/s
load([storage_dir 'fK_G_GHK.mat']);
P_f_K = G_GHK/F * 1e12; % Unit pL/s

%% Set up the vectors
load([storage_dir 'f_gate_parameters.mat']);
params_gate = params_vec;

alpha = params_gate(1)*1e3; % unit s^-1
beta = params_gate(3)*1e3; % unit s^-1

kf = [P_f_Na/x_f_channel; ... % R_GHK
    P_f_K/x_f_channel; ... % R_GHK
    alpha];    

kr = [P_f_Na/x_f_channel; ... % R_GHK
    P_f_K/x_f_channel; ... % R_GHK
    beta];    

k = [kf;kr];
W = [ones(size(N,2),1);
    W_i;W_e;ones(4,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.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
zf = params_gate(2);
zr = params_gate(4);

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
    
    % declare and calculate mu, v, ode
    
    % 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