- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2021-06-04 16:32:29+12:00
- Desc:
- New BG parameters
- Permanent Source URI:
- https://models.physiomeproject.org/workspace/674/rawfile/a51cfe49277237dc75e18131fa08ca3089c799c1/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
include_type2_reactions = true;
use_type2_constraints = false;
use_type2_constraints = use_type2_constraints * include_type2_reactions;
%% 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
if include_type2_reactions
stoich_path = [data_dir filesep 'all_forward_matrix.txt'];
else
stoich_path = [data_dir filesep 'all_noType2_forward_matrix.txt'];
end
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
if include_type2_reactions
stoich_path = [data_dir filesep 'all_reverse_matrix.txt'];
else
stoich_path = [data_dir filesep 'all_noType2_reverse_matrix.txt'];
end
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 of kappa + K
num_cols = sum(stoich_file_data{1}{1} == ',')+1; % num of reactions
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 GsAC and Gs/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
% knowns: K_m [=] fmol/L
% kcat [=] per_sec
K1_m = 1.03e12;
k1cat = 0.2;
K2_m = 3.15e+11;
k2cat = 8.5;
K3_m = 8.60e+11;
k3cat = 1e-16;
K4_m = 1.30e+09;
k4cat = 5;
K5_m = 3.00e+10;
K6_m = 4.00e+08;
K7_m = 4.40e+10;
if false
% use k1m == k2p (GUESS. Not a new equilibrium constraint)
% k1p = 2kcat/Km
% use detailed balance to find remaining k2m
% parameter
% for j = 1:4
% end
const = 1e-50;
t1mult = 1e0;
t2mult = 1e8;
% R1 a/b
k1ap = 2*k1cat/K1_m/t1mult;
k1am =const; %k1cat;
k1bp =const; %k1cat;
% R2 a/b
k2ap = 2*k2cat/K2_m/t1mult;
k2am =const; %k2cat;
k2bp =const; %k2cat;
% R3 a/b
k3ap = 2*k3cat/K3_m/t1mult;
k3am =const; %k3cat;
k3bp =const; %k3cat;
% R4 a/b
k4ap = 2*k4cat/K4_m/t1mult;
k4am =const; %k4cat;
k4bp =const; %k4cat;
if include_type2_reactions
% R5
k5p =const; %K5_m/t2mult;
k5m =const; %K5_m/t2mult;
% R6
k6p =const; %K6_m/t2mult;
k6m =const; %K6_m/t2mult;
% R7
k7p =const; %K7_m/t2mult;
k7m =const; %K7_m/t2mult;
end
else
% R1 a/b
k1ap = 1; % 59 dessauer_1997
k1am = k1ap*K1_m - k1cat;
k1bp = 1; %k1cat;
% R2 a/b
k2ap = 1; % iancu
k2am = k2ap*K2_m - k2cat;
k2bp = k2cat;
% R3 a/b
k3ap = 1; % Guess - like Gsa
k3am = k3ap*K3_m - k3cat;
k3bp = k3cat;
% R4 a/b
k4ap = 1e6; % Guess
k4am = k4ap*K4_m - k4cat;
k4bp = k4cat;
if include_type2_reactions
const = 1;
% R5
k5p =100; %K5_m/t2mult;
k5m =100; %K5_m/t2mult;
% R6
k6p =const; %K6_m/t2mult;
k6m =const; %K6_m/t2mult;
% R7
k7p =const; %K7_m/t2mult;
k7m =const; %K7_m/t2mult;
% k5p = 1;
% k5m = K5_m*k5p;
% % R6
% k6p = 1;
% k6m = K6_m*k6p;
% % R7
% k7p = 1;
% k7m = K7_m*k7p;
end
end
% Calculate remaining parameter using detailed balance
k1bm = (k1ap*k1bp)/k1am;
k2bm = (k2ap*k2bp)/k2am;
k3bm = (k3ap*k3bp)/k3am;
k4bm = (k4ap*k4bp)/k4am;
if include_type2_reactions == false
k5p = [];
k6p = [];
k7p = [];
k5m = [];
k6m = [];
k7m = [];
end
%% 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);
if include_type2_reactions == false
kappa = [kappa;0;0;0];
K = [K;100;100;100;100];
end
save([output_dir 'all_params.mat'],'kappa','K','k_kinetic');
rxnID = {'1a','1b','2a','2b','3a','3b','4a','4b','5','6','7'};
% disp('kappa ');
for j = 1:length(kappa)
fprintf('var kappa_%s: fmol_per_sec {init: %g};\n',rxnID{j},kappa(j));
end
% disp('K');
Kname = {'ATP','cAMP','AC','AC_ATP','Gs_AC','Gs_AC_ATP','FSK_AC','FSK_AC_ATP','PDE', ...
'PDE_cAMP','5AMP','IBMX','PDEinh','Gs','FSK'};
for j = 1:length(K)
fprintf('var K_%s: per_fmol {init: %g};\n',Kname{j},K(j));
end
disp(newline)