#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
%file pulse_generator.py
%model of synthetic pulse generating system
%adapted from Basu et al. (2004) PNAS 101 6355-6360
%Figure 7.28

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



#assign parameter values

aG = 80 # muM/min
bG = 0.07 # /min 
KC= 0.008 # muM
aC= 0.5 # muM/min
bC = 0.3 # /min 
k1 = 0.5 # /muM^3 /min
k2 = 0.02 # /min 
KR= 0.02 # muM
RT = 0.5 # muM



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


#assign initial condition.  
#state vector is [G C R]
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
    
    G=S[0]
    C=S[1]
    R=S[2]
    
    
    if t < 0:
        A=0
    else:
        A=10
    
    dS[0]=- bG*G + aG*((R/KR)/(1+(R/KR) + (C/KC)**2 + (R/KR)*(C/KC)**2))
    dS[1]=- bC*C + aC*R/(KR+R)
    dS[2]=- k2*R + A**2*k1*np.power([RT-2*R],2)

    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()