#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
% file predator_prey.py
% model of predator-prey network
% from Balagadde et al. (2008) Molecular Systems Biology 4
% problem 7.8.20 (b)


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



#assign parameter values

k2=0.4
cm=1e5
d1=1
K1=10
beta=2
D=0.2
gamma1=0.1
gamma2=0.1
del1=0.017
del2=0.11
k1=0.8
d2=0.3

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


#assign initial condition.  
#state is [c1 c2 a1 a2]
Sinit=[10000, 10000, 1, 1]


 
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(4) #generate a list to store derivatives
    
    c1=S[0]
    c2=S[1]
    a1=S[2]
    a2=S[3]
    
    dS[0]= k1*c1*(1-(c1+c2)/cm)-d1*c1*K1/(K1+a2**beta) - D*c1
    dS[1]= k2*c2*(1-(c1+c2)/cm)-d2*c2*a1**beta/(K1+a1**beta) - D*c2
    dS[2]= gamma1*c1-(del1+D)*a1
    dS[3]= gamma2*c2-(del2+D)*a2
    
    
    
    return dS


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

#plot simulation
plt.figure() #generate figure 7.17A
plt.plot(times, S[:,0], label="c1", linewidth=2)
plt.plot(times, S[:,1], label="c2", linewidth=2)
plt.xlabel("Time")
plt.ylabel("size")
plt.legend()