Location: BG_composite_M2 @ a937e7e137dc / MATLAB2 / find_BG_parameters_M2.m

Author:

Shelley Fong <s.fong@auckland.ac.nz>

Date:

2021-09-09 16:29:05+12:00

Desc:

Calling new submodules from cellml

Permanent Source URI:

https://models.cellml.org/workspace/700/rawfile/a937e7e137dc3338f8fe5c4fa5505ca4cbe7f2dd/MATLAB2/find_BG_parameters_M2.m

% find bond-graph parameters for muscarinic module cell with multiple modules
% present:
%  LRGbinding_M2; GiProtein;

% based on Pan cardiac AP 2018 code

clear;
% clc;
close all;

%% Set directories
current_dir = cd;
Idx_backslash = find(current_dir == filesep);
parent_dir = current_dir(1:Idx_backslash(end-1));
main_dir = current_dir(1:Idx_backslash(end));
output_dir = [current_dir filesep 'output' filesep];

%% Define constants
R = 8.314; % unit J/mol/K
T = 310;
F = 96485;

%% Define volumes (unit pL)
V_myo = 34.4;
V_e = 5.182;     % external volume
V_SR = V_myo*0.035; % SR volume
V_di = V_myo*0.0539; % diadic volume
V = struct();
V.V_myo = V_myo;
V.V_SR = V_SR;
V.V_di = V_di;
V.V_e = V_e;

%% Load stoichiometric matrices and kinetic rate constants
array_names = {'LRGbinding'; 'GiProtein'};
num_subsystems = length(array_names);
sys_struct = cell(num_subsystems,1);
rxnIDs = {};
Knames = {};
keepKname = {};
for i_system = 1:num_subsystems
    sys_name = array_names{i_system};
    sys_dir = [current_dir '\BG_' sys_name '\'];
    cd(sys_dir)
    sys_struct{i_system}.name = sys_name;
    forward_mat_path = ['data\all_forward_matrix.txt'];
    reverse_mat_path = ['data\all_reverse_matrix.txt'];
    N_f = load_matrix(forward_mat_path);
    N_r = load_matrix(reverse_mat_path);
    sys_struct{i_system}.N_f = N_f;
    sys_struct{i_system}.N_r = N_r;
    
%     disp(array_names{i_system});
    dims = struct();
    dims.num_rows = size(N_f,1);
    dims.num_cols = size(N_f,2);
    I = eye(dims.num_cols);
    M = [I transpose(N_f); I transpose(N_r)];
    [k_kinetic, N_cT, K_C, W] = kinetic_parameters(M, true, dims, V);
    sys_struct{i_system}.kfkr = k_kinetic;
    sys_struct{i_system}.Kc = K_C;
    sys_struct{i_system}.N_cT = N_cT;
%         sys_struct{i_system}.W = W(dims.num_cols+1:end);
    fID = fopen('data\rxnID.txt', 'r');
    rxnID = textscan(fID,'%s', 'delimiter','\n');
    fclose(fID);
    rxnIDs = [rxnIDs; rxnID{1}];

    fID = fopen('data\Kname.txt', 'r');
    Kname = textscan(fID,'%s', 'delimiter','\n');
    fclose(fID);
    Knames = [Knames; Kname{1}]; 
end

cd(current_dir)
% relations between submodule to whole module

sys_struct{1}.I_vec = 1:6;                  % LRGbinding
sys_struct{2}.I_vec = [2 3 5 4 6 7:9];      % GiProtein
num_rows = sys_struct{size(sys_struct,1)}.I_vec(end);

N_f = [];
N_r = [];

for i_system = 1:num_subsystems
%     disp(array_names{i_system});
    T = calcT(sys_struct{i_system}.I_vec,num_rows);

    sys_struct{i_system}.T = T;
    
    N_f = [N_f T*sys_struct{i_system}.N_f];
    N_r = [N_r T*sys_struct{i_system}.N_r];
end

N_fT = transpose(N_f);
N_rT = transpose(N_r);

N = N_r - N_f;
N_T = N_rT - N_fT;

num_cols = size(N,2);
I = eye(num_cols);

M = [I N_fT; I N_rT];
M_rref = rref(M);

%% Set up the vectors for kinetic rate constants
%  cAMP; LRGbinding; B1AR; PKA; PLB; Inhib1; GsProtein

kf = [];
kr = [];
W = ones(size(N,2),1);
for i_system=1:num_subsystems
    nrx = length(sys_struct{i_system}.kfkr)/2;
%     disp(array_names{i_system});
%     disp(nrx)
    kf = [kf, sys_struct{i_system}.kfkr(1:nrx)'];
    kr = [kr, sys_struct{i_system}.kfkr(nrx+1:end)'];
%     W = [W; sys_struct{i_system}.W]; 
end
k_kinetic = [kf, kr]';
W = [ones(size(N,2),1);V_myo*ones(6,1);V_myo*ones(3,1)];

lambdaW = exp(pinv(M)*log(k_kinetic));
lambda = lambdaW./W;
kappa = lambda(1:size(N,2));
K = lambda(size(N,2)+1:end);

save([output_dir 'whole_saucerman03_BG_parameters.mat'],'kappa','K','k_kinetic');

%% Checks
N_rref = rref(N);
R_mat = null(N,'r');

k_est = exp(M*log(lambdaW));
diff = (k_kinetic - k_est)./k_kinetic;
error = sum(abs(diff));

K_eq = kf./kr;
zero_est = R_mat'*K_eq';


%% print outputs
for ik=1:length(kappa)
    fprintf('var kappa_%s: fmol_per_sec {init: %g, pub: out};\n',rxnIDs{ik},kappa(ik));
end
Kunique = {};
for ik=1:length(Knames)
    if ~any(strcmp(Kunique,Knames{ik}))
        Kunique = [Kunique; Knames{ik}];
    end
end
for ik = 1:length(Kunique)
    fprintf('var K_%s: per_fmol {init: %g, pub: out};\n',Kunique{ik},K(ik));
end
%% FUNCTION SCRIPTS
function matrix = load_matrix(stoich_path)
    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;
    matrix = 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
            matrix(i_row,i_col) = str2double(line_split{i_col});
        end
    end
end


function T = calcT(I_vec,num_rows)
    num_cols = length(I_vec);
    T = zeros(num_rows,num_cols);
    for i = 1:num_cols
        T(I_vec(i),i) = 1;
    end
end