#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 12 17:37:29 2022
%file gprotein_pathway.m
%Model of g-protein signalling pathway
%adapted from Yi et al. 2003, PNAS 100, pp. 10764-10769 
%Figure 6.5


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

#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.5B)
dose_response_trigger=0;

#parameter assignment

krl=2*1e6;
krlm=1e-2;
krs=4
krd0=4*1e-4
krd1=2*1e-2


kg1=1; 


kga=1e-5;
kgd1=0.11;


#Total populations
gt=1e4;
rt=4000;

lt=10^(-9)




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


#set initial conditions The state vector is S=[RL, G, Ga]
Sinit=[0,10000,0]

#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=[0,0,0] #generate a list to store derivatives
    #locally define state variables:
    rl=S[0]
    g=S[1]
    ga=S[2]
    #conservations:
    r=rt-rl
    gbg=gt-g
    gd=gt-g-ga
    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 > 100 and t< 700:
    ltlocal=10^(-9)
   

    #concentration kinetics:
    ddt_rl= krl*ltlocal*r - krlm*rl
    ddt_g= -kga*rl*g + kg1*gd*gbg
    ddt_ga= kga*rl*g - kgd1*ga

    dS[0]=ddt_rl
    dS[1]=ddt_g
    dS[2]=ddt_ga

    return dS



#OPTIONS = odeset('RelTol',1e-6,'AbsTol',1e-9, 'refine',5);

S=odeint(dSdt_original, Sinit, times, rtol=1e-18, atol=1e-18) #run simulation


#plot simulation
plt.figure() #generate figure
plt.plot(times, S[:,2], label="Active Galpha-GTP (Ga)", 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.5B

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