#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 12 17:37:29 2022
%Model of apoptosis signalling pathway
%modified from Eissing et al. (2004) Journal of Biological Chemistry 279, pp.
%36892-36897
%Figure 6.16

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



#assign parameter values

k1=507
k2=3.9e-3
k3=1e-5 
k4=5.8e-3
k5=5e-4
k6=0.21
k7=81.9
k8=3.9e-3; 
k9=5.8e-6
k10=5.8e-3
k11=5e-4
k12=0.21
k13=40
k14=1e-3
k15=464
k16=1.16e-2
k17=3e-4
k18=1.16e-2;
k19=1.73e-2

input=0


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


#assign initial condition.  These values were taken from a previous
#simulation run to steady state with input=0
#state vector S=[C8 C8star C3 C3star BAR IAP C8starBAR C3starIAP]
Sinit=[ 1e5*1.30, 1e5*0, 1e5*0.21, 1e5*0, 1e5*0.4, 1e5*0.4, 1e5*0, 1e5*0]
 
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(8) #generate a list to store derivatives
    
    C8=S[0]
    C8star=S[1]
    C3=S[2]
    C3star=S[3]
    BAR=S[4]
    IAP=S[5]
    C8starBAR=S[6]
    C3starIAP=S[7]
    
    #assign time-varying profile for parameter 'input'
    input_local=input;
    if t>100 and t<100+1100:
       input_local=200

    
    dS[0]=k1 - k2*C8 - k3*(C3star+input_local)*C8
    dS[1]=k3*(C3star+input_local)*C8- k4*C8star - k5*C8star*BAR + k6*C8starBAR
    dS[2]=k7-k8*C3 - k9*C8star*C3
    dS[3]=k9*C8star*C3 - k10*C3star-k11*C3star*IAP+k12*C3starIAP
    dS[4]=k13 - k5*C8star*BAR + k6*C8starBAR - k14*BAR
    dS[5]=k15 - k11*C3star*IAP+k12*C3starIAP - (k16+k17*C3star)*IAP
    dS[6]=k5*C8star*BAR - k6*C8starBAR - k18*C8starBAR
    dS[7]=k11*C3star*IAP - k12*C3starIAP - k19*C3starIAP

    return dS


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

#plot simulation
plt.figure() #generate figure
plt.plot(times, S[:,1], label="C8*", linewidth=2)
plt.plot(times, S[:,3], label="C3*", linewidth=2)
plt.xlabel("Time")
plt.ylabel("concentration")
plt.legend()