# Gs protein module. From Dupont Sneyd: simple G protein model.
# Gs is associated with B1AR proteins

#     return (k_kinetic, N_cT, K_C, W) kinetic parameters, constraints, and vector of volumes in each
# compartment (pL) (1 if gating variable, or in element corresponding to
# kappa)

import numpy as np 

def kinetic_parameters(M, include_type2_reactions, dims, V):
    # 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

    num_cols = dims['num_cols']
    num_rows = dims['num_rows']

    bigNum = 1e3
    fastKineticConstant = bigNum
    smallReverse = fastKineticConstant/(pow(bigNum,2))
    c_G = 3.83 # uM
    
    kbindm = fastKineticConstant   # 1/s
    kbindp = kbindm / 0.285        # 1/uM.s
    kActp = 16/c_G                 # 1/uM.s
    kActm = smallReverse           # 1/uM.s            
    kHydp = 0.8                    # 1/s      
    kHydm = smallReverse           # 1/s
    kReassocp = 1.21e3             # 1/uM.s
    kReassocm = kReassocp/bigNum   # 1/s
    # phosphorylation by NDPK (nucleoside diphosphate kinase) - omitting MM enzyme
    kPiPhosp = fastKineticConstant
    kPiPhosm = smallReverse

    # CLOSED LOOP involving G - aGTP - aGDP - G
    # use detailed balance to find kReasocm
    if True:
        kReassocm = kActp*kHydp*kReassocp/(kActm*kHydm)
    
    
    k_kinetic = [
        kbindp, kActp, kHydp, kReassocp, kPiPhosp,
        kbindm, kActm, kHydm, kReassocm, kPiPhosm
        ]

    # CONSTRAINTS
    N_cT = np.zeros([2,len(M[1])]) 
    
    # Reaction BIND: equal portions of L + derivatives (L == LR)    **SMALL_ERROR**
    N_cT[0][num_cols] = -1
    N_cT[0][num_cols + 3] = 1
    # or LR.G = aGTP.betaGamma
    if False:                                                          
        N_cT[1][num_cols + 2] = 1
        N_cT[1][num_cols + 3] = 1
        N_cT[1][num_cols + 4] = -1
        N_cT[1][num_cols + 5] = -1
    
    # [a-GTP] + [a-GDP] = [beta.gamma]                              **SMALL_ERROR**
    if True: 
        N_cT[1][num_cols + 5] = 1   # beta_gamma
        N_cT[1][num_cols + 4] = -1  # a_GTP
        N_cT[1][num_cols + 6] = -1  # a_GDP
    
    # Gibbs free energy of L + R binding                            **MED_ERROR**
    if False:
        N_cT[0][num_cols + 3] = 1   # RL
        N_cT[0][num_cols] = -1  # L
        N_cT[0][num_cols + 1] = -1  # R
        G_0_bind = -45.1872 # kJ/mol
        R = 8.314
        T = 310
        K_bind = np.exp(G_0_bind/(R*T))*10^6
        K_C = [K_bind]
    
    
    K_C = [1, 1]

    # volume vector
    W = list(np.append([1] * num_cols, [V['V_myo']] * num_rows))

    return (k_kinetic, N_cT, K_C, W)