#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 12 17:37:29 2022
% file lac_operon.py
% model of lac operon in E. coli
% from Santillan et al. (2007) Biophys. J. 92, pp. 3830-3842
% Figure 7.7




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



#assign parameter values

delta_M=0.48
delta_Y=0.03
delta_L=0.02
a1=0.29
K2=2.92*1e6
RToverK1=213.2
c1=18.8
kL=6*1e4
KML=680
kg=3.6*1e3
KMg=7*1e5
Le=0

flag=1



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


#set initial conditions State vector is  [M Y L]
Sinit=[0.01, 0.1, 0]
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(3) #generate a list to store derivatives
    
    M=S[0]
    Y=S[1]
    L=S[2]
    
    
    #dependent states
    B=Y/4
    A=L
    
    Lelocal=Le
    #time-varying input profile
    if flag==1:
        Lelocal=0
        if t<500:
            Lelocal=0.00
        elif t<1000:
            Lelocal=50
        elif t<1500:
            Lelocal=100
        elif t<2000:
            Lelocal=150
  
    
    dS[0]=  a1*1/(1+RToverK1*(K2/(K2+A))**4) - delta_M*M;
    dS[1]= c1*M - delta_Y*Y;
    dS[2]= kL*Y*Lelocal/(KML+Lelocal) - 2*kg*B*L/(KMg+L) - delta_L*L;
  

    return dS




#plot transient Figure 7.7A

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

#generate input profile for plotting
t_timecourse=[0,499.99, 500, 999.99, 1000, 1499.99, 1500, 1999.99, 2000, 2500 ];
Le_timecourse=[0, 0, 50, 50, 100, 100, 150, 150, 0, 0];



plt.figure() #generate figure
plt.plot(times, S[:,1]/4, label="", linewidth=2)
plt.plot(t_timecourse, Le_timecourse,label="", linewidth=2)
plt.xlabel("Time (min)")
plt.ylabel("beta-galactosidase concentration (molecules per cell)")
plt.legend()
plt.xlim(0,2500)
plt.ylim(0,200)




#generate dose-response  Figure 7.7B
flag=2;

#dynamics for modified model 
def dSdt_modified(S,t):
    dS=np.zeros(3) #generate a list to store derivatives
    
    M=S[0]
    Y=S[1]
    L=S[2]
    
    
    #dependent states
    B=Y/4
    A=L
    
    Enz=40;
  
    
    dS[0]=  a1*1/(1+RToverK1*(K2/(K2+A))**4) - delta_M*M;
    dS[1]= c1*M - delta_Y*Y;
    dS[2]= kL*4*Enz*Le/(KML+Le) - 2*kg*Enz*L/(KMg+L) - delta_L*L;
  

    return dS


#run simulations to steady state over a range of fixed inputs
N=180;ArithmeticErrorcompYVec=np.zeros(N);
compLeVec=np.zeros(N)
compYVec=np.zeros(N)
pertYVec=np.zeros(N)
pertLeVec=np.zeros(N)
for i in range(N):
    Le=120*(i-1)/N;
    S=odeint(dSdt_original, Sinit, times)
    compLeVec[i]=Le;
    compYVec[i]=S[len(times)-1,1]
    S=odeint(dSdt_modified, Sinit, times)
    pertLeVec[i]=Le;
    pertYVec[i]=S[len(times)-1,1]


#generate Figure 7.7B`
plt.figure() #generate figure
plt.plot(compLeVec,compYVec/4, label="original", linewidth=2)
plt.plot(pertLeVec,pertYVec/4,label="modified", linewidth=2)
plt.xlabel("External lactose concentration (\muM)")
plt.ylabel("beta-galactosidase concentration (molecules per cell)")
plt.legend()
plt.xlim(0,100)
plt.ylim(0,100)