#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
%file band_detector.py
%model of synthetic band detector system
%adapted from Basu et al. (2005) Nature 434 177-193.
%Figure 7.30

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



#assign parameter values

aG = 2 # muM/min
bG = 0.07 # /min 
KL= 0.8 # muM
aL1= 1 # muM/min
aL2= 1 # muM/min
KC= 0.008 # muM
bC = 0.07 # /min 
k1 = 0.5 # /muM^3 /min
k2 = 0.02 # /min 
KR= 0.01 # muM
RT = 0.5 # muM
bL = 0.02 # /min 
aC= 1 # muM/min
A=0.1



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


#assign initial condition.  
#state vector is [G L C R]
Sinit=[ 0,0,0,0]
 
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(4) #generate a list to store derivatives
    
    G=S[0]
    L=S[1]
    C=S[2]
    R=S[3]
    
    
    dS[0]=-bG*G + aG*(1/(1+(L/KL)**2))
    dS[1]=-bL*L + aL1/(1+(C/KC)**2) + aL2*R/(KR+R)
    dS[2]=-bC*C + aC*R/(KR+R)
    dS[3]=-k2*R + A**2*k1*(RT-2*R)**2

    return dS



N=90  #set grid-size 
Amesh=np.zeros(N)
gmesh=np.zeros(N)
cmesh=np.zeros(N)
lmesh=np.zeros(N)
 
 #over the grid of inputs, simulate and record steady-state values
for i in range(N):
     i
     A=10**(-4 + 4*(i-1)/N);
     S=odeint(dSdt_original, Sinit, times) #run simulation
     Amesh[i] = A;
     gmesh[i]=S[len(times)-1,0]
     lmesh[i]=S[len(times)-1,1]
     cmesh[i]=S[len(times)-1,2]

 
 #produce figure 7.30
plt.figure()
plt.semilogx(Amesh, gmesh, label="GFP", linewidth=2)
plt.semilogx(Amesh, lmesh, label="LacI", linewidth=2)
plt.semilogx(Amesh, cmesh,label="cI", linewidth=2)
plt.xlabel('Extracellular AHL concentration (\mu M)')
plt.ylabel('Intracellular concentration (\mu M)')
plt.legend()