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


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



#assign parameter values

e0=1
e1=1.5
e2=2
S0=1
km1=1
km2=1
km3=1
km4=1
km5=1
v1=1
v2=1
v3=2
v4=2
v5=0.5



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


#set initial conditions State vector is [s1, s2, s3]
Sinit=[0,0];
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=[0,0] #generate a list to store derivatives
    
    s1=S[0]
    s2=S[1]
    dS[0]=e0*(v1*S0-v2*s1)/(1+S0/km1+s1/km2)- e1*(v3*s1-v4*s2)/(1+s1/km3+s2/km4)
    dS[1]=e1*(v3*s1-v4*s2)/(1+s1/km3+s2/km4) - e2*(v5*s2)/(1+s2/km5)
    return dS


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

#plot simulation
plt.figure() #generate figure
plt.plot(times, S[:,0], label="s1", linewidth=2)
plt.plot(times, S[:,1], label="s2", linewidth=2)
plt.xlabel("Time")
plt.ylabel("concentration")
plt.legend()