#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 12 17:37:29 2022
% file calcium_oscillations.m
% model of calcium-induced calcium release in hepatocytes
% adapted from Othmer and Tang (1997) in 'Case studies in mathematical
% modelling', Othmer et al. eds, Prentice Hall
% Figure 6.18



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



#assign parameter values

gamma0=0.1
gamma1=20.5
p1=8.5
p2=0.065
k1=12
k2=15
k3=1.8
km1=8
km2=1.65
km3=0.21
Cs=8.37
vr=0.185
I=1




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


#set initial conditions  State vector is S = [ C R RI RIC RICC]
Sinit=[0,1,0,0,0];
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(5) #generate a list to store derivatives
    
    C=S[0]
    R=S[1]
    RI=S[2]
    RIC=S[3]
    RICC=S[4]
    
    Ilocal=I
    
    #if mode=2, set the input parameter according to a time-varying profile:

    if mode==2:
        Ilocal=0
        if t>20:
            Ilocal=0.7
        if t>60:
            Ilocal=1.2
        if t>90:
            Ilocal=4

    
    dS[0]=vr*(gamma0+gamma1*RIC)*(Cs-C) - (p1*C**4)/(p2**4+C**4)
    dS[1]=-k1*Ilocal*R+km1*RI
    dS[2]=-(km1+k2*C)*RI+ k1*Ilocal*R + km2*RIC
    dS[3]=-(km2+k3*C)*RIC + k2*C*RI + km3*RICC
    dS[4]=k3*C*RIC - km3*RICC
    
    return dS


#plot simulation: Figure6.18A

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

plt.figure() #generate figure
plt.plot(times, S[:,0], label="C", linewidth=2)
plt.plot(times, S[:,3], label="RIC", linewidth=2)
plt.plot(times, S[:,4], label="RICC", linewidth=2)
plt.xlabel("Time")
plt.ylabel("concentration")
plt.legend()
plt.xlim(0,25)
plt.ylim(0,2.5)


#plot simulation: Figure6.18B
# mode 2: varying input
mode=2

#set initial condition
Sinit=[0, 1 ,0 ,0 ,0];
times=np.linspace(0,120,1000)

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

#generate Figure 6.18B
plt.figure() #generate figure
plt.plot(times, S[:,0], label="cytosolic [$Ca^{2+}$] ($\mu$ M)", linewidth=2)
plt.xlabel("Time (s)")
plt.ylabel("concentration")
plt.legend()
plt.xlim(0,120)
plt.ylim(0,2.5)