% GPCR module. From Dupont Sneyd, simple GPCR model.

function [k_kinetic, N_cT, K_C] = kinetic_parameters(M, include_type2_reactions, num_cols)
    % Set the kinetic rate constants.
    % original model had reactions that omitted enzymes as substrates e.g. BARK
    % convert unit from 1/s to 1/uM.s by dividing by conc of enzyme
    % all reactions were irreversible, made reversible by letting kr ~= 0

    
    % CONVERT TO fM 
    bigNum = 1e3;
    fastKineticConstant = bigNum;
    smallReverse = fastKineticConstant/(bigNum^2);
    u_to_f = 1e9;
    c_G = 3.83*u_to_f; % fM
    
    kbindm = fastKineticConstant; % 1/s
    kbindp = kbindm / (0.285*u_to_f); % 1/fM.s
    kActp = 16/c_G;            % 1/fM.s
    kActm = smallReverse;          % 1/fM.s            
    kHydp = 0.8;            % 1/s      
    kHydm = smallReverse;           % 1/s
    kReassocp = 1.21e3/u_to_f;     % 1/fM.s
    kReassocm = kReassocp/bigNum;          % 1/s
    
    % No closed loop, so no detailed balance
    
    k_kinetic = [
        kbindp, kActp, kHydp, kReassocp...
        kbindm, kActm, kHydm, kReassocm
        ]';

    % CONSTRAINTS
    N_cT = zeros(1,size(M,2)); 
    % Reaction BIND: equal portions of L + derivatives (L == LR)    **SMALL_ERROR**
    N_cT(1,num_cols + 1) = -1;
%     N_cT(1,num_cols + 2) = -1;    % this must be present, like Mg reaction in SERCA
    N_cT(1,num_cols + 4) = 1;
    if false
        % Reaction BIND: equal portions of R == LR
        % error 10(2): this reaction is already in use?
        N_cT(2,num_cols + 2) = -1;
        N_cT(2,num_cols + 4) = 1;
    end
    % Reaction HYD: substrate aGTP is in eqlm with product aGDP
    % error 10(2): a-G is present in other reactions
    if false
        N_cT(2,num_cols + 5) = 1;
        N_cT(2,num_cols + 7) = -1;
    end
    % Reaction ACT:
    %  a-GTP == betaGamma                                           **SMALL_ERROR**
    if false                                                          
        N_cT(2,num_cols + 5) = 1;
        N_cT(2,num_cols + 6) = -1;
    end
    % or LR.G = aGTP.betaGamma
    if false                                                          
        N_cT(2,num_cols + 3) = 1;
        N_cT(2,num_cols + 4) = 1;
        N_cT(2,num_cols + 5) = -1;
        N_cT(2,num_cols + 6) = -1;
    end
    % [a-GTP] + [a-GDP] = [beta.gamma]                              **SMALL_ERROR**
    if true % could be acceptable !!!!!!!!!!!!!!!!!!!!!!!!!!
        N_cT(2,num_cols + 6) = 1;   % beta_gamma
        N_cT(2,num_cols + 5) = -1;  % a_GTP
        N_cT(2,num_cols + 7) = -1;  % a_GDP
    end
    % Gibbs free energy of L + R binding                            **MED_ERROR**
    if true
        N_cT(3,num_cols + 4) = 1;   % RL
        N_cT(3,num_cols + 1) = -1;  % L
        N_cT(3,num_cols + 2) = -1;  % R
        G_0_bind = -45.1872; % kJ/mol
        R = 8.314;
        T = 310;
        K_bind = exp(G_0_bind/(R*T))*10^6;
    end
    
    K_C = [1 1 K_bind];
return 











