#%%
#Libraries
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({'font.size': 13})
from matplotlib.widgets import Slider, Button

#%%
# Explicit equations
mu1max =  .9  
mu2max =  .6  
K1 =  .1  
K2 =  .5  
alpha1 =  .0001  
alpha2 =.0001  
beta1 =  .05  
beta2 =  .05  
Ycs1 =  .4  
Ycs2 =  .4
def ODEfun(Yfuncvec, t, mu1max, mu2max, K1, 
           K2, alpha1, alpha2, beta1, beta2,
           Ycs1, Ycs2): 
    Cc= Yfuncvec[0]
    Cs1= Yfuncvec[1]
    Cs2= Yfuncvec[2]
    e1=Yfuncvec[3]
    e2=Yfuncvec[4]

    e1max =  alpha1/( mu1max+ beta1)
    e2max =  alpha2/( mu2max+ beta2)
    E1 =  e1/e1max
    E2 =  e2/e2max
    mu1 =  mu1max*E1*Cs1/(K1+Cs1)
    mu2 =  mu2max*E2*Cs2/(K2+Cs2)
    u1 =  mu1/(mu1+mu2)
    u2 =  mu2/(mu1+mu2)
    mumax = np.where(mu1 > mu2, mu1, mu2) 
    v1 =  mu1/ mumax
    v2 =  mu2/mumax
    # Differential equations
    dCcdt =(mu1*v1+mu2*v2)*Cc 
    dCs1dt = -mu1*v1*Cc/Ycs1  
    dCs2dt = -mu2*v2*Cc/Ycs2  
    de1dt = alpha1*Cs1/(K1+Cs1)*u1 - beta1*e1 - (mu1*v1+mu2*v2)*e1  
    de2dt = alpha2*Cs2/(K2+Cs2)*u2 - beta2*e2 - (mu1*v1+mu2*v2)*e2  
    return np.array([dCcdt,dCs1dt,dCs2dt,de1dt,de2dt])
 

tspan = np.linspace(0, 10, 1000) # Range for the independent variable
y0 = np.array([0.1,4,20,8.3e-5,1e-6]) # Initial values for the dependent variables

#%%
fig, (ax1,ax2) = plt.subplots(2,1)
fig.suptitle("""LEP-PRS-9-6 Two Substrate""", x = 0.2, y=0.98, fontweight='bold')
plt.subplots_adjust(left  = 0.5)
fig.subplots_adjust(wspace=0.25,hspace=0.3)
sol = odeint(ODEfun, y0, tspan, (mu1max, mu2max, K1, 
                                 K2, alpha1, alpha2, beta1, beta2,
                                 Ycs1, Ycs2))
Cc = sol[:,0]
Cs1= sol[:,1]
Cs2= sol[:,2]
e1 = sol[:,3]
e2 = sol[:,4]

p1, p2, p3 = ax1.plot(tspan, Cc, tspan, Cs1,tspan, Cs2)
ax1.legend(['$C_C(g/L)$','$C_{S1}(g/L)$','$C_{S2}(g/L)$'], loc='best')
ax1.set_xlabel('time (hr)', fontsize='medium')
ax1.set_ylabel('Concentration (M)', fontsize='medium')
#ax1.set_ylim(0,20)
ax1.set_xlim(0,4)
ax1.grid()

p4, p5 = ax2.plot(tspan, e1,tspan, e2)
ax2.legend(["$e_1(g_{enzyme}/g_{cell mass})$", "$e_2(g_{enzyme}/g_{cell mass})$"], loc='best')
ax2.set_xlabel('time (hr)', fontsize='medium')
ax2.set_ylabel('Concentration (M)', fontsize='medium')
ax2.set_ylim(0,0.0012)
ax2.set_xlim(0,4)
ax2.grid()


ax1.text(-4.7, -28,'Differential Equations'
         '\n'
         r'$\dfrac{dC_C}{dt} = (\mu_1 v_1 + \mu_2 v_2) C_C$'
                  '\n \n'
         r'$\dfrac{dC_{S1}}{dt} = \dfrac{-\mu_1 v_1 C_C}{Y_{c/s1}}$'
                  '\n \n'
         r'$\dfrac{dC_{S2}}{dt} = \dfrac{-\mu_2 v_2 C_C}{Y_{c/s2}}$'
                  '\n \n'                  
         r'$\dfrac{de_1}{dt} = \dfrac{\alpha_1 C_{S1}u_1}{K_1+C_{S1}} - \beta_1 e_1 - (\mu_1 v_1 + \mu_2 v_2)e_1$'
                  '\n \n'     
         r'$\dfrac{de_2}{dt} = \dfrac{\alpha_2 C_{S2}u_2}{K_2+C_{S2}} - \beta_2 e_2 - (\mu_1 v_1 + \mu_2 v_2)e_2$'
                  '\n \n'                  
         'Explicit Equations'
                  '\n'
         r'$e_{1max} =  \dfrac{\alpha_1}{mu_{1max}+ \beta_1}$'
         '\n'      
         r'$e_{2max} =  \dfrac{\alpha_2}{mu_{2max}+ \beta_2}$'
         '\n'          
         r'$E_1 =  \dfrac{e_1}{e_{1max}}$'
         '\n'
         r'$E_2 =  \dfrac{e_2}{e_{2max}}$'
         '\n'         
         r'$\mu_1 =  \dfrac{\mu_{1max}E_1C_{S1}}{K_1+C_{S1}}$'
         '\n'         
         r'$\mu_2 =  \dfrac{\mu_{2max} E_2 C_{S2}}{K_2+C_{S2}}$'
         '\n'          
         r'$u_1 =  \dfrac{\mu_1}{\mu_1 + \mu_2}$'
         '\n'          
         r'$u_2 =  \dfrac{\mu_2}{\mu_1 + \mu_2}$'
         '\n'          
         r'$If (\mu_1 > \mu_2) then (\mu_{max} = \mu_1) else (\mu_{max} = \mu_2)$'    
         '\n'          
         r'$v_1 =  \dfrac{\mu_1}{\mu_{max}}$'
         '\n'          
         r'$v_2 =  \dfrac{\mu_2}{\mu_{max}}$'
         
         , ha='left', wrap = True, fontsize=12,
        bbox=dict(facecolor='none', edgecolor='black', pad=15), fontweight='bold')

#%%
# Slider
axcolor = 'black'
ax_Ycs1 = plt.axes([0.28, 0.8, 0.12, 0.02], facecolor=axcolor)
ax_alpha1 = plt.axes([0.28, 0.75, 0.12, 0.02], facecolor=axcolor)
ax_beta1 = plt.axes([0.28, 0.7, 0.12, 0.02], facecolor=axcolor)
ax_K1 = plt.axes([0.28, 0.65, 0.12, 0.02], facecolor=axcolor)
ax_mumax1 = plt.axes([0.28, 0.6, 0.12, 0.02], facecolor=axcolor)
ax_Ycs2 = plt.axes([0.28, 0.55, 0.12, 0.02], facecolor=axcolor)
ax_alpha2 = plt.axes([0.28, 0.5, 0.12, 0.02], facecolor=axcolor)
ax_beta2 = plt.axes([0.28, 0.45, 0.12, 0.02], facecolor=axcolor)
ax_K2 = plt.axes([0.28, 0.4, 0.12, 0.02], facecolor=axcolor)
ax_mumax2 = plt.axes([0.28, 0.35, 0.12, 0.02], facecolor=axcolor)

sYcs1 = Slider(ax_Ycs1, r'$Y_{c/s1}$', 0.05, 0.8, valinit=0.4, valfmt = "%1.2f")
salpha1 = Slider(ax_alpha1, r'$\alpha_1 (\frac{g_{enzyme1}}{g_{cell mass}.s})$', 5e-6, 0.001, valinit=0.0001, valfmt = "%1.4f")
sbeta1 = Slider(ax_beta1, r'$\beta_1 (hr^{-1})$', 0.01, 0.2, valinit=0.05, valfmt = "%1.2f")
sK1 = Slider(ax_K1, r'$K_1 (\frac{g}{L})$',0.02, 0.3, valinit= 0.1, valfmt = "%1.2f")
smumax1 = Slider(ax_mumax1, r'$\mu_{1max} (hr^{-1})$', 0.2, 1.5, valinit=0.9, valfmt = "%1.2f")
sYcs2 = Slider(ax_Ycs2, r'$Y_{c/s2}$', 0.05, 1.5, valinit=0.4, valfmt = "%1.2f")
salpha2 = Slider(ax_alpha2, r'$\alpha_2 (\frac{g_{enzyme2}}{g_{cell mass}.s})$', 5e-6, 0.001, valinit=0.0001, valfmt = "%1.4f")
sbeta2 = Slider(ax_beta2, r'$\beta_2 (hr^{-1})$', 0.01, 0.5, valinit=0.05, valfmt = "%1.2f")
sK2 = Slider(ax_K2, r'$K_2 (\frac{g}{L})$',0.02, 1, valinit= 0.5, valfmt = "%1.2f")
smumax2 = Slider(ax_mumax2, r'$\mu_{2max} (hr^{-1})$', 0.2, 2, valinit=0.6, valfmt = "%1.2f")


mu1max =  .9  
mu2max =  .6  
K1 =  .1  
K2 =  .5  
alpha1 =  .0001  
alpha2 =.0001   
beta1 =  .05  
beta2 =  .05  
Ycs1 =  .4  
Ycs2 =  .4

def update_plot2(val):
    Ycs1 = sYcs1.val
    alpha1 = salpha1.val
    beta1 =sbeta1.val
    K1 = sK1.val
    mu1max = smumax1.val
    Ycs2 = sYcs2.val
    alpha2 = salpha2.val
    beta2 =sbeta2.val
    K2 = sK2.val
    mu2max = smumax2.val    
    sol = odeint(ODEfun, y0, tspan, (mu1max, mu2max, K1, 
                                     K2, alpha1, alpha2, beta1, beta2,
                                     Ycs1, Ycs2))
    Cc = sol[:,0]
    Cs1= sol[:,1]
    Cs2= sol[:,2]
    e1 = sol[:,3]
    e2 = sol[:,4]
    p1.set_ydata(Cc)
    p2.set_ydata(Cs1)
    p3.set_ydata(Cs2)
    p4.set_ydata(e1)
    p5.set_ydata(e2)    
    fig.canvas.draw_idle()


sYcs1.on_changed(update_plot2)
salpha1.on_changed(update_plot2)
sbeta1.on_changed(update_plot2)
sK1.on_changed(update_plot2)
smumax1.on_changed(update_plot2)
sYcs2.on_changed(update_plot2)
salpha2.on_changed(update_plot2)
sbeta2.on_changed(update_plot2)
sK2.on_changed(update_plot2)
smumax2.on_changed(update_plot2)

resetax = plt.axes([0.29, 0.85, 0.09, 0.04])
button = Button(resetax, 'Reset variables', color='cornflowerblue', hovercolor='0.975')

def reset(event):
    sYcs1.reset()
    salpha1.reset()
    sbeta1.reset()
    sK2.reset()
    smumax2.reset()
    sYcs2.reset()
    salpha2.reset()
    sbeta2.reset()
    sK2.reset()
    smumax2.reset()    
button.on_clicked(reset)