% Based on Pan cardiac action potential BG modelling 2018
% only LCC channel is retained.

clear;
clc;
close all;

%% Set directories
current_dir = cd;
Idx_backslash = find(current_dir == filesep);
data_dir = [current_dir filesep 'data' filesep];
code_dir = [current_dir filesep 'code' filesep];
output_dir = [current_dir filesep 'output' filesep];
storage_dir = [current_dir filesep 'storage' filesep];

%% Define constants
R = 8.314; % unit J/mol/K
T = 310;
F = 96485;
% 
N_A = 6.022e23;
A_cap = 1.534e-4; % Unit cm^2

x_LCC = 50000/N_A*1e15; % unit fmol

%% Define volumes (unit pL)
W_i = 38;
W_e = 5.182;

% %% Import IV parameters

% L-type calcium channel
% Load P_GHK (unit pL/s)
load([storage_dir 'LCC_P_GHK_Ca.mat']);
P_LCC_Ca = P_GHK; % Unit pL/s
load([storage_dir 'LCC_P_GHK_K.mat']);
P_LCC_K = P_GHK; % Unit pL/s

% 
load([storage_dir 'LCC_d_parameters.mat']);
params_d = params_vec;

load([storage_dir 'LCC_f_parameters.mat']);
params_f = params_vec;

%% Load stoichiometric matrices
stoich_path = [data_dir 'LCC_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 'LCC_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);

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

M_T = transpose(M);
M_T_rref = rref(M_T);

%% Set up the vectors

% L-type calcium channel
alpha_d0_bg = params_d(1)*1e3; % unit s^-1
beta_d0_bg = params_d(3)*1e3; % unit s^-1

alpha_f1_0_bg = params_f(1)*1e3; % unit s^-1
beta_f1_0_bg = params_f(3)*1e3; % unit s^-1

alpha_f2_0_bg = params_f(5)*1e3; % unit s^-1
beta_f2_0_bg = params_f(7)*1e3; % unit s^-1

K_f3_0 = beta_f1_0_bg*alpha_f2_0_bg/alpha_f1_0_bg/beta_f2_0_bg; % dimensionless
rate_f3 = 1e5; % unit s^-1

KmCa = 0.6e-3; % Unit mM
rate_fCa = 1e5; % unit s^-1

kf_LCC = [P_LCC_Ca/x_LCC; ... % R_GHK_Ca1
    P_LCC_Ca/x_LCC; ... % R_GHK_Ca2
    P_LCC_K/x_LCC; ... % R_GHK_K1
    P_LCC_K/x_LCC; ... % R_GHK_K2
    alpha_d0_bg; ... % Rd000
    alpha_d0_bg; ... % Rd010
    alpha_d0_bg; ... % Rd020
    alpha_d0_bg; ... % Rd001
    alpha_d0_bg; ... % Rd011
    alpha_d0_bg; ... % Rd021
    alpha_f1_0_bg; ... % Rf1_000
    alpha_f1_0_bg; ... % Rf1_100
    alpha_f1_0_bg; ... % Rf1_001
    alpha_f1_0_bg; ... % Rf1_101
    alpha_f2_0_bg; ... % Rf2_000
    alpha_f2_0_bg; ... % Rf2_100
    alpha_f2_0_bg; ... % Rf2_001
    alpha_f2_0_bg; ... % Rf2_101
    K_f3_0*rate_f3; ... % Rf3_010
    K_f3_0*rate_f3; ... % Rf3_110
    K_f3_0*rate_f3; ... % Rf3_011
    K_f3_0*rate_f3; ... % Rf3_111
    rate_fCa; ... % RfCa000
    rate_fCa; ... % RfCa100
    rate_fCa; ... % RfCa010
    rate_fCa; ... % RfCa110
    rate_fCa; ... % RfCa020
    rate_fCa]; % RfCa120

kr_LCC = [P_LCC_Ca/x_LCC; ... % R_GHK_Ca1
    P_LCC_Ca/x_LCC; ... % R_GHK_Ca2
    P_LCC_K/x_LCC; ... % R_GHK_K1
    P_LCC_K/x_LCC; ... % R_GHK_K2
    beta_d0_bg; ... % Rd000
    beta_d0_bg; ... % Rd010
    beta_d0_bg; ... % Rd020
    beta_d0_bg; ... % Rd001
    beta_d0_bg; ... % Rd011
    beta_d0_bg; ... % Rd021
    beta_f1_0_bg; ... % Rf1_000
    beta_f1_0_bg; ... % Rf1_100
    beta_f1_0_bg; ... % Rf1_001
    beta_f1_0_bg; ... % Rf1_101
    beta_f2_0_bg; ... % Rf2_000
    beta_f2_0_bg; ... % Rf2_100
    beta_f2_0_bg; ... % Rf2_001
    beta_f2_0_bg; ... % Rf2_101
    rate_f3; ... % Rf3_010
    rate_f3; ... % Rf3_110
    rate_f3; ... % Rf3_011
    rate_f3; ... % Rf3_111
    rate_fCa/KmCa; ... % RfCa000
    rate_fCa/KmCa; ... % RfCa100
    rate_fCa/KmCa; ... % RfCa010
    rate_fCa/KmCa; ... % RfCa110
    rate_fCa/KmCa; ... % RfCa020
    rate_fCa/KmCa]; % RfCa120

% Overall
k = [kf_LCC;kr_LCC];
W = [ones(size(N,2),1);W_i;W_e;W_i;W_e;ones(12,1)]; % kappa Cai Cae Ki Ke [gates]

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

save([output_dir 'LCC_parameters.mat'],'kappa','K','params_f');

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

K_eq = kf_LCC./kr_LCC;
zero_est = R_mat'*K_eq;

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

%% output for cellML
if false
    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);

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