- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2021-08-04 16:31:51+12:00
- Desc:
- Finding BG parameters from kinetic, version 1 (incomplete)
- Permanent Source URI:
- https://models.cellml.org/workspace/6d6/rawfile/21e434d1a778b748dbd26e9e8ad894d2fe47c825/MATLAB/find_BG_parameters.m
% This script calculates the bond graph parameters for all reactions of the
% a given module. Specify the directory.
% based on SERCA model of Pan et al, which is based on Tran et al. (2009).
% Parameters calculated in module's directory, by using the kinetic
% parameters and stoichiometric matrix.
% 19 Jul INITIAL CONDITIONS
% given dq/dt = A.q + c where A is a square matrix, use null(A) to find ICs
% of q.
% A is constructed from linearising equations through lm, and
% contains BG parameters
% K kappa are contained in c, but disregarding this
% as this is not dependent on q.
clear;
clc;
close all;
%% booleans
write_parameters_file = true;
include_constraints = true;
first3Konly = true;
sep_PKACI = true;
%% Set directories
current_dir = pwd;
Idx_backslash = find(current_dir == filesep);
data_dir = [current_dir filesep 'data' filesep];
output_dir = [current_dir filesep 'output' filesep];
%% separate PKACI into PKACI1 and PKACI3 according to rxn
if sep_PKACI
matstr = 'sepPKACI';
else
matstr = 'all';
end
%% Load forward matrix
if ~first3Konly
stoich_path = [data_dir sprintf('%s_forward_matrix.txt',matstr)];
else
stoich_path = [data_dir sprintf('first3_%s_forward_matrix.txt', matstr)];
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 ~first3Konly
stoich_path = [data_dir sprintf('%s_reverse_matrix.txt',matstr)];
else
stoich_path = [data_dir sprintf('first3_%s_reverse_matrix.txt', matstr)];
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_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];
addpath(current_dir);
addpath(data_dir);
[k_kinetic, N_cT, K_C] = kinetic_parameters(M, ~first3Konly, num_cols); %%%%%%%%%%%%%%%%
if ~include_constraints
N_cT = [];
end
M = [M; N_cT];
% Calculate bond graph constants from kinetic parameters
% Note: units of kappa are fmol/s, units of K are fmol^-1
k = transpose([k_kinetic' K_C]);
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'));
% Check this is a thermo consistent nullspace
%Zt*ln(K) = 0
if false
for j = 1:length(k_kinetic)/2
K_eq(j) = k_kinetic((j*2)-1)/k_kinetic(j*2);
end
zspace = transpose(Z)*log([K_eq; K_C]);
end
%% Save bond graph parameters
if ~first3Konly
kappa = lambda(1:num_cols);
K = lambda(num_cols+1:end);
if sep_PKACI
matstr = '_sepPKACI';
else
matstr = '';
end
fID = fopen([data_dir 'rxnID.txt'], 'r');
rxnID = textscan(fID,'%s', 'delimiter','\n');
fclose(fID);
fID = fopen([data_dir sprintf('Kname%s.txt', matstr)], 'r');
Kname = textscan(fID,'%s');
fclose(fID);
Knamed = struct();
for k = 1:length(Kname{1})
% assign name to K values. store as struct for output in .mat.
sname = Kname{1}{k};
Knamed.(sname) = K(k);
end
kappa_named = struct();
for k = 1:length(rxnID{1})
% assign name to K values. store as struct for output in .mat.
sname = rxnID{1}{k};
kappa_named.(['r' sname]) = kappa(k);
end
if write_parameters_file
save([output_dir 'all_params.mat'],'kappa','K','k_kinetic', 'Knamed','kappa_named'); %
end
for j = 1:length(kappa)
fprintf('var kappa_%s: fmol_per_sec {init: %g};\n',rxnID{1}{j},kappa(j));
end
for j = 1:length(K)
fprintf('var K_%s: per_fmol {init: %g};\n',Kname{1}{j},K(j));
end
if false
if false
qname = {'L','LR','LRG','LR_BARK','LRG_BARK','R','RG','B1d','B1p','PKACI_B1p',...
'B1tot','PKACI','BARK','G','GaT','GBy','GaD','GsGDP','Pi'};
qval = [0.0001672850,0.0000016970,0.0001046605,0.0000016970,0.0001046605,0.0002197247,...
0.0000254602,0.0000000000,...
0.0000219260,.0000219260,0.0004579000,0.0022120882,0.0228000000,0.1453052187,...
0.0009519000,0.0009762200,0.0000244948,0.0000000000,570]';
for j=1:length(qname)
fprintf('var q_%s: fmol {init: %g};\n',qname{j}, qval(j));
end
end
for j=1:length(K)
fprintf('var q_%s: fmol {init: 1e-13};\n',Kname{1}{j});
end
for j=1:length(kappa)
fprintf('var v%s: fmol_per_sec;\n',rxnID{1}{j})
end
for j=1:length(K)
fprintf('var mu_%s: J_per_mol;\n',Kname{1}{j});
end
for j=1:length(K)
fprintf('mu_%s = R*T*ln(K_%s*q_%s);\n',Kname{1}{j},Kname{1}{j},Kname{1}{j});
end
for j=1:length(kappa)
fprintf('v%s = kappa_%s*exp(mu_a/(R*T));\n',rxnID{1}{j}, rxnID{1}{j})
end
for j=1:length(K)
fprintf('ode(q_%s, time) = p;\n',Kname{1}{j});
end
end
disp(newline)
end
