% B1AR module

function [k_kinetic, N_cT] = 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 
    % oncentration of BARK = 0.6uM (crude approx from litsearch, for GRK)
    bigNum = 1e6;
    fastKineticConstant = bigNum;
    smallReverse = fastKineticConstant/(bigNum^2);
    u_to_f = 1e9;
    cBARK = 0.6*u_to_f;            %fM
    cPKACI = 0.0587898478276673*u_to_f;        %fM
    
    % creation of LRG - Type 2 reactions
    % [K_C, K_L, K_R, K_L]
    vK_mSig = [33*u_to_f, 0.285*u_to_f, 1/(0.062*u_to_f), 1/(0.285*u_to_f)]; % fM
    % for detailed balance, let Km3 = 1/(Km1*Km2*Km4). Source paper values are
    % thermo inconsistent
    vK_mSig(2) = 1/(vK_mSig(1)*vK_mSig(3)*vK_mSig(4));
    
    vkSigm = zeros(1,4);
    vkSigp = zeros(1,4);
    for i= 1:4
        vkSigm(i) = fastKineticConstant;    %1/s for Sig1, Sig2. 1/fM.s for Sig 3, Sig4
        vkSigp(i) = vkSigm(i)/vK_mSig(i);   %1/s
    end
    % detailed balance for last reverse ksigm(4)
    vkSigm(4) = vkSigp(1)*vkSigp(2)*vkSigp(3)*vkSigp(4)/vkSigm(1)*vkSigm(2)*vkSigm(3);
%     old_kSigm4 = vkSigp(1)*vkSigm(2)*vkSigm(3)*vkSigp(4)/vkSigm(1)*vkSigp(2)*vkSigp(3);


    kB11p = 1.1e-3 / cBARK;      % 1/fM.s
    kB11m = 2.2e-3;             % 1/s
    kB12p = fastKineticConstant;%1/s
    kB12m = smallReverse;       % 1/fM.s
    kB21p = kB11p;              % 1/fM.s
    kB21m = kB11m;              % 1/s
    kB22p = fastKineticConstant;%1/s
    kB22m = smallReverse;       % 1/fM.s
    kP1p = 2.2e-3/ cPKACI;          % 1/fM.s (kp_kam)
    kP1m = smallReverse;        % 1/s
    kP2p = fastKineticConstant; % 1/s
    kP2m = 3.6e-3/u_to_f;   % 1/fM.s (kp_kap)
    kAct1p = 16;            % 1/s
    kAct1m = smallReverse;          % 1/fM.s            
    kAct2m = 16;            % 1/s
    kAct2p = 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
    
    % detailed balance for kP2p
    kP2p = kP2m*kP1m/kP1p;
    
    k_kinetic = [
        kB11p, kB12p, kB21p, kB22p, kP1p, kP2p, kAct1p, kAct2p, kHydp, kReassocp, vkSigp...
        kB11m, kB12m, kB21m, kB22m, kP1m, kP2m, kAct1m, kAct2m, kHydm, kReassocm, vkSigm
        ]';

    % CONSTRAINTS
    N_cT = zeros(2,size(M,2)); 
    % substrate LR/LRG is in eqlm with product B1ARd
    N_cT(1,num_cols + 1) = 1;
    N_cT(1,num_cols + 2) = -1;
    N_cT(2,num_cols + 5) = 1;
    N_cT(2,num_cols + 2) = -1;
    % substrate B1ARp is in eqlm with product B1ARtot
    N_cT(3,num_cols + 7) = 1;
    N_cT(3,num_cols + 8) = -1;


return 