#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 12 17:37:29 2022
% file camp_signaling.py
% model of cAMP signaling pathway in Dictyostelium
% described in Maeda et al (2004) Science 304 pp. 875-878. 
% problem 6.8.15




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



#assign parameter values

k1=2.0
k2=0.9
k3=2.5
k4=1.5
k5=0.6
k6=0.8
k7=1.0
k8=1.3
k9=0.3
k10=0.8
k11=0.7
k12=4.9
k13=23.0
k14=4.5



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


#set initial conditions State vector is 
#MKKK, MKKKP, MKK, MKKP, MKKPP, MK, MKP, MKPP
Sinit=[1,1,1,1,1,1,1]
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(7) #generate a list to store derivatives
    
    ACA=S[0]
    PKA=S[1]
    ERK2=S[2]
    REGA=S[3]
    CAMPi=S[4]
    CAMPe=S[5]
    CAR1=S[6]
    
    dS[0]=k1*CAR1 - k2*ACA*PKA;
    dS[1]=k3*CAMPi-k4*PKA;
    dS[2]=k5*CAR1-k6*ERK2*PKA;
    dS[3]=k7-k8*REGA*ERK2;
    dS[4]=k9*ACA-k10*REGA*CAMPi;
    dS[5]=k11*ACA-k12*CAMPe;
    dS[6]=k13*CAMPe-k14*CAR1;


    return dS


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

#plot transient
plt.figure() #generate figure
plt.plot(times, S[:,4], label="CAMPi", linewidth=2)
plt.plot(times, S[:,1], label="PKA", linewidth=2)
plt.plot(times, S[:,2], label="ERK2", linewidth=2)
plt.xlabel("Time")
plt.ylabel("concentration")
plt.legend()