# This Code is part of the publication
# Development of an experimental setup and a validated model for a highly integrated biocatalytic process
#
# Before use: adjust file path in line 18-21 to Experimental_data_SuppInfo/Miniplant from SuppInfo

import math

from numpy import linspace
from scipy.integrate import odeint
import numpy as np
import pandas as pd
from sklearn.metrics import mean_squared_error
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
from openpyxl import load_workbook

#--------------------------------------------Parameters
data_exp_aq = pd.read_excel('Exp 1.xlsx', sheet_name='Aqueous')
data_exp_org = pd.read_excel('Exp 1.xlsx', sheet_name='Organic')
data_exp_cond = pd.read_excel('Exp 1.xlsx', sheet_name='Conditions')
data_exp_error = pd.read_excel('Exp 1.xlsx', sheet_name='Error')

# run time
t_run = pd.to_numeric(data_exp_org.iloc[:, 0], errors='coerce').dropna().max()

# experimental error (relative)
error = data_exp_error.iloc[0, 13]/100

#reaction kinetics
# set 1: kinetic parameters, set 2: standard deviation high velocity, set3: standard deviation low velocity
vmax_sets = [
    {'Lip,for': 0.00036, 'ADH': 4.82, 'FDH': 5.17},
    {'Lip,for': 0.00038, 'ADH': 126.86, 'FDH': 6.72},
    {'Lip,for': 0.00034, 'ADH': 1e-5, 'FDH': 3.62},
]

K_sets = [
    {'eq': 25.89, 'mCOH': 4.52, 'mCAC': 10.37, 'mCCI': 0.156, 'iCAC': 1.14, 'iCOH': 4.064,
     'mCAL_ADH': 4.15, 'mNADH_ADH': 0.2, 'iCOH_ADH': 1.52, 'mNAD_FDH': 0.053, 'mFor_FDH': 112.4, 'iNADH_FDH': 0.027},

    {'eq': 1090, 'mCOH': 3.92, 'mCAC': 7.78, 'mCCI': 0.2, 'iCAC': 0.85, 'iCOH': 9.52,
     'mCAL_ADH': 1e-5, 'mNADH_ADH': 10e-5, 'iCOH_ADH': 18.4, 'mNAD_FDH': 0.0080, 'mFor_FDH': 23.05, 'iNADH_FDH': 0.061},

    {'eq': 1e-5, 'mCOH': 5.12, 'mCAC': 12.95, 'mCCI': 0.11, 'iCAC': 1.43, 'iCOH': 1e-20,
     'mCAL_ADH': 119.82, 'mNADH_ADH': 1.59, 'iCOH_ADH': 1e-5, 'mNAD_FDH': 0.098, 'mFor_FDH': 201.72, 'iNADH_FDH': 1e-5},
]

vLip_list = []
vADH_list = []
vFDH_list = []
c_mix_list=[]
result1, result2, result3 = [], [], []


#molare Mass
M={}
M['COH']=134.18
M['CAL']=132.16
M['CAC']=148.16

#Löslichkeit in Wasser [mM]
L_aq={}
L_aq['COH']= 1.8/M['COH']*1000
L_aq['CAL']= 1/M['CAL']*1000
L_aq['CAC']= 0.5/M['CAC']*1000



#patition coefficients
P={}
P['CAL'] = 8.4              # partition coefficient CAL
P['COH'] = 5.6              # partition coefficient COH


#--------------------------------------------Variables

#Volume and volume flow
Vol={}
Vol['aq']=data_exp_cond.iloc[14, 1]          # Volume aqueous phase (Tank+Modul) [l]
Vol['org']=data_exp_cond.iloc[15, 1]         # Volume organic phase [l]
Vol['Modul']=1.2        #Volume in Module [l]
Vol['aq tank']=Vol['aq']-Vol['Modul'] #Volume in aqueous tank [l]
Vol_flow={}
Vol_flow['aq']=data_exp_cond.iloc[10, 1]/1000     # Volume flow aqueous phase [l/min]
Vol_flow['org']=data_exp_cond.iloc[11, 1]/1000     # Volume flow organic phase [l/min]
Vol_flow['HPLC']=data_exp_cond.iloc[9, 1]/1000      # Volume flow of HPLC Pump into membrane module
Vol_flow['Back']=(Vol_flow['aq']-Vol_flow['HPLC'])  #overflow from Centrifuge directly to aqueous phase tank

#enzyme amount
m={}
m['Lipase'] = data_exp_cond.iloc[17, 1]
m['ADH']    = data_exp_cond.iloc[18, 1]
m['FDH']    = data_exp_cond.iloc[19, 1]

#initial concentrations in organic and aqueous phase [mM]
c_org_int={}
c_org_int['COH']=data_exp_cond.iloc[1, 2]
c_org_int['CAC']=data_exp_cond.iloc[2, 2]
c_org_int['CAL']=data_exp_cond.iloc[0, 2]
c_org_int['CCI']=data_exp_cond.iloc[3, 2]

c_aq_int={}
c_aq_int['COH']  =data_exp_cond.iloc[1, 1]
c_aq_int['CAC']  =data_exp_cond.iloc[2, 1]
c_aq_int['CAL']  =data_exp_cond.iloc[0, 1]
c_aq_int['CCI']  =data_exp_cond.iloc[3, 1]
c_aq_int['NADH'] =data_exp_cond.iloc[4, 1]
c_aq_int['NAD']  =data_exp_cond.iloc[5, 1]
c_aq_int['For']  =data_exp_cond.iloc[6, 1]

c_mem_int={}
c_mem_int['COH']  =data_exp_cond.iloc[1, 3]
c_mem_int['CAC']  =data_exp_cond.iloc[2, 3]
c_mem_int['CAL']  =data_exp_cond.iloc[0, 3]
c_mem_int['CCI']  =data_exp_cond.iloc[3, 3]
c_mem_int['NADH'] =data_exp_cond.iloc[4, 3]
c_mem_int['NAD']  =data_exp_cond.iloc[5, 3]
c_mem_int['For']  =data_exp_cond.iloc[6, 3]

#Partition coefficient Acid
def P_CAC(c,t):
    P_CAC  =  0.0002*c[1]**2-0.0081*c[1]+0.1511
    return (P_CAC)

#mixed concentrations in centrifuge
def mixcentri(c,t):
    P_Ac = P_CAC(c,t)
    c_orga_mixcentri_COH = (c[0] * Vol_flow['org'] + c[4] *Vol_flow['aq']) / (Vol_flow['org']+Vol_flow['aq']/P['COH'])
    c_orga_mixcentri_CAC = (c[1] * Vol_flow['org'] + c[5] *Vol_flow['aq']) / (Vol_flow['org']+Vol_flow['aq']/P_Ac)
    c_orga_mixcentri_CAL = (c[2] * Vol_flow['org'] + c[6] *Vol_flow['aq']) / (Vol_flow['org']+Vol_flow['aq']/P['CAL'])
    c_orga_mixcentri_CCI = c[3]
    c_aq_mixcentri_COH   = (c[0] * Vol_flow['org'] + c[4] *Vol_flow['aq']) / (Vol_flow['org']*P['COH']+Vol_flow['aq'])
    c_aq_mixcentri_CAC   = (c[1] * Vol_flow['org'] + c[5] *Vol_flow['aq']) / (Vol_flow['org']*P_Ac+Vol_flow['aq'])
    c_aq_mixcentri_CAL   = (c[2] * Vol_flow['org'] + c[6] *Vol_flow['aq']) / (Vol_flow['org']*P['CAL']+Vol_flow['aq'])
    c_aq_mixcentri_CCI   = 0
    c_mixcentri     = [c_orga_mixcentri_COH, c_orga_mixcentri_CAC,c_orga_mixcentri_CAL,c_orga_mixcentri_CCI, c_aq_mixcentri_COH,c_aq_mixcentri_CAC, c_aq_mixcentri_CAL, c_aq_mixcentri_CCI]
    return (c_mixcentri)

# reaction kinetics
def reaction_velocity_Lipase(c,t):
    c_mixcentri =mixcentri(c,t)
    v=vmax['Lip,for'] * m['Lipase'] * (c_mixcentri[0] * c_mixcentri[1] - c_mixcentri[3] / K['eq']) / (
                    K['iCAC'] * K['mCOH'] + K['mCOH'] * c_mixcentri[1] + K['mCAC'] * c_mixcentri[0] + c_mixcentri[1] * c_mixcentri[0]
                    + K['iCAC']*K['mCOH']*c_mixcentri[3]/K['mCCI']+K['iCAC']*K['mCOH']*c_mixcentri[3]/K['mCCI']/K['iCOH'])
    return (v)

def reaction_velocity_ADH(c,t):
    v=vmax['ADH']*c[15]*c[17]/(K['mCAL_ADH']*K['mNADH_ADH']*(1+c[12]/K['iCOH_ADH'])+K['mNADH_ADH']*c[15]+K['mCAL_ADH']*c[17]*(1+c[12]/K['iCOH_ADH'])+c[15]*c[17])*m['ADH']
    return (v)

def reaction_velocity_FDH(c,t):
    v=vmax['FDH']*c[9]*c[10]/(K['mNAD_FDH']*K['mFor_FDH']*(1+c[8]/K['iNADH_FDH'])+K['mFor_FDH']*c[9]+K['mNAD_FDH']*c[10]*(1+c[8]/K['iNADH_FDH'])+c[9]*c[10])*m['FDH']
    return (v)



def molar_balance(c, t, vmax,K):
    global vLip_list, vADH_list, vFDH_list

    for i in range(0, 19):
        if c[i] < 0:
            c[i] = 0


    dNprecCAL=0
    dNprecCALmem=0
    dNprecCOH=0
    dNprecCOHmem=0
    vLip=reaction_velocity_Lipase(c,t)
    vADH=reaction_velocity_ADH(c,t)
    vFDH=reaction_velocity_FDH(c,t)

    vLip_list.append(vLip)
    vADH_list.append(vADH)
    vFDH_list.append(vFDH)

    from math import e
    NADHoxi = c[8]*(1-0.66*e**(-0.00179*t))
    NADHoxiMem = c[17] * (1 - 0.66 * e ** (-0.00179 * t))

    c_mix=mixcentri(c,t)
    c_mix_list.append(c_mix)

    dcorgaCOHdt = -vLip/Vol['org'] + (c_mix[0]-c[0])*Vol_flow['org']/Vol['org']
    dcorgaCACdt = -vLip/Vol['org'] + (c_mix[1]-c[1])*Vol_flow['org']/Vol['org']
    dcorgaCALdt = (c_mix[2]-c[2])*Vol_flow['org']/Vol['org']
    dcorgaCCIdt = vLip/Vol['org']

    dcaqCOHdt = (c[13]-c[4])*Vol_flow['HPLC']/Vol['Modul'] + (c_mix[4]-c[4])*Vol_flow['Back']/(Vol['aq tank'])-dNprecCOH
    dcaqCACdt = (c_mix[5]-c[5])*Vol_flow['aq']/Vol['aq']
    dcaqCALdt = (c[15]-c[6])*Vol_flow['HPLC']/Vol['Modul']+(c_mix[6]-c[6])*Vol_flow['Back']/(Vol['aq tank'])-dNprecCAL
    dcaqCCIdt = (c_mix[7]-c[7])*Vol_flow['aq']/Vol['aq']
    dcaqNADHdt= (c[17]-c[8])*Vol_flow['HPLC']/Vol['Modul'] -NADHoxi
    dcaqNADdt = (c[18]-c[9])*Vol_flow['HPLC']/Vol['Modul'] +NADHoxi
    dcaqFordt = (c[19]-c[10])*Vol_flow['HPLC']/Vol['Modul']

    dcmemCOHdt = vADH / Vol['Modul'] + (c_mix[4] - c[13]) * Vol_flow['HPLC'] / Vol['Modul']-dNprecCOHmem
    dcmemCACdt = (c_mix[5] - c[14]) * Vol_flow['HPLC'] / Vol['Modul']
    dcmemCALdt = -vADH / Vol['aq'] + (c_mix[6] - c[15]) * Vol_flow['HPLC'] / Vol['Modul']-dNprecCALmem
    dcmemCCIdt = (c_mix[7] - c[16]) * Vol_flow['aq'] / Vol['aq']
    dcmemNADHdt = vFDH / Vol['Modul'] - vADH / Vol['Modul'] - NADHoxiMem + (c[8] - c[17]) * Vol_flow['HPLC'] / Vol['Modul']
    dcmemNADdt = vADH / Vol['Modul'] - vFDH / Vol['aq'] + NADHoxiMem + (c[9] - c[18]) * Vol_flow['HPLC'] / Vol['Modul']
    dcmemFordt = -vFDH / Vol['Modul'] + (c[10] - c[19]) * Vol_flow['HPLC'] / Vol['Modul']

    dcdt    = [dcorgaCOHdt, dcorgaCACdt, dcorgaCALdt, dcorgaCCIdt, dcaqCOHdt, dcaqCACdt, dcaqCALdt, dcaqCCIdt, dcaqNADHdt, dcaqNADdt, dcaqFordt, dNprecCAL, dNprecCOH, dcmemCOHdt, dcmemCACdt, dcmemCALdt, dcmemCCIdt, dcmemNADHdt, dcmemNADdt, dcmemFordt, dNprecCALmem, dNprecCOHmem]

    return (dcdt)


#---------------------simulation

t = linspace(0,t_run, int(t_run+1))
c_0 = [c_org_int['COH'], c_org_int['CAC'], c_org_int['CAL'], c_org_int['CCI'], c_aq_int['COH'] , c_aq_int['CAC'] , c_aq_int['CAL'] , c_aq_int['CCI'], c_aq_int['NADH'], c_aq_int['NAD'], c_aq_int['For'], Prec_int['CAL'],Prec_int['COH'], c_mem_int['COH'] , c_mem_int['CAC'] , c_mem_int['CAL'] , c_mem_int['CCI'], c_mem_int['NADH'], c_mem_int['NAD'], c_mem_int['For'], Prec_mem_int['CAL'],Prec_mem_int['COH']]
# c is concentration vector for organic (c[0]-c[3]) and aqueous tank (c[4]-c[10]) and percipitation (c[11] and c[12]) and membrane modul(c[13]-c[19])
t = linspace(0,int(round(t_run)), int(round(t_run))+1)

def run_simulation(vmax, K):
    return odeint(molar_balance, c_0, t, args=(vmax, K))


vLip_array = np.array(vLip_list)
vADH_array = np.array(vADH_list)
vFDH_array = np.array(vFDH_list)

# concentration indices
#0 COH org
#1 CAC org
#2 CAL org
#3 CCI org
#4 COH aq
#5 CAC aq
#6 CAL aq
#7 CCI aq
#8 NADH aq
#9 NAD aq
#10 FOR aq
#11 CAL perc
#12 COH perc
#13 COH mem
#14 CAC mem
#15 CAL mem
#16 CCI mem
#17 NADH mem
#18 NAD mem
#19 FOR mem

#COH = Vol['aq tank']*c[:,4]+Vol['Modul']*c[:,13]+Vol['org']*c[:,0]+Vol['org']*c[:,3]

for i, (vmax, K) in enumerate(zip(vmax_sets, K_sets)):
    result = run_simulation(vmax, K)

    if i == 0:
        result1 = result
    elif i == 1:
        result2 = result
    elif i == 2:
        result3 = result

    print(f"Simulation {i+1} abgeschlossen.")

t_exp_aq = pd.to_numeric(data_exp_aq.iloc[:, 0], errors='coerce').dropna()
c_aq_exp_ald = pd.to_numeric(data_exp_aq.iloc[:, 1], errors='coerce').dropna()
c_aq_exp_alc = pd.to_numeric(data_exp_aq.iloc[:, 2], errors='coerce').dropna()

t_exp_org = pd.to_numeric(data_exp_org.iloc[:, 0], errors='coerce').dropna()
c_org_exp_alc = pd.to_numeric(data_exp_org.iloc[:, 2], errors='coerce').dropna()
c_org_exp_ac = pd.to_numeric(data_exp_org.iloc[:, 3], errors='coerce').dropna()
c_org_exp_est = pd.to_numeric(data_exp_org.iloc[:, 4], errors='coerce').dropna()

#---------------------plot


# Compute experimental error
error_ald = error * c_aq_exp_ald
error_alc = error * c_aq_exp_alc

t_h=t/60
t_exp_org_h=t_exp_org/60
t_exp_aq_h=t_exp_aq/60
f1 =plt.figure(1)
plt.plot(t_h,result1[:,6], 'black', label = 'aldehyde', linestyle = 'dashed')
plt.fill_between(t_h, result2[:,6], result3[:,6], color='grey', alpha=0.3)
plt.errorbar(t_exp_aq_h, c_aq_exp_ald, yerr=error_ald, fmt='.', color='black', label='aldehyde exp', linestyle='none')
plt.errorbar(t_exp_aq_h, c_aq_exp_alc, yerr=error_alc, fmt='.', color='limegreen', label='alcohol exp', linestyle='none')
plt.plot(t_h,result1[:,4], 'limegreen', label = 'alcohol', linestyle = 'dashed')
plt.fill_between(t_h, result2[:,4], result3[:,4], color='limegreen', alpha=0.3)
plt.xlabel('time / h', fontsize=14)
plt.ylabel('concentration / mM', fontsize=14)
plt.legend(loc='center right',frameon=True,
    facecolor='white',
    edgecolor='black',
    framealpha=1.0,)
plt.xlim(0,11)
plt.ylim(0,1.5)



error_alc_org = error * c_org_exp_alc
error_ac_org = error * c_org_exp_ac
error_est_org = error * c_org_exp_est

fig, ax1 = plt.subplots()

ax2 = ax1.twinx()
ax1.plot(t_h,result1[:,0], 'limegreen', label = 'alcohol', linestyle = 'dashed')
ax1.fill_between(t_h, result2[:,0], result3[:,0], color='limegreen', alpha=0.3)
ax1.plot(t_h,result1[:,1], 'yellow', label = 'acid', linestyle = 'dashed')
ax1.fill_between(t_h, result2[:,1], result3[:,1], color='yellow', alpha=0.3)
#ax1.plot(t,result1[:,2], 'blue', label = 'aldehyde', linestyle = 'dashed')
#plt.fill_between(t, result2[:,2], result3[:,2], color='blue', alpha=0.3, label='Spanne Result 1–3')
ax1.errorbar(t_exp_org_h,c_org_exp_alc,yerr=error_alc_org, fmt='.', color='limegreen', label = 'alcohol exp', linestyle = 'none')
ax1.errorbar(t_exp_org_h,c_org_exp_ac,yerr=error_ac_org, fmt='.', color= 'yellow', label = 'acid exp', linestyle = 'none')


ax2.plot(t_h,result1[:,3], 'red', label = 'ester', linestyle = 'dashed')
#ax2.plot(t,result2[:,3], 'red', label = 'ester', linestyle = 'dashed')
#ax2.plot(t,result3[:,3], 'red', label = 'ester', linestyle = 'dashed')
plt.fill_between(t_h, result2[:,3], result3[:,3], color='red', alpha=0.3)
ax2.errorbar(t_exp_org_h,c_org_exp_est, yerr=error_est_org, fmt='.', color='red', label = 'ester exp', linestyle = 'none')
#f2.plot(t,c[:,4], 'darkgray', label = 'alcohol')
ax1.set_xlabel('time / h', fontsize=14)
ax1.set_ylabel('concentration alcohol / mM', fontsize=14)
ax1.set_ylim([0, 7])
ax1.set_xlim([0, 11])
ax2.set_ylabel('concentration ester / mM', fontsize=14)
ax2.set_ylim([0, 0.7])
handles1, labels1 = ax1.get_legend_handles_labels()
handles2, labels2 = ax2.get_legend_handles_labels()
handles_ordered = [handles1[0], handles1[1], handles2[0]] + handles1[2:] + handles2[1:]
labels_ordered  = [labels1[0],  labels1[1],  labels2[0]]  + labels1[2:]  + labels2[1:]
combined_legend = ax2.legend(
    handles_ordered, labels_ordered,
    loc='lower right',
    frameon=True,
    facecolor='white',
    edgecolor='black',
    framealpha=1.0)

combined_legend.set_zorder(10)

plt.show()

# Determination of model values
model_interp_func_aq_alc = interp1d(t, result1[:, 4], kind='linear', bounds_error=False, fill_value='extrapolate')
model_interp_func_aq_ald = interp1d(t, result1[:, 6], kind='linear', bounds_error=False, fill_value='extrapolate')
model_interp_func_org_alc = interp1d(t, result1[:, 0], kind='linear', bounds_error=False, fill_value='extrapolate')
model_interp_func_org_ac = interp1d(t, result1[:, 1], kind='linear', bounds_error=False, fill_value='extrapolate')
model_interp_func_org_est = interp1d(t, result1[:, 3], kind='linear', bounds_error=False, fill_value='extrapolate')

c_model_interp_aq_alc = model_interp_func_aq_alc(t_exp_aq)
c_model_interp_aq_ald = model_interp_func_aq_ald(t_exp_aq)
c_model_interp_org_alc = model_interp_func_org_alc(t_exp_org)
c_model_interp_org_ac = model_interp_func_org_ac(t_exp_org)
c_model_interp_org_est = model_interp_func_org_est(t_exp_org)

# RMSE calculation
rmse_aq_alc = np.sqrt(mean_squared_error(c_aq_exp_alc, c_model_interp_aq_alc))
rmse_aq_ald = np.sqrt(mean_squared_error(c_aq_exp_ald, c_model_interp_aq_ald))
rmse_org_alc = np.sqrt(mean_squared_error(c_org_exp_alc, c_model_interp_org_alc))
rmse_org_ac = np.sqrt(mean_squared_error(c_org_exp_ac, c_model_interp_org_ac))
rmse_org_est = np.sqrt(mean_squared_error(c_org_exp_est, c_model_interp_org_est))

# normalized RSME
per_aq_alc = 100*rmse_aq_alc/np.mean(c_aq_exp_alc)
per_aq_ald = 100*rmse_aq_ald/np.mean(c_aq_exp_ald)
per_org_alc = 100*rmse_org_alc/np.mean(c_org_exp_alc)
per_org_ac = 100*rmse_org_ac/np.mean(c_org_exp_ac)
per_org_est = 100*rmse_org_est/np.mean(c_org_exp_est)

print(f"RMSE for alcohol in aqueous phase: {rmse_aq_alc:.4f} mM, normalized RSME : {per_aq_alc:.4f} %")
print(f"RMSE for aldehyde in aqueous phase: {rmse_aq_ald:.4f} mM, normalized RSME : {per_aq_ald:.4f} %")
print(f"RMSE for alcohol in organic phase: {rmse_org_alc:.4f} mM, normalized RSME : {per_org_alc:.4f} %")
print(f"RMSE for acid in organic phase: {rmse_org_ac:.4f} mM, normalized RSME : {per_org_ac:.4f} %")
print(f"RMSE for ester in organic phase: {rmse_org_est:.4f} mM, normalized RSME : {per_org_est:.4f} %")


