- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2021-07-22 12:36:57+12:00
- Desc:
- Removing unnecessary files
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/674/rawfile/d29401fb021b702a4deac586d4d7ec0eb7717f24/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
write_parameters_file = false;
include_FSK = true;
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
%% 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
if include_FSK == false
fsk_rxn_ids = [5,6,11];
N_f(:,fsk_rxn_ids) = [];
N_r(:,fsk_rxn_ids) = [];
N_f([7,8,15],:) = [];
N_r([7,8,15],:) = [];
num_cols = num_cols - length(fsk_rxn_ids);
end
N_fT = transpose(N_f);
N_rT = transpose(N_r);
%% Calculate stoichiometric matrix
% I matrix to align with placement of kappa down the column.
% x-axis of stoich matrix (R1a, R1b etc) coincides with the kp km of that
% kinetic reaction
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 five_AMP
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 GsTAC and GsT/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;
% initialise arrays
vkap = zeros(4,1);
vkam = zeros(4,1);
vkbp = zeros(4,1);
vkbm = zeros(4,1);
vkm2 = zeros(3,1);
vkp2 = zeros(3,1);
vK_m = [K1_m,K2_m,K3_m,K4_m,K5_m,K6_m,K7_m];
vkcat = [k1cat,k2cat,k3cat,k4cat];
if true
kap_mult = [0.93, 0.1, 1, 1]; % finding rates from literature?
fastKineticConstant = 1e6; %1e-12;
for i=1:4
vkap(i) = kap_mult(i)*fastKineticConstant;
vkam(i) = vkap(i)*vK_m(i) - vkcat(i);
vkbp(i) = vkcat(i);
vkbm(i) = (vkap(i)*vkbp(i))/vkam(i);
end
if include_type2_reactions
for i= 1:3
vkm2(i) = fastKineticConstant;
vkp2(i) = vkm2(i)/vK_m(i+4);
end
end
elseif false
% k1p = 2kcat/Km
% use detailed balance to find remaining k2m parameter
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
% Calculate remaining parameter using detailed balance
k1bm = (k1ap*k1bp)/k1am;
k2bm = (k2ap*k2bp)/k2am;
k3bm = (k3ap*k3bp)/k3am;
k4bm = (k4ap*k4bp)/k4am;
end
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));
if include_FSK
k_kinetic = [vkap(1) vkbp(1) ...
vkap(2) vkbp(2) ...
vkap(3) vkbp(3) ...
vkap(4) vkbp(4) ...
vkp2(:)' ...
vkam(1) vkbm(1) ...
vkam(2) vkbm(2) ...
vkam(3) vkbm(3) ...
vkam(4) vkbm(4) ...
vkm2(:)']';
else
k_kinetic = [vkap(1) vkbp(1) ...
vkap(2) vkbp(2) ...
vkap(4) vkbp(4) ...
vkp2(1:2)' ...
vkam(1) vkbm(1) ...
vkam(2) vkbm(2) ...
vkam(4) vkbm(4) ...
vkm2(1:2)']';
end
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
if write_parameters_file
save([output_dir 'all_params.mat'],'kappa','K','k_kinetic');
end
if include_FSK
rxnID = {'1a','1b','2a','2b','3a','3b','4a','4b','5','6','7'};
Kname = {'ATP','cAMP','AC','AC_ATP','GsT_AC','GsT_AC_ATP','FSK_AC','FSK_AC_ATP','PDE', ...
'PDE_cAMP','five_AMP','IBMX','PDEinh','GsT','FSK'};
else
rxnID = {'1a','1b','2a','2b','4a','4b','5','6'};
Kname = {'ATP','cAMP','AC','AC_ATP','GsT_AC','GsT_AC_ATP','PDE', ...
'PDE_cAMP','five_AMP','IBMX','PDEinh','GsT'};
end
for j = 1:length(kappa)
fprintf('var kappa_%s: fmol_per_sec {init: %g};\n',rxnID{j},kappa(j));
end
for j = 1:length(K)
fprintf('var K_%s: per_fmol {init: %g};\n',Kname{j},K(j));
end
if include_FSK == false
fprintf('var kappa_3a: fmol_per_sec {init: 0};\n');
fprintf('var kappa_3b: fmol_per_sec {init: 0};\n');
fprintf('var kappa_7: fmol_per_sec {init: 0};\n');
fprintf('var K_FSK: per_fmol {init: 1e3};\n');
fprintf('var K_FSK_AC: per_fmol {init: 1e3};\n');
fprintf('var K_FSK_AC_ATP: per_fmol {init: 1e3};\n');
end
disp(newline)
%% check fit against kinetic model for each reaction in pathway
% 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;
vkin_1 = [];
