#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 12 17:37:29 2022
%ptwo_component_pathway.py
%generates figure 6.3


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



#assign parameter values
#kinetic parameters
k1=5
km1=1
k2=6
k3=2

#concentration totals
RT=2
LT=3
PT=8

#trigger=0 to simulate time-series with varying ligand (Figure 6.3A)
#trigger set to 1 below to simulate multiple steady states for
#dose-response curve (Figure 6.3B)
dose_response_trigger=0;



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


#set initial conditions State vector is [s1, s2, s3]
Sinit=[0,PT];
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=[0,0] #generate a list to store derivatives
    
    RL=S[0]
    P=S[1]
    LTlocal=LT
    
    #if the trigger is zero, set a time-varying ligand profile
    #(otherwise, the ligand profile will be passed from the main function)
    if dose_response_trigger < 1:
        LTlocal=0
        if t> 1 and t < 3:
            LTlocal=3
    #determine conserved species
    R=RT-RL
    Pstar=PT-P
    #assign dynamics
    RLdot=k1*R*LTlocal - km1*RL
    Pdot=-k2*RL*P+k3*Pstar
    dS[0]=RLdot
    dS[1]=Pdot
    return dS




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

#determine conserved species for plotting
R=RT-S[:,0];
Pstar=PT-S[:,1];


#plot simulation
plt.figure() #generate figure
plt.plot(times, Pstar, label="Active Response Protein (P*)", linewidth=2)
plt.plot(times, S[:,0], label="Receptor-Ligand Complex (RL)", linewidth=2)
plt.plot([0,0.99999,1,2.99999,3,5], [0,0,3,3,0,0], 'k:', label="Total Ligand (L_T)", linewidth=2)
plt.xlabel("Time")
plt.ylabel("concentration")
plt.legend()


#generate figure 6.3B

#generate dose response curve

#turn off time-varying lingand profile
dose_response_trigger=1


#generate mesh 
N=41
LTmesh=np.zeros(N)
Pstarmesh=np.zeros(N)
RLmesh=np.zeros(N)

#run multiple simulations to steady state
times=np.linspace(0, 12, 1000)
for j in range(N):

    LT=1*(j-1)/(N-1)
    S=odeint(dSdt_original, Sinit, times) #run simulation

    Pstar=PT-S[:,1];
    LTmesh[j]=LT
    Pstarmesh[j]=Pstar[len(times)-1]
    RLmesh[j] = S[len(times)-1,0]


#generate figure
plt.figure()
plt.plot(LTmesh, Pstarmesh, 'k', label="Active Response Protein (P*)",linewidth=2)
plt.plot(LTmesh, RLmesh, 'k-.', label="Receptor-Ligand Complex (RL)", linewidth=2)
plt.xlabel("Total Ligand Concentration (L_T) (arbitrary units)")
plt.ylabel("Steady-state concentration (arbitrary units)")
plt.xlim(0,1)
plt.ylim(0,9)