# Kinetic rate constants for composite model of FCU-AC 
# they must be computed together because of the presence of new potential closed loops

# order: {'cAMP','LRGbinding_B1AR','LRGbinding_M2','GsProtein','GiProtein'}

import numpy as np

def grand_kinetic_parameters(M, include_type2_reactions, dims, V, noLRG):

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

    # shared large/small constants
    bigNum = 1e5
    fastKineticConstant = bigNum
    smallReverse = 1e0
    
# M2 composite: LRG_M2 + GiProtein
    
    K_M2sig1 = 30    # uM  Kc
    K_M2sig2 = 0.16  # uM  Kh
    K_M2sig4 = 11    # uM  Kl
    k_M2sig1p = fastKineticConstant
    k_M2sig1m = k_M2sig1p*K_M2sig1
    k_M2sig2p = fastKineticConstant
    k_M2sig2m = k_M2sig2p*K_M2sig2
    k_M2sig4p = fastKineticConstant
    k_M2sig4m = k_M2sig4p*K_M2sig4
    k_M2sig3p = fastKineticConstant
    k_M2sig3m = k_M2sig1m*k_M2sig2m*k_M2sig3p*k_M2sig4p/(k_M2sig1p*k_M2sig2p*k_M2sig4m) # interim?

    k_Doff1p = fastKineticConstant
    k_Doff1m = smallReverse
    k_Ton1p = fastKineticConstant
    k_Ton1m = smallReverse
    k_Doff2p = fastKineticConstant
    k_Doff2m = smallReverse
    k_Ton2p = fastKineticConstant
    k_Ton2m = smallReverse
    k_M2Act1p = 2.5                   # 1/s
    k_M2Act1m = smallReverse          # 1/s            
    k_M2Act2p = 0.05                  # 1/s
    k_M2Act2m = smallReverse          # 1/s            
    k_M2Hydp = 0.8                    # 1/s      
    k_M2Hydm = smallReverse           # 1/s
    k_M2Reassocp = 1.21e3             # 1/uM.s
    k_M2Reassocm = k_M2Reassocp/fastKineticConstant   # 1/s
    # phosphorylation of GDP by NDPK (nucleoside diphosphate kinase) - omitting MM enzyme

    # Loop1: find kact1m using detailed balance
    k_M2Act1m = (k_M2Act1p*k_M2sig1p)/k_M2sig1m
    # Loop2: find kact1m using detailed balance
    k_M2Act2m = (k_M2Act2p*k_M2sig3p)/k_M2sig3m

    # Big Loop: find k_M2sig3p or kact2m using detailed balance
    # *** unnecessary ***
    k_M2sig3m = k_M2Act1p*k_M2sig2m*k_M2sig3p*k_M2sig4p/(k_M2Act1m*k_M2sig2p*k_M2sig4m)
    
    # CLOSED LOOP involving G - aGTP - aGDP - G
    # use detailed balance to find kReasocm with either Act (as they have
    # same equilibrium constant
    if True:
        k_M2Reassocm = k_M2Act1p*k_M2Hydp*k_M2Reassocp/(k_M2Act1m*k_M2Hydm)
    
    
# TOTAL    
    #  order: {'cAMP','LRGbinding_B1AR','LRGbinding_M2','GsProtein','GiProtein'}
# forwaard first, then reverse!
    k_kinetic = [k_M2sig1p, k_M2sig2p, k_M2sig3p, k_M2sig4p,# LRG_M2
        k_Doff1p, k_Ton1p, k_M2Act1p, k_Doff2p, k_Ton2p, k_M2Act2p, k_M2Hydp, k_M2Reassocp, #GiProtein
        k_M2sig1m, k_M2sig2m, k_M2sig3m, k_M2sig4m,
        k_Doff1m, k_Ton1m, k_M2Act1m, k_Doff2m, k_Ton2m, k_M2Act2m, k_M2Hydm, k_M2Reassocm
        ]


    # CONSTRAINTS
    N_cT = []
    if False:
        N_cT = np.zeros([1,size(M,2)]) 
        # substrate B1ARp is in eqlm with product B1ARtot
        N_cT[0][num_cols] = 1
        N_cT[0][num_cols + 1] = -1
    
    K_C = np.ones([1,len(N_cT)])

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

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