#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 12 17:37:29 2022
%qmichaelis_menten.py
%Figure 3.3 Michaelis-Menten kinetics


@author: bingalls
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint



#assign parameter values
k1=30
km1=1
k2=10
et=1


#set time_grid for simulation
t_min=0; t_max=1; dt=0.001
times=np.arange(t_min, t_max+dt, dt) #generate time-grid list

###Full model#####
#set initial conditions State vector is [s, c, p]
S0=[5,0,0];
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=[0,0,0] #generate a list to store derivatives
    dS[0]=-k1*S[0]*(et-S[1]) + km1*S[1]
    dS[1]=-km1*S[1] + k1*S[0]*(et-S[1])-k2*S[1]
    dS[2]=k2*S[1]
    return dS


S=odeint(dSdt_original, S0, times) #run simulation

#plot simulation
plt.figure() #generate figure
plt.plot(times, S[:,0], label="s", linewidth=2)
plt.plot(times, S[:,1], label="c", linewidth=2)
plt.plot(times, S[:,2], label="p", linewidth=2)
plt.plot(times, et-S[:,1], label="e", linewidth=2)
plt.xlabel("Time")
plt.ylabel("concentration")
plt.legend()




###Approximate model#####
#set initial conditions for approximate model
Sred0=[5];



#declare right-hand-side for approximate model
def dSdt_approx(S,t):
    dS=[0] #generate a list to store derivatives
    dS=[-k1*k2*et*S[0]/(km1+k2+k1*S[0])]
    return dS

S_red=odeint(dSdt_approx, Sred0, times) #run simulation


#plot simulation
plt.figure() #generate figure
plt.plot(times, S[:,0], label="s (full model)", linewidth=2)
plt.plot(times, S[:,2], label="p (full model)", linewidth=2)
plt.plot(times, S_red[:,0], label="s (reduced model)", linewidth=2)
plt.plot(times, Sred0-S_red[:,0], label="p (reduced model)", linewidth=2)
plt.xlabel("Time")
plt.ylabel("concentration")
plt.legend()