#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 12 17:37:29 2022
% file mapk_pathway.m
% model of MAPK signalling pathway
% from Kholodenkho (2000) Euro. J. Biochem. 267 pp 1583-1588
% problem 6.8.8

function mapk_pathway
%from Kholodenkho, 2000



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



#assign parameter values

KA=0
KI=0
V1=2.5
KM1=10
V2=0.25
KM2=8
kcat3=0.025
KM3=15
kcat4=0.025
KM4=15
V5=0.75
KM5=15
V6=0.75
KM6=15
kcat7=0.025
KM7=15
kcat8=0.025
KM8=15
V9=0.5
KM9=15
V10=0.5
KM10=15



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


#set initial conditions State vector is 
#MKKK, MKKKP, MKK, MKKP, MKKPP, MK, MKP, MKPP
Sinit=[100, 0, 300, 0, 0, 300, 0, 0]
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(8) #generate a list to store derivatives
    
    MKKK=S[0]
    MKKKP=S[1]
    MKK=S[2]
    MKKP=S[3]
    MKKPP=S[4]
    MK=S[5]
    MKP=S[6]
    MKPP=S[7]
    
    r1=V1*MKKK*(1+KA*MKPP)/((1+(KI*MKPP))*(KM1 + MKKK))
    r2=V2*MKKKP/(KM2+MKKKP)
    r3=kcat3*MKKKP*MKK/(KM3+MKK)
    r4=kcat4*MKKKP*MKKP/(KM4+MKKP)
    r5=V5*MKKPP/(KM5+MKKPP)
    r6=V6*MKKP/(KM6+MKKP)
    r7=kcat7*MKKPP*MK/(KM7+MK)
    r8=kcat8*MKKPP*MKP/(KM8+MKP)
    r9=V9*MKPP/(KM9+MKPP)
    r10=V10*MKP/(KM10+MKP)


    dS[0]=r2-r1
    dS[1]=r1-r2
    dS[2]=r6-r3
    dS[3]=r3+r5-r4-r6
    dS[4]=r4-r5
    dS[5]=r10-r7
    dS[6]=r7+r9-r8-r10
    dS[7]=r8-r9

    return dS


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

#plot transient
plt.figure() #generate figure
plt.plot(times, S[:,7], label="MKPP", linewidth=2)
plt.xlabel("Time")
plt.ylabel("concentration")
plt.legend()




#generate dose-reponse.  

M=50;

input=np.zeros(M)
output=np.zeros(M)

times=np.linspace(0, 8000, 1000) 


for i in range(M):
    V1=0.00+(i/M)*0.4;
    S=odeint(dSdt_original, Sinit, times) 
    input[i]=V1;
    output[i]=S[len(times)-1,7];


plt.figure()
plt.plot(input, output, 'k', linewidth=2)
plt.title('dose-response')