Location: BG_Inhib1 @ 994fd1a745dc / MATLAB / find_BG_parameters.m

Author:
Shelley Fong <s.fong@auckland.ac.nz>
Date:
2021-08-05 09:45:51+12:00
Desc:
Initial commit
Permanent Source URI:
https://models.physiomeproject.org/workspace/6d6/rawfile/994fd1a745dc1c37db281e5e010a534eb88b04f6/MATLAB/find_BG_parameters.m
Compare... Download Source 

% 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.

% 5 AUG PLB separated into PLB and Inhib1 submodules
% (parameter error was large with both submodules together)

% 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;

%% 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];


%% Load forward matrix
stoich_path = [data_dir '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


%% Load reverse matrix
stoich_path = [data_dir '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 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, 1, 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

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 sprintf('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)