Location: BG_Na_background @ 4e37e053bf80 / matlab_parameter_fitting / K_kappa_simplechannel.m

Author:

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

Date:

2024-11-08 11:33:15+13:00

Desc:

Adding exposure files and fixing cellml

Permanent Source URI:

https://models.physiomeproject.org/workspace/82b/rawfile/4e37e053bf80a5eda4b9d523b59348e9bc21db47/matlab_parameter_fitting/K_kappa_simplechannel.m

% finding K, kappa for simple non gated channels like Na_background.

% Steps:
% 1. find alpha and beta rate constants for entry/exit of ion
% 2. use PSO to fit kf and kr to alpha and beta (assume zf and zr = 1) for
% GATES
% 3. Convert G_GHK to Permeability 
% 3. Solve Ln(k) = M.Ln(lambda) where
%             lambda = [kappa ; K ]


V = -120:60;
W_i = 25.8;
W_e = 3.5182;
z = 1; % for Ca, z is 2  
N_A = 6.022e23;
x_Kp_channel = 725/N_A*1e15; % unit fmol
x_Nab_channel = x_Kp_channel; % number of channels in myocyte [=] fmol
F = 96485;

load('nab_G_GHK.mat');
% P_GHK = []; % 1/mM.s

%% 1
alpha = 10;
beta = 0;

%% 2
if false
    error_func_alpha = @(params) square_error(alpha - calc_alpha_simple(params,V/1000,z)); % fit to alpha - k.exp(zFV/RT) : do lsqminerror.
    error_func_beta = @(params) square_error(beta - calc_beta_simple(params,V/1000,z));

    error_func = @(params) error_func_alpha(params) + error_func_beta(params); % + error_func_gss(params) + error_func_tau(params);

    A = [];
    b = [];
    Aeq = [];
    beq = [];
    lb = [0; -10; 0; -10];
    ub = [Inf; 10; Inf; 10];

    options_unc = optimoptions('fminunc','MaxFunEvals',10000);
    options_ps = optimoptions('particleswarm','UseParallel',false,'HybridFcn',@fmincon);

    [params_vec,fval,exitflag,output] = particleswarm(error_func,2,lb,ub,options_ps);
end

%% 3 convert to permeability
P_Nab = G_GHK/F * 1e12; % Unit pL/s
channel_k = P_Nab/x_Nab_channel;

%% 4
% stoichiometric matrices
% include kf and kr in matrices
N_f = [1; 0]; % [Nai, Nao]
N_r = [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];

% solve lambda
% don't need params_vec: only add more reactions if they are present e.g.
% gate. Cab does not have a gate.
kf = [  channel_k]; 
kr = [  channel_k];


k = [kf;kr];
W = [1;W_i;W_e];

lambdaW = exp(pinv(M)*log(k));
lambda = lambdaW./W;
kappa = lambda(1:size(N,2));
K = lambda(size(N,2)+1:end);


save(['K_kappa.mat'],'kappa','K');

% Checks for kappa and K
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));