% Kinetic rate constants for composite model of FCU-AC 
% they must be computed together because of the presence of new potential closed loops

% order: {'cAMP','LRGbinding_B1AR','LRGbinding_M2','GsProtein','GiProtein'}

function [k_kinetic, N_cT, K_C, W] = grand_kinetic_parameters(M, include_type2_reactions, dims, V, noLRG)

    num_cols = dims.num_cols;
    num_rows = dims.num_rows;

    % shared large/small constants
    bigNum = 1e5;
    fastKineticConstant = bigNum;
    smallReverse = 1e0;
    
% M2 composite: LRG_M2 + GiProtein
    
    K_M2sig1 = 30;    % uM  Kc
    K_M2sig2 = 0.16;  % uM  Kh
    K_M2sig4 = 11;    % uM  Kl
    k_M2sig1p = fastKineticConstant;
    k_M2sig1m = k_M2sig1p*K_M2sig1;
    k_M2sig2p = fastKineticConstant;
    k_M2sig2m = k_M2sig2p*K_M2sig2;
    k_M2sig4p = fastKineticConstant;
    k_M2sig4m = k_M2sig4p*K_M2sig4;
    k_M2sig3p = fastKineticConstant;
    k_M2sig3m = k_M2sig1m*k_M2sig2m*k_M2sig3p*k_M2sig4p/(k_M2sig1p*k_M2sig2p*k_M2sig4m); % interim?
    
    k_M2Act1p = 2.5;                   % 1/s
    k_M2Act1m = smallReverse;          % 1/s            
    k_M2Act2p = 0.05;                  % 1/s
    k_M2Act2m = smallReverse;          % 1/s            
    k_M2Hydp = 0.8;                    % 1/s      
    k_M2Hydm = smallReverse;           % 1/s
    k_M2Reassocp = 1.21e3;             % 1/uM.s
    k_M2Reassocm = k_M2Reassocp/fastKineticConstant;   % 1/s
    
    % Loop1: find kact1m using detailed balance
    k_M2Act1m = (k_M2Act1p*k_M2sig1p)/k_M2sig1m;
    % Loop2: find kact1m using detailed balance
    k_M2Act2m = (k_M2Act2p*k_M2sig3p)/k_M2sig3m;

    % Big Loop: find k_M2sig3p or kact2m using detailed balance
    % *** unnecessary ***
    k_M2sig3m = k_M2Act1p*k_M2sig2m*k_M2sig3p*k_M2sig4p/(k_M2Act1m*k_M2sig2p*k_M2sig4m);
    
    % CLOSED LOOP involving G - aGTP - aGDP - G
    % use detailed balance to find kReasocm with either Act (as they have
    % same equilibrium constant
    if true
        k_M2Reassocm = k_M2Act1p*k_M2Hydp*k_M2Reassocp/(k_M2Act1m*k_M2Hydm);
    end
    
% TOTAL    
    %  order: {'LRGbinding_M2','GiProtein'}
% forward first, then reverse!
    k_kinetic = [...
        k_M2sig1p, k_M2sig2p, k_M2sig3p, k_M2sig4p...% LRG_M2
        k_M2Act1p, k_M2Act2p, k_M2Hydp, k_M2Reassocp... %GiProtein
        k_M2sig1m, k_M2sig2m, k_M2sig3m, k_M2sig4m...
        k_M2Act1m, k_M2Act2m, k_M2Hydm, k_M2Reassocm
        ]';
%     k_kinetic = [...
%         k_M2sig1p, k_M2sig2p, k_M2sig3p, k_M2sig4p...% LRG_M2
%         k_M2sig1m, k_M2sig2m, k_M2sig3m, k_M2sig4m...
%         ]';
%     k_kinetic = [...
%         k_M2Act1p, k_M2Act2p, k_M2Hydp, k_M2Reassocp... %GiProtein
%         k_M2Act1m, k_M2Act2m, k_M2Hydm, k_M2Reassocm
%         ]';
    
    % CONSTRAINTS
    N_cT = [];
    if false
        N_cT = zeros(1,size(M,2)); 
        % substrate B1ARp is in eqlm with product B1ARtot
        N_cT(1,num_cols + 1) = 1;
        N_cT(1,num_cols + 2) = -1;
    end
    K_C = ones(1,size(N_cT,1));

    % volume vector
    W = [ones(num_cols,1); V.V_myo*ones(num_rows,1)];

return 