#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 12 17:37:29 2022
%problem3_7_5.py
%Problem 3.7.5 part (a)


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



#assign parameter values
v0=2
Vm1=9
Vm2=12
Vm3=15
Km1=1
Km2=0.4
Km3=3


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

###Full model#####
#set initial conditions State vector is [s1, s2, s3]
S0=[0.3,0.2,0.1];
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=[0,0,0] #generate a list to store derivatives
    dS[0]=v0 - Vm1*S[0]/(Km1+S[0]) 
    dS[1]=Vm1*S[0]/(Km1+S[0])-Vm2*S[1]/(Km2+S[1])
    dS[2]=Vm2*S[1]/(Km2+S[1])-Vm3*S[2]/(Km3+S[2])
    return dS


S=odeint(dSdt_original, S0, 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.plot(times, S[:,2], label="s3", linewidth=2)
plt.xlabel("Time")
plt.ylabel("concentration")
plt.legend()