#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
% file deterministic_vilar_oscillator.m
% deterministic simulation of genetic oscillator
% from Vilar et al. (2002) PNAS 99 5988-5992.
% problem 7.8.27

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



#assign parameter values

eps=0.1
ga=250
ba=5
Ka=0.5
alpha0=0.1
da=1
kc=200
gr=50
br=10
Kr=1
dr=eps*da


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


#assign initial condition.  
#state vector S=[A R C]
Sinit=[5, 10, 35]
 
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(3) #generate a list to store derivatives
    
    A=S[0]
    R=S[1]
    C=S[2]
    
    dS[0]=ba*(ga/ba)*(alpha0+(A/Ka))/(1+(A/Ka)) - da*A - kc*A*R
    dS[1]=br*(gr/br)*((A/Kr))/(1+(A/Kr)) - dr*R - kc*A*R + da*C
    dS[2]=kc*A*R - da*C
   

    return dS


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

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