- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2021-09-06 10:31:23+12:00
- Desc:
- Bug fixes
- Permanent Source URI:
- https://models.physiomeproject.org/workspace/6cc/rawfile/da2b6203aca96bb233245d94ca870d00f6ea09d0/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.
% return W from kinetic_parameters
clear;
clc;
close all;
%% booleans
write_parameters_file = true;
include_constraints = true;
include_type2_reactions = 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];
if contains(current_dir, 'beta1') && false
matstr = '_withR_LR_scheme4';
else
matstr = '';
end
%% Define volumes
V_myo = 34.4; % pL
V = struct();
V.V_myo = V_myo;
%% Load forward matrix
if include_type2_reactions
stoich_path = [data_dir sprintf('all_forward_matrix%s.txt',matstr)];
else
stoich_path = [data_dir '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;
dims = struct();
dims.num_cols = num_cols;
dims.num_rows = num_rows;
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 sprintf('all_reverse_matrix%s.txt',matstr)];
else
stoich_path = [data_dir '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_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, W] = kinetic_parameters(M, include_type2_reactions, dims, V); %%%%%%%%%%%%%%%%
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]);
lambdaW = exp(pinv(M) * log(k));
% Check that kinetic parameters are reproduced by bond graph parameters
k_sub = exp(M*log(lambdaW));
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
lambda = lambdaW./W;
kappa = lambda(1:num_cols);
K = lambda(num_cols+1:end);
fID = fopen([data_dir 'rxnID.txt'], 'r');
rxnID = textscan(fID,'%s', 'delimiter','\n');
fclose(fID);
fID = fopen([data_dir 'Kname.txt'], '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)
