- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2021-05-31 18:20:34+12:00
- Desc:
- increased guess of forward kinetic rate constant in all reactions
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/674/rawfile/b22d4b553a41f8ed56c3e4a790af928d6bbe8213/MATLAB/all_linalg_parameters.m
% This script calculates the bond graph parameters for all reactions of the
% cAMP pathway, based on SERCA model of Pan et al, which is based
% on Tran et al. (2009). Parameters calculated by using the kinetic
% parameters and stoichiometric matrix.
clear;
clc;
close all;
%% booleans
use_type2_constraints = false;
%% Set directories
current_dir = cd;
Idx_backslash = find(current_dir == filesep);
data_dir = [current_dir filesep 'data' filesep];
output_dir = [current_dir filesep 'output' filesep];
%% Load forward matrix
stoich_path = [data_dir filesep 'all_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 filesep 'all_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];
% show substrate ATP is in eqlm with product cAMP
if use_type2_constraints
N_cT = zeros(5,size(M,2)); % one row each per constraint
else
N_cT = zeros(2,size(M,2));
end
N_cT(1,num_cols + 1) = 1;
N_cT(1,num_cols + 2) = -1;
% another constraint for equilibrium between cAMP and 5AMP
N_cT(2,num_cols + 2) = 1;
N_cT(2,num_cols + 11) = -1;
if use_type2_constraints
% equilibrium between PDEinh and PDE/IBMX
N_cT(3,num_cols + 9) = 1;
N_cT(3,num_cols + 12) = 1;
N_cT(3,num_cols + 13) = -1;
% equilibrium between GsaGTPAC and GsaGTP/AC
N_cT(4,num_cols + 3) = 1;
N_cT(4,num_cols + 14) = 1;
N_cT(4,num_cols + 5) = -1;
% equilibrium between FSKAC and FSK/AC
N_cT(5,end) = 1;
N_cT(5,num_cols + 3) = 1;
N_cT(5,num_cols + 7) = -1;
end
M = [M; N_cT];
%% Set the kinetic rate constants
kguess = 1e8;
% R1 a/b
K_m = 1.03e12; % [=] fmol/L
kcat = 0.2; % [=] per_sec
k1ap = K_m; % Guess
k1am = k1ap*K_m - kcat;
k1bp = kcat;
% R2 a/b
K_m = 3.15e+11; % [=] fmol/L
kcat = 8.5; % [=] per_sec
k2ap = K_m; % Guess
k2am = k2ap*K_m - kcat;
k2bp = kcat;
% R3 a/b
K_m = 8.60e+11; % [=] fmol/L
kcat = 1e-16; % [=] per_sec
k3ap = K_m; % Guess
k3am = k3ap*K_m - kcat;
k3bp = kcat;
% R4 a/b
K_m = 1.30e+09; % [=] fmol/L
kcat = 5; % [=] per_sec
k4ap = K_m; % Guess
k4am = k4ap*K_m - kcat;
k4bp = kcat;
% R5
K_m = 3.00e+10; % [=] fmol/L
k5p = K_m;
k5m = K_m*k5p;
% R6
K_m = 4.00e+08; % [=] fmol/L
k6p = K_m;
k6m = K_m*k6p;
% R7
K_m = 4.40e+10; % [=] fmol/L
k7p = K_m;
k7m = K_m*k7p;
% Calculate remaining parameter using detailed balance
k1bm = (k1ap*k1bp)/k1am;
k2bm = (k2ap*k2bp)/k2am;
k3bm = (k3ap*k3bp)/k3am;
k4bm = (k4ap*k4bp)/k4am;
%% Calculate bond graph constants from kinetic parameters
% Note: units of kappa are fmol/s, units of K are fmol^-1
Kc = ones(1,size(N_cT,1));
% Kc = [];
k_kinetic = [k1ap k1bp ...
k2ap k2bp ...
k3ap k3bp ...
k4ap k4bp ...
k5p ...
k6p ...
k7p ...
k1am k1bm ...
k2am k2bm ...
k3am k3bm ...
k4am k4bm ...
k5m ...
k6m ...
k7m]';
k = transpose([k_kinetic' Kc]);
lambda = exp(pinv(M) * log(k));
% Check that kinetic parameters are reproduced by bond graph parameters
k_sub = exp(M*log(lambda));
diff = (k - k_sub)./k;
error = sum(abs(diff));
% Check that there is a detailed balance constraint
Z = transpose(null(transpose(M),'r'));
%% Save bond graph parameters
kappa = lambda(1:num_cols);
K = lambda(num_cols+1:end);
save([output_dir 'all_params.mat'],'kappa','K','k_kinetic');
disp('kappa ');
for j = 1:length(kappa)
fprintf('%g\n',kappa(j));
end
disp('K');
for j = 1:length(K)
fprintf('%g\n',K(j));
end
disp(newline)
