Location: The cardiac Na+ /K+ ATPase: An updated, thermodynamically consistent model @ 6baaa1d0bfa1 / Code / MATLAB / code / bg_parameters.m

Author:

Michael Pan <mpan92@gmail.com>

Date:

2020-08-27 08:39:13+10:00

Desc:

Fixed units

Permanent Source URI:

https://models.physiomeproject.org/workspace/4ce/rawfile/6baaa1d0bfa1a0772c316ab59631b619df2601f4/Code/MATLAB/code/bg_parameters.m

% This script loads the kinetic parameters from fitting the Na/K ATPase
% kinetic model, and solves for the bond graph parameters.

%clear;
clc;

%% Directory names
current_dir = cd;
Idx_backslash = find(current_dir == filesep);
main_dir = current_dir(1:Idx_backslash(end));
data_dir = [main_dir 'data' filesep];
output_dir = [main_dir 'output' filesep];
storage_dir = [main_dir 'storage' filesep];

%% Load forward matrix
stoich_path = [data_dir 'forward_matrix.txt'];
stoich_file_id = fopen(stoich_path,'r');

stoich_file_data = textscan(stoich_file_id,'%s','delimiter','\n');
fclose(stoich_file_id);

num_rows = length(stoich_file_data{1});
num_cols = sum(stoich_file_data{1}{1} == ',')+1;

N_f = zeros(num_rows,num_cols);

for i_row = 1:num_rows
    line_str = stoich_file_data{1}{i_row};
    line_split = regexp(line_str,',','split');
    
    for i_col = 1:num_cols
        N_f(i_row,i_col) = str2double(line_split{i_col});
    end
end

N_fT = transpose(N_f);

%% Load reverse matrix
stoich_path = [data_dir 'reverse_matrix.txt'];
stoich_file_id = fopen(stoich_path,'r');

stoich_file_data = textscan(stoich_file_id,'%s','delimiter','\n');
fclose(stoich_file_id);

num_rows = length(stoich_file_data{1});
num_cols = sum(stoich_file_data{1}{1} == ',')+1;

N_r = zeros(num_rows,num_cols);

for i_row = 1:num_rows
    line_str = stoich_file_data{1}{i_row};
    line_split = regexp(line_str,',','split');
    
    for i_col = 1:num_cols
        N_r(i_row,i_col) = str2double(line_split{i_col});
    end
end

N_rT = transpose(N_r);

%% Calculate stoichiometric matrix
N = N_r - N_f;
N_T = N_rT - N_fT;

I = eye(num_cols);

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

% syms k12p k23p k22ap k34p k33ap k45p k56p k67p k78p k89p k910p k101p kATPp;
% syms k12m k23m k22am k34m k33am k45m k56m k67m k78m k89m k910m k101m pATPm;
% 
% k = transpose([k12p k23p k22ap k34p k33ap k45p k56p k67p k78p k89p k910p k101p kATPp k12m k23m k22am k34m k33am k45m k56m k67m k78m k89m k910m k101m pATPm]);
% 
% Mk = [M log(k)];
% Mk_rref = rref(Mk(1:M_rank,:));

M_T = transpose(M);
M_T_rref = rref(M_T);

% row = N_T(1,:) + N_T(2,:) + N_T(4,:) + N_T(6,:) + N_T(7,:) + N_T(8,:) + N_T(9,:) + N_T(10,:) + N_T(11,:) + N_T(12,:) - N_T(13,:);

%% Define volumes using data from Luo and Rudy (unit pL)
W_i = 38;
W_e = 5.182;
W_total = W_i + W_e;

%% Calculate the bond graph constants 
load([storage_dir 'kinetic_fitting_results.mat']);
params = params_vec;

fast_kinetic_constant = 1e6;

K_d_Nai = params(3);
k_3b_p = fast_kinetic_constant;
k_3b_m = fast_kinetic_constant*(0.5*K_d_Nai);

k_4b_p = fast_kinetic_constant;
k_4b_m = fast_kinetic_constant*(2*K_d_Nai);

K_d_Nai_0 = params(1);
k_5b_m_0 = fast_kinetic_constant;
k_5b_p_0 = fast_kinetic_constant/K_d_Nai_0;

k_6b_p = params(7);
k_6b_m = k_6b_p/6.3;

k_7b_p = params(8);
k_7b_m = params(11);

K_d_Nae_0 = params(2);
k_8b_m_0 = fast_kinetic_constant;
k_8b_p_0 = fast_kinetic_constant*K_d_Nae_0;

K_d_Nae = params(4);
k_9b_m = fast_kinetic_constant;
k_9b_p = fast_kinetic_constant*2*K_d_Nae;

k_10b_m = fast_kinetic_constant;
k_10b_p = fast_kinetic_constant*0.5*K_d_Nae;

K_d_Ke = params(5);
k_11b_p = fast_kinetic_constant;
k_11b_m = fast_kinetic_constant*0.5*K_d_Ke;

k_12b_p = fast_kinetic_constant;
k_12b_m = fast_kinetic_constant*2*K_d_Ke;

k_13b_p = params(9);
k_13b_m = params(12);

K_d_MgATP = params(6);
k_14b_p = fast_kinetic_constant;
k_14b_m = fast_kinetic_constant*K_d_MgATP;

k_15b_p = params(10);
k_15b_m = params(13);

% Calculate remaining parameter
G_MgATP_0 = 11900;
RT = 8.314*310;
K_MgATP = exp(-G_MgATP_0/RT)*10^6;
K_d_Ki = sqrt((k_3b_m*k_4b_m*k_5b_m_0*k_6b_m*k_7b_m*k_8b_m_0*k_9b_m*k_10b_m*k_11b_m*k_12b_m*k_13b_m*k_14b_m*k_15b_m*K_MgATP)/(k_3b_p*k_4b_p*k_5b_p_0*k_6b_p*k_7b_p*k_8b_p_0*k_9b_p*k_10b_p*k_11b_p*k_12b_p*k_13b_p*k_14b_p*k_15b_p));
% K_d_Ki = 255.13;

k_1b_m = fast_kinetic_constant;
k_1b_p = fast_kinetic_constant*2*K_d_Ki;

k_2b_m = fast_kinetic_constant;
k_2b_p = fast_kinetic_constant*0.5*K_d_Ki;

k = transpose([k_1b_p k_2b_p k_3b_p k_4b_p k_5b_p_0 k_6b_p k_7b_p k_8b_p_0 ...
    k_9b_p k_10b_p k_11b_p k_12b_p k_13b_p k_14b_p k_15b_p ...
    k_1b_m k_2b_m k_3b_m k_4b_m k_5b_m_0 k_6b_m k_7b_m k_8b_m_0 ...
    k_9b_m k_10b_m k_11b_m k_12b_m k_13b_m k_14b_m k_15b_m]);

% Check if the thermodynamic condition holds (not enforced by M at the
% moment)
K_MgATP_est = prod(k(1:15))/prod(k(16:30));
G_MgATP_0_est = -RT*log(K_MgATP_est/1e6);

%% Add extra constraints for sodium, potassium and MgATP hydrolysis
K_Na = 1;
M_Na = zeros(1,size(M,2));
M_Na(33) = 1;
M_Na(34) = -1;

K_K = 1;
M_K = zeros(1,size(M,2));
M_K(31) = 1;
M_K(32) = -1;

M_MgATP = zeros(1,size(M,2));
M_MgATP(35) = 1;
M_MgATP(36) = -1;
M_MgATP(37) = -1;
M_MgATP(38) = -1;

k_mod = [k; K_Na; K_K; K_MgATP];
M_mod = [M; M_Na; M_K; M_MgATP];

    
%% Calculate bond graph constants
% Note: units of kappa are mmol/s, units of K are mmol^-1
lambdaW = exp(pinv(M_mod) * log(k_mod));
W = [ones(30,1); W_i; W_e; W_i; W_e; W_i; W_i; W_i; W_i];

lambda = lambdaW./W;

k_check = exp(M*log(lambdaW));
kappa = lambda(1:num_cols);
K = lambda(num_cols+1:end);

diff = sum(abs((k_mod - exp(M_mod*log(lambdaW)))./k_mod));

%% Print variable assignments for Octave
kappa_names = {'kappa_1';'kappa_2';'kappa_3';'kappa_4';'kappa_5';'kappa_6';...
    'kappa_7';'kappa_8';'kappa_9';'kappa_10';'kappa_11';'kappa_12';'kappa_13';...
    'kappa_14';'kappa_15'};
K_names = {'k_1';'k_2';'k_3';'k_4';'k_5';'k_6';'k_7';'k_8';'k_9';'k_10';...
    'k_11';'k_12';'k_13';'k_14';'k_15';'k_ki';'k_ke';'k_nai';'k_nae';...
    'k_mgatp';'k_mgadp';'k_p';'k_h'};

for i = 1:length(kappa_names)
    fprintf([kappa_names{i} ' = ' num2str(kappa(i)) ';\n']);
end

for i = 1:length(K_names)
    fprintf([K_names{i} ' = ' num2str(K(i)) ';\n']);
end

%% Print LaTeX output for table
fprintf('\nKinetic parameters:\n');
tabular_kinetics = cell(15,2);
tabular_kinetics{1,1} = '$k_1^+$';
tabular_kinetics{1,2} = ['$' num2str(k_6b_p) '\\si{s^{-1}}$'];
tabular_kinetics{2,1} = '$k_1^-$';
tabular_kinetics{2,2} = ['$' num2str(k_6b_m) '\\si{s^{-1}}$'];
tabular_kinetics{3,1} = '$k_2^+$';
tabular_kinetics{3,2} = ['$' num2str(k_7b_p) '\\si{s^{-1}}$'];
tabular_kinetics{4,1} = '$k_2^-$';
tabular_kinetics{4,2} = ['$' num2str(k_7b_m) '\\si{s^{-1}}$'];
tabular_kinetics{5,1} = '$k_3^+$';
tabular_kinetics{5,2} = ['$' num2str(k_13b_p) '\\si{s^{-1}}$'];
tabular_kinetics{6,1} = '$k_3^-$';
tabular_kinetics{6,2} = ['$' num2str(k_13b_m) '\\si{mM^{-2}s^{-1}}$'];
tabular_kinetics{7,1} = '$k_4^+$';
tabular_kinetics{7,2} = ['$' num2str(k_15b_p) '\\si{s^{-1}}$'];
tabular_kinetics{8,1} = '$k_4^-$';
tabular_kinetics{8,2} = ['$' num2str(k_15b_m/1e6) '\\times 10^{6}\\si{s^{-1}}$'];
tabular_kinetics{9,1} = '$K_{\\text{d,Na}_i}^0$';
tabular_kinetics{9,2} = ['$' num2str(K_d_Nai_0) '\\si{mM}$'];
tabular_kinetics{10,1} = '$K_{\\text{d,Na}_e}^0$';
tabular_kinetics{10,2} = ['$' num2str(K_d_Nae_0) '\\si{mM}$'];
tabular_kinetics{11,1} = '$K_{\\text{d,Na}_i}$';
tabular_kinetics{11,2} = ['$' num2str(K_d_Nai) '\\si{mM}$'];
tabular_kinetics{12,1} = '$K_{\\text{d,Na}_e}$';
tabular_kinetics{12,2} = ['$' num2str(K_d_Nae) '\\si{mM}$'];
tabular_kinetics{13,1} = '$K_{\\text{d,K}_i}$';
tabular_kinetics{13,2} = ['$' num2str(K_d_Ki) '\\si{mM}$'];
tabular_kinetics{14,1} = '$K_{\\text{d,K}_e}$';
tabular_kinetics{14,2} = ['$' num2str(K_d_Ke) '\\si{mM}$'];
tabular_kinetics{15,1} = '$K_{\\text{d,MgATP}}$';
tabular_kinetics{15,2} = ['$' num2str(K_d_MgATP) '\\si{mM}$'];

for i_line = 1:size(tabular_kinetics,1)
    fprintf([tabular_kinetics{i_line,1} ' & ' tabular_kinetics{i_line,2} ...
         '\\\\ \n']);
end

fprintf('\nBond graph parameters:\n');
reaction_names = {'R1';'R2';'R3';'R4';'R5';'R6';'R7';'R8';'R9';'R10'; ...
    'R11';'R12';'R13';'R14';'R15'};
kappa_var_names = {'\\kappa_1';'\\kappa_2';'\\kappa_3';'\\kappa_4';...
    '\\kappa_5';'\\kappa_6';'\\kappa_7';'\\kappa_8';'\\kappa_9';...
    '\\kappa_{10}';'\\kappa_{11}';'\\kappa_{12}';'\\kappa_{13}';...
    '\\kappa_{14}';'\\kappa_{15}'};
species_names = {'\\text{P}_1';'\\text{P}_2';'\\text{P}_3';'\\text{P}_4';...
    '\\text{P}_5';'\\text{P}_6';'\\text{P}_7';'\\text{P}_8';'\\text{P}_9'; ...
    '\\text{P}_{10}';'\\text{P}_{11}';'\\text{P}_{12}';'\\text{P}_{13}';...
    '\\text{P}_{14}';'\\text{P}_{15}';'\\text{K}_\\text{i}^+';'\\text{K}_\\text{e}^+';...
    '\\text{Na}_\\text{i}^+';'\\text{Na}_\\text{e}^+'; '\\text{MgATP}';...
    '\\text{MgADP}';'\\text{P}_\\text{i}';'\\text{H}^+'};
K_var_names = {'K_1';'K_2';'K_3';'K_4';'K_5';'K_6';'K_7';'K_8';'K_9';'K_{10}';...
    'K_{11}';'K_{12}';'K_{13}';'K_{14}';'K_{15}';'K_\\text{Ki}';'K_\\text{Ke}';...
    'K_\\text{Nai}';'K_\\text{Nae}';'K_\\text{MgATP}';'K_\\text{MgADP}';...
    'K_\\text{Pi}';'K_\\text{H}'};

for i = 1:length(kappa_names)
    fprintf([reaction_names{i} ' & $' kappa_var_names{i} '$ & ' num2str(kappa(i)) '\\\\ \n' ]);
end

for i = 1:length(K_names)
    fprintf(['$' species_names{i} '$ & $' K_var_names{i} '$ & ' num2str(K(i)) '\\\\ \n' ]);
end

%% Save results
save([storage_dir 'NaK_bg_params.mat'],'kappa','K','params');