% PLB module
% all rxn directions align with kinetic Saucerman model (Km unchanged,
% aside from dimensional scaling)

function [k_kinetic, N_cT, K_C] = kinetic_parameters(M, include_all_reactions, num_cols)
    % Set the kinetic rate constants.
    % all reactions are reversible. no closed loops.
    % cAMP binds to R subunit one at a time
    
    % CONVERT TO fM 
    bigNum = 1e3;
    fastKineticConstant = bigNum;
    smallReverse = fastKineticConstant/(bigNum^2);
    u_to_f = 1e9;    
    
    cPP2A = 0.224; % uM (Bhalla Iyengar a on PMR)    but not sure what cell type this is

    % reaction IDs
%         PLBph1, PLBph2, 
%         PLBd1 (reverse), PLBd2 (reverse)
%         PKAi1 (reverse), PKAi2 (reverse)
%         PP2ai1, PP2ai2
%         inh

    k_pka_plb = 54; %: per_sec 54
    Km_pka_plb = 21; %: uM 21
    k_pp1_plb = 8.5; % per_sec 8.5
    Km_pp1_plb = 7; % uM 7
    k_pka_i1 = 60; %: per_sec 60
    Km_pka_i1 = 1; %: uM 1
    Vmax_pp2a_i1 = 14; %: uM_per_sec 14
    Km_pp2a_i1 = 1; %: uM 1
    Ki_inhib1 = 1e-3; %: uM 1e-3

    vKm = u_to_f*[Km_pka_plb, Km_pp1_plb, Km_pka_i1, Km_pp2a_i1];
    vkcat = [k_pka_plb, ...
        k_pp1_plb, ...
        k_pka_i1, ...
        Vmax_pp2a_i1/cPP2A];

    N = length(vKm); 
    k1m = zeros(1,N);
    k1p = zeros(1,N);
    k2m = zeros(1,N);
    k2p = zeros(1,N);

%     iforward = [1,4];
    for i=1:N % forward rxns
%         if find(i==iforward)
            k2p(i) = vkcat(i);
            k1m(i) = fastKineticConstant; % 1/s
            k1p(i) = (k1m(i) + k2p(i)) / (vKm(i)*u_to_f);
            k2m(i) = k1p(i)*k2p(i)/k1m(i); % detailed balance
%         else
%             k1m(i) = vkcat(i);
%             k2p(i) = fastKineticConstant; % 1/s
%             k2m(i) = (k2p(i) + k1m(i)) / (vKm(i)*u_to_f);
%             k1p(i) = k2m(i)*k1m(i)/k2p(i); % detailed balance
%         end
    end
    km_inh = fastKineticConstant;
    kp_inh = km_inh / (Ki_inhib1*u_to_f); % type 2
        
    if include_all_reactions
        k_kinetic = [
            k1p, k2p, kp_inh, k1m, k2m, km_inh
            ]';
    else
        irr = [1 2 3];
        k_kinetic = [
            k1p(irr), k2p(irr), k1m(irr), k2m(irr)
            ]';
    end

    % CONSTRAINTS
    N_cT = [];
%     N_cT = zeros(1,size(M,2)); 
% %     % Reaction i: [PKA:PKI] = [C][PKI] at SS   big error. Not isolated reaction
%     % repeat for type 1 and 2
%     if false
%         N_cT(1,num_cols + 2) = 1;
%         N_cT(1,num_cols + 6) = 1;
%         N_cT(1,num_cols + 8) = -1;
%     end
%     % Gibbs free energy of L + R binding                            **MED_ERROR**
%     if false
%         N_cT(3,num_cols + 4) = 1;   % ARC
%         N_cT(3,num_cols + 1) = -1;  % cAMP
%         N_cT(3,num_cols + 3) = -1;  % RC
%         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 = ones(1,size(N_cT,1));
return 











