#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 12 17:37:29 2022
%Model of E. coli chemotaxis signalling pathway
%Figure 6.14


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



#assign parameter values

k1=200
k2=1
k3=1
km1=1
km2=1
km3=1 
k5=0.05
km5=0.005
k4=1
km4=1
KM1=1
KM2=1
cheR=1
L=20





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



#assign initial condition species vector S=[Am AmL A AL B BP]
#these values were determined by running a previous simulation to steady
#state with L=20.
Sinit=[0.0360,    1.5593,    0.0595 ,   0.3504 ,   0.7356 ,   0.2644]

#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(6) #generate a list to store derivatives
    
    #assign state variables
    Am=S[0]
    AmL=S[1]
    A=S[2]
    AL=S[3]
    B=S[4]
    BP=S[5]
    
    Llocal=L
    #set up time-varying ligan profile
    if (t>10):
        Llocal=40
    if (t>30):
        Llocal=80

    
    dS[0]=km1*cheR - (k1*BP*Am)/(KM1 + Am) - k3*Am*Llocal + km3*AmL
    dS[1]=km2*cheR - (k2*BP*AmL)/(KM2 + AmL) + k3*Am*Llocal - km3*AmL
    dS[2]=-km1*cheR + (k1*BP*Am)/(KM1 + Am) - k4*A*Llocal + km4*AL
    dS[3]=-km2*cheR + (k2*BP*AmL)/(KM2 + AmL) + k4*A*Llocal - km4*AL
    dS[4]=-(k5*Am*B) +  (km5*BP)
    dS[5]=(k5*Am*B) -  (km5*BP)
    return dS


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

#set up ligand profile for plotting
lt_toplot=[0, 10, 10, 30, 30, 50]
l_toplot=[0.012, 0.012, 0.014, 0.014, 0.018, 0.018]


#plot simulation
plt.figure() #generate figure
plt.plot(times, S[:,0], linewidth=2)
plt.plot(lt_toplot, l_toplot, label="ligand level", linewidth=2)
plt.xlabel("Time (arbitrary units)")
plt.ylabel("Concentration of active CheA [Am] (arbitrary units)")
plt.legend()
plt.xlim(0,50)
plt.ylim(0.01, 0.04)