#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
% file population_control.py
% model of population control network
% from You et al. (2004) Nature 428 pp. 868-871 
% problem 7.8.20 (a)


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



#assign parameter values

k=0.97
Nm=1.24e9
d=0.004
ke=5
de=2
va=4.8e-7
da=0.64

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


#assign initial condition.  
#state is [N, E, A]
Sinit=[1,1,1]

 
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(3) #generate a list to store derivatives
    
    N=S[0]
    E=S[1]
    A=S[2]
    
    dS[0]= k*N*(1-N/Nm)-d*E*N
    dS[1]= ke*A-de*E
    dS[2]= va*N-da*A
    
    
    return dS


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

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