# Kraeutler Reduced Hill Toy module

# Return kinetic parameters, constraints, and vector of volumes in each compartment (pL) 

# 8 reactions - 4 MM schemes of 2 reactions each
    # Re_AC1        
    # Re_AC2
    # Re_EC1
    # Re_EC2
    # Re_BD1
    # Re_BD2
    # Re_AND_E1
    # Re_AND_E2

# assume an OR gate is the same as a 1 node

# assume reactions follow classical Hill equation for (1) A*+B <-> complex
# then complex instantaneously breaks into (2) complex -> A* + B*
# the derivations for Km using Kraeutler parameters apply to reaction (1)

# apply stoichiometry of nS + E -> ES_n -> nP + E
# n = 1.4, hardcoded in stoich matrices

import numpy as np

def kinetic_parameters(M, dims, V):
    # Set the kinetic rate constants
    
    num_cols = dims['num_cols']
    num_rows = dims['num_rows']
    
    num_rx = 5
    K_Kraeut = 1.37695668165159
    E0 = [1]*num_rx # normalised enzyme concentrations
    Ymax = [1]*num_rx # A,B,C,D,E
    W = 1
    tau = 1
    EC50 = 0.5
    n = 1.4
    beta = (pow(EC50,n) - 1)/(2*pow(EC50,n) - 1)

    # units are dimensionless for molar quantities (all normalised)
    K1_m = K_Kraeut        #AC
    K2_m = K_Kraeut        #EC
    K3_m = K_Kraeut        #BD
    K4_m = K_Kraeut        #AND_CE

    # special parameters found for f_inh reaction that is fitted to a f_act reaction, which is in the form of MM
    K5_m = 0.5373114832002015        #AND_DE (inh fitted to f_act MM curve)
    beta_fit = 0.8682115872125091
    n_fit = -4.625919239202265

    Vmax = [0]*num_rx
    Vmax[0] = W*beta*Ymax[2]/tau
    Vmax[1] = W*beta*Ymax[2]/tau
    Vmax[2] = W*beta*Ymax[3]/tau
    Vmax[3] = W*beta*Ymax[4]/tau    # AND gate with same Ymax as E (for E1)

    Vmax[4] = W*beta_fit*Ymax[4]/tau    #W*beta*YmaxE/tau        # AND gate
    vkcat = []
    for iv, v in enumerate(Vmax):
        vkcat.append(v/E0[iv])  # kcat = Vmax/E0. NOT fastKineticConstant as it is NOT transient
    
    # initialise arrays
    vkap = np.zeros(num_rx)
    vkam = np.zeros(num_rx)
    vkbp = np.zeros(num_rx)
    vkbm = np.zeros(num_rx)
    vK_m = [K1_m,K2_m,K3_m,K4_m, K5_m]

    fastKineticConstant = 1e6 # 1/s
    smallReverse = 1e-9

    for i in range(num_rx):
        vkap[i] = fastKineticConstant
        vkam[i] = vkap[i]*vK_m[i] - vkcat[i]
        vkbp[i] = vkcat[i]
        vkbm[i] = (vkap[i]*vkbp[i])/vkam[i]    

    # add negative feedback reactions for every activated species [C, D, E], apart from stimuli [A, B]
    tau = 1
    k_negFB_p = [1/tau]*3
    k_negFB_m = [smallReverse]*3

    # Calculate bond graph constants from kinetic parameters
    # Note: units of kappa are fmol/s, units of K are fmol^-1
    k_kinetic = list(np.column_stack((vkap, vkbp)).flatten()) + k_negFB_p + list(np.column_stack((vkam, vkbm)).flatten()) + k_negFB_m

    # CONSTRAINTS
    N_cT = []
    K_C = []

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

    return (k_kinetic, N_cT, K_C, W)