#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
%file Goodwin.py
%Goodwin oscillator model
%from Goodwin (1965) Adv. in. Enz. Reg. 3 pp. 425-428.
%Figure 7.17

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



#assign parameter values

a1=360
kappa1=43
k1=1
b1=1
alpha1=1
beta1=0.6
gamma1=1
delta1=0.8
n=12


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


#assign initial condition.  
#state vector S=[X Y Z]
Sinit=[ 0,0,0]
 
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(3) #generate a list to store derivatives
    
    X=S[0]
    Y=S[1]
    Z=S[2]
    
    dS[0]=a1/(kappa1+k1*Z**n) - b1*X
    dS[1]=alpha1*X - beta1*Y
    dS[2]=gamma1*Y - delta1*Z

    return dS


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

#plot simulation
plt.figure() #generate figure 7.17A
plt.plot(times, S[:,0], label="mRNA (X)", linewidth=2)
plt.plot(times, S[:,1], label="enzyme (Y)", linewidth=2)
plt.plot(times, S[:,2], label="metabolite (Z)", linewidth=2)
plt.xlabel("Time")
plt.ylabel("concentration")
plt.legend()

#generate figure 7.17B

plt.figure() 


ax = plt.axes(projection='3d')

# Data for a three-dimensional line
zline = S[:,0]
xline = S[:,1]
yline = S[:,2]
ax.plot3D(xline, yline, zline, 'gray')

ax.view_init(-160, 160)