#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 12 17:37:29 2022
%parabidopsis_branch.py
%Problem 5.6.7


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



#assign parameter values

kcatCGS=30  # /sec
KmCGScys=460 # muM
KmCGSPhser=15000 # muM
KiCGSPi=2000 # muM
 
JPhser=0.3 # muM/sec
 
Cys=250 # muM
CGS=0.7 # muM
ADOMET=20 # muM
TS=5 # muM




#set time_grid for simulation
times=np.linspace(0, 5000, 1000) #generate time-grid list


#set initial conditions State vector is [s1, s2, s3]
Sinit=[150];
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=[0] #generate a list to store derivatives
    Phser=S[0]

    kcatappCGS=kcatCGS/(1+KmCGScys/Cys)

    KmCGSapp=(KmCGSPhser/(1+KmCGScys/Cys))

    kTS=4.9*1e-6+ 5.6*1e-4*np.power(ADOMET,2.9)/(np.power(32,2.9)+np.power(ADOMET,2.9))

    vThr=TS*kTS*Phser
 
    vcystathionine=kcatappCGS*CGS*Phser/(KmCGSapp+Phser)

    dS=[JPhser - vThr-vcystathionine]

    return dS


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

# to plot fluxes
kcatappCGS=kcatCGS/(1+KmCGScys/Cys)

KmCGSapp=(KmCGSPhser/(1+KmCGScys/Cys))

kTS=(1/11)*(5.4*1e-5+ 6.2*1e-3*np.power(ADOMET,2.9)/(np.power(32,2.9)+np.power(ADOMET,2.9)))

vThr=TS*kTS*S
 
vcystathionine=kcatappCGS*CGS*S/(KmCGSapp+S);
 

#plot simulation
plt.figure() #generate figure
plt.plot(times, vThr, label="vThr", linewidth=2)
plt.plot(times, vcystathionine, label="vcystathionine", linewidth=2)
plt.xlabel("Time")
plt.ylabel("concentration")
plt.legend()