# Merged GPCR module, which is a combination of GsProtein and LRGbinding_B1AR.
# following Saucerman and Iancu: act1 and act2 with LR
# and LRG as substrates (G is already bound, so there is only one substrate
# to each act reaction)

# 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)

# 14Oct21: adding phosphate binding reaction to RG and LRG as separate
# reactions. Only RG_Pi and LRG_Pi proceed to Act reactions. 

# 23nov21 removed RG + L -> LRG (sig2)
# 24nov21: merged Gs and GPCR binding
# 24nov21: adding receptor internalisation after aGTP has dissociated from the RG/LRG complex

#19nov24: merge inactive R to active. Remove Rswitch reaction.
#9 dec remove internalisation aspect (but retain Rtag)
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 = 1e2
    fastKineticConstant = bigNum
    smallReverse = 1/bigNum
    medReverse = np.sqrt(smallReverse)

    k_Rswitchp = fastKineticConstant
    k_Rswitchm = k_Rswitchp*0.01  # assume that only 1% of receptors are constitutionally active
    k_LRswitchp = fastKineticConstant
    k_LRswitchm = smallReverse

    KRc = 30    # uM  Kc
    KRh = 0.16  # uM  Kh
    KRL = 11    # uM  Kl
    # below modification: get in its of mM
    if False:
        KRc *= 1e-3
        KRh *= 1e-3
        KRL *= 1e-3
    KRr = KRL/(KRc*KRh)
    RCLpow = 1.75
    kRcp = pow(fastKineticConstant,RCLpow)
    kRcm = kRcp*KRc
    kRLp = pow(fastKineticConstant,RCLpow)
    kRLm = kRLp*KRL
    kRrp = pow(fastKineticConstant,RCLpow)
    # find kRLm using detailed balance
    # kRrm = kRcm*kRhm*kRrp*kRLp/(kRcp*kRhp*kRLm)
    kRrm = kRrp * KRr
    # kRrm = smallReverse

    kAct1p = 0.05                 # 1/s
    kAct1m = smallReverse #medReverse           # 1/s
    kAct2p = 2.5                 # 1/s
    kAct2m = smallReverse #medReverse           # 1/s
    kHydp = 0.8                    # 1/s      
    kHydm = smallReverse #medReverse           # 1/s
    kReassocp = 1.21e0 # 1.21e3             # 1/uM.s
    kReassocm = smallReverse #medReverse   # 1/s

    # modify
    if False:
        # kAct1p *= 100
        kAct2p *= 1e4
        # kHydp *= 100

    # RECEPTOR INTERNALISATION
    k_interRp = fastKineticConstant
    k_interRm = smallReverse
    k_interLRp = fastKineticConstant
    k_interLRm = smallReverse

    # ensure that the closed loop formed by Act1 & Act2 obey detailed balance
    # CLOSED LOOP involving G - aGTP - aGDP - G
    # use detailed balance to find kReasocm with either Act (as they have
    # same equilibrium constant
    if False:
        kAct2m = kAct1m * kAct2p / kAct1p
        kReassocm = kRcp*kAct1p*kHydp*kReassocp/(kRcm*kAct1m*kHydm)

    k_kinetic = [
        # kRcp, kRrp, kRLp, kAct1p, kAct2p, kHydp, kReassocp,
        # kRcm, kRrm, kRLm, kAct1m, kAct2m, kHydm, kReassocm
        k_Rswitchp,k_LRswitchp,kRcp, kRrp, kRLp, kAct1p, kAct2p, kHydp, kReassocp, k_interRp, k_interLRp,
        k_Rswitchm,k_LRswitchm,kRcm, kRrm, kRLm, kAct1m, kAct2m, kHydm, kReassocm, k_interRm, k_interLRm
        ]

    # CONSTRAINTS
    N_cT = []
    K_C = []


    # [a-GTP] + [a-GDP] = [beta.gamma]                              **SMALL_ERROR**
    if False:
        N_cT = np.zeros(len(M[0]))
        N_cT[num_cols + 6] = 1   # beta_gamma
        N_cT[num_cols + 5] = -1  # a_GTP
        N_cT[num_cols + 7] = -1  # a_GDP
        K_C = [1]

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

    return (k_kinetic, [N_cT], K_C, W)