## kinetic parameters for the LCC channel, adapted from Pan et al
## in the format of the other channels in the BG library

## hardcode the parameters from Pan's optimisation for ionic permeability and gate variables directly
import numpy as np 

def kinetic_parameters(M, include_type2_reactions, dims, V):

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

    # ## Import IV parameters
    x_LCC = 50000 / V['N_A'] * 1e15  # unit fmol
    x_LCC = V['numChannels']/V['N_A']*1e15 # unit fmol

    # L-type calcium channel
    # Load P_GHK (unit pL/s)
    # load([storage_dir 'LCC_P_GHK_Ca.mat'])
    P_LCC_Ca = 28.247076000000000 # P_GHK_Ca Unit pL/s
    # load([storage_dir 'LCC_P_GHK_K.mat'])
    P_LCC_K = 0.022204650000000  # P_GHK_K Unit pL/s
    # load([storage_dir 'LCC_d_parameters.mat'])
    params_d = [0.48676190823203597,
2.1404117552482584,
0.09802386860112802,
-2.1404117552482584
    ]

    # load([storage_dir 'LCC_f_parameters.mat'])
    params_f = [0.04688708435236662,	9.60767238499644,	0.9284601827359618,	8.534912099303016,	0.003047491353055938,	-2.4856670387239808,	0.045833689404789546,	-0.4710148075619314]

    # L-type calcium channel
    alpha_d0_bg = params_d[0]*1e3 # unit s^-1
    beta_d0_bg = params_d[2]*1e3 # unit s^-1

    alpha_f1_0_bg = params_f[0]*1e3 # unit s^-1
    beta_f1_0_bg = params_f[2]*1e3 # unit s^-1

    alpha_f2_0_bg = params_f[4]*1e3 # unit s^-1
    beta_f2_0_bg = params_f[6]*1e3 # unit s^-1

    K_f3_0 = beta_f1_0_bg*alpha_f2_0_bg/alpha_f1_0_bg/beta_f2_0_bg # dimensionless
    rate_f3 = 1e5 # unit s^-1

    KmCa = 0.6e-3 # Unit mM
    rate_fCa = 1e5 # unit s^-1

    kf_LCC = [P_LCC_Ca/x_LCC, # R_GHK_Ca1
        P_LCC_Ca/x_LCC, # R_GHK_Ca2
        P_LCC_K/x_LCC, # R_GHK_K1
        P_LCC_K/x_LCC, # R_GHK_K2
        alpha_d0_bg, # Rd000
        alpha_d0_bg, # Rd010
        alpha_d0_bg, # Rd020
        alpha_d0_bg, # Rd001
        alpha_d0_bg, # Rd011
        alpha_d0_bg, # Rd021
        alpha_f1_0_bg, # Rf1_000
        alpha_f1_0_bg, # Rf1_100
        alpha_f1_0_bg, # Rf1_001
        alpha_f1_0_bg, # Rf1_101
        alpha_f2_0_bg, # Rf2_000
        alpha_f2_0_bg, # Rf2_100
        alpha_f2_0_bg, # Rf2_001
        alpha_f2_0_bg, # Rf2_101
        K_f3_0*rate_f3, # Rf3_010
        K_f3_0*rate_f3, # Rf3_110
        K_f3_0*rate_f3, # Rf3_011
        K_f3_0*rate_f3, # Rf3_111
        rate_fCa, # RfCa000
        rate_fCa, # RfCa100
        rate_fCa, # RfCa010
        rate_fCa, # RfCa110
        rate_fCa, # RfCa020
        rate_fCa] # RfCa120

    kr_LCC = [P_LCC_Ca/x_LCC, # R_GHK_Ca1
        P_LCC_Ca/x_LCC, # R_GHK_Ca2
        P_LCC_K/x_LCC, # R_GHK_K1
        P_LCC_K/x_LCC, # R_GHK_K2
        beta_d0_bg, # Rd000
        beta_d0_bg, # Rd010
        beta_d0_bg, # Rd020
        beta_d0_bg, # Rd001
        beta_d0_bg, # Rd011
        beta_d0_bg, # Rd021
        beta_f1_0_bg, # Rf1_000
        beta_f1_0_bg, # Rf1_100
        beta_f1_0_bg, # Rf1_001
        beta_f1_0_bg, # Rf1_101
        beta_f2_0_bg, # Rf2_000
        beta_f2_0_bg, # Rf2_100
        beta_f2_0_bg, # Rf2_001
        beta_f2_0_bg, # Rf2_101
        rate_f3, # Rf3_010
        rate_f3, # Rf3_110
        rate_f3, # Rf3_011
        rate_f3, # Rf3_111
        rate_fCa/KmCa, # RfCa000
        rate_fCa/KmCa, # RfCa100
        rate_fCa/KmCa, # RfCa010
        rate_fCa/KmCa, # RfCa110
        rate_fCa/KmCa, # RfCa020
        rate_fCa/KmCa] # RfCa120

    k_kinetic = kf_LCC+kr_LCC

    # CONSTRAINTS
    N_cT = []
    K_C = []

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

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