# PLC module, translated from Bhalla and Iyengar 1999

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

# Based on Pan 2018 cardiac AP

import numpy as np

def kinetic_parameters(M, include_type2_reactions, dims, V):
    # Set the kinetic rate constants

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

    fkc = 1e6
    smr = 1e-3

    # constants for reaction2 and reaction3 which are Michaelis-Menten
    Km = [5, 19.84166667] # [=] mM
    Vmax = [48, 10] # [=] 1/s
    E0 = [0.25, 1]  # [=] mM    #from cellML kinetic file
    kcat = [Vmax[i]/E0[i] for i in range(2)]

    k0p = 0.15  # R0 [=] 1/s
    k1p = 2.5  # R1 [=] 1/s
    k2ap = fkc  # R2a [=] 1/mM.s
    k2bp = kcat[0]  # R2b [=] 1/s
    k3ap = fkc  # R3a [=] 1/mM.s
    k3bp = kcat[1]  # R3b [=] 1/s
    k4p = 0.000005  # R4 [=] 1/mM.s
    k5p = 0.0000042  # R5 [=] 1/mM.s
    k6p = 0.00005  # R6 [=] 1/mM.s
    k7m = 0.000042  # R7 [=] 1/s   # reaction is reversed in kinetic version, so swap km and kp
    k8p = 0.0133  # R8 [=] 1/s

    k0m = smr          # R0 [=] 1/s
    k1m =  smr     # R1 [=] 1/s
    k2am =  (k2ap*Km[0]) - kcat[0]  # R2a [=] 1/s
    k2bm =  k2ap*kcat[0]/k2am  # R2b [=] 1/mM2.s
    k3am =  (k3ap*Km[1]) - kcat[1]  # R3a [=] 1/s
    k3bm =  k3ap*kcat[1]/k3am  # R3b [=] 1/mM2.s
    k4m = 1  # R4 [=] 1/s
    k5m = 1  # R5 [=] 1/s
    k6m = 1  # R6 [=] 1/s
    k7p = 1  # R7 [=] 1/mM.s
    k8m = smr  # R8 [=] 1/s

    # rate constants for reactions 0 - 8
    # note: closed loop present, but parameters are adjusted from kinetic model to meet detailed balance
    # closed loop for k7m and k8m. Different sign on k4 to make loop anticlockwise

    # k7p = k4m*k5p*k6p*k7m/(k4p*k5m*k6m) # error is ~zero when this commented out
    k8m = k4m*k5p*k6p*k8p/(k4p*k5m*k6m)

    # Calculate bond graph constants from kinetic parameters
    # Note: units of kappa are fmol/s, units of K are fmol^-1
    kf = [k0p,  # R0 [=] 1/s
          k1p,  # R1 [=] 1/s
          k2ap,  # R2a [=] 1/mM.s
          k2bp,  # R2b [=] 1/s
          k3ap,  # R3a [=] 1/mM.s
          k3bp,  # R3b [=] 1/s
          k4p,  # R4 [=] 1/mM.s
          k5p,  # R5 [=] 1/mM.s
          k6p,  # R6 [=] 1/mM.s
          k7p,  # R7 [=] 1/mM.s
          k8p,  # R8 [=] 1/s
          ]

    kr = [
        k0m,
        k1m,
        k2am,
        k2bm,
        k3am,
        k3bm,
        k4m,  # R4 [=] 1/s
        k5m,  # R5 [=] 1/s
        k6m,  # R6 [=] 1/s
        k7m,  # R7 [=] 1/s
        k8m  # R8 [=] 1/s
    ]

    k_kinetic = kf + kr

    # CONSTRAINTS
    N_cT = []
    K_C = []

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

    return (k_kinetic, N_cT, K_C, W)