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

N_fT = transpose(N_f);

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

N_rT = transpose(N_r);

%% Calculate stoichiometric matrix
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 GsAC and Gs/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; 

if false
    % use k1m == k2p (GUESS. Not a new equilibrium constraint)
    % k1p = 2kcat/Km
    % use detailed balance to find remaining k2m
    % parameter
    
%     for j = 1:4
%     end
    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
else
    % R1 a/b
    k1ap = 0.1; % 59 dessauer_1997
    k1am = k1ap*K1_m - k1cat;
    k1bp = k1cat;
    % R2 a/b
    k2ap = 1; % iancu
    k2am = k2ap*K2_m - k2cat;
    k2bp = k2cat;
    % R3 a/b
    k3ap = 1; % Guess - like Gsa
    k3am = k3ap*K3_m - k3cat;
    k3bp = k3cat;
    % R4 a/b
    k4ap = 3e0; % Guess 1e6
    k4am = k4ap*K4_m - k4cat;
    k4bp = k4cat;
    if include_type2_reactions
        const = 1;
        % R5
        k5p =100; %K5_m/t2mult;
        k5m =100; %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;

%         k5p = 1;
%         k5m = K5_m*k5p;
%         % R6
%         k6p = 1;
%         k6m = K6_m*k6p;
%         % R7
%         k7p = 1;
%         k7m = K7_m*k7p;
    end
end

% Calculate remaining parameter using detailed balance
k1bm = (k1ap*k1bp)/k1am;
k2bm = (k2ap*k2bp)/k2am;
k3bm = (k3ap*k3bp)/k3am;
k4bm = (k4ap*k4bp)/k4am;

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));
% Kc = [];
k_kinetic = [k1ap k1bp ...
    k2ap k2bp ...
    k3ap k3bp ...
    k4ap k4bp ...
    k5p ...
    k6p ...
    k7p ...
    k1am k1bm ...
    k2am k2bm ...
    k3am k3bm ...
    k4am k4bm ...
    k5m ...
    k6m ...
    k7m]';
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

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

rxnID = {'1a','1b','2a','2b','3a','3b','4a','4b','5','6','7'};
% disp('kappa ');
for j = 1:length(kappa)
    fprintf('var kappa_%s: fmol_per_sec {init: %g};\n',rxnID{j},kappa(j));
end
% disp('K');
Kname = {'ATP','cAMP','AC','AC_ATP','Gs_AC','Gs_AC_ATP','FSK_AC','FSK_AC_ATP','PDE', ...
'PDE_cAMP','five_AMP','IBMX','PDEinh','Gs','FSK'};
for j = 1:length(K)
    fprintf('var K_%s: per_fmol {init: %g};\n',Kname{j},K(j));
end

disp(newline)



