#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""

% file phage_lambda.m\py
% model of phage lambda decision switch
% Figure 7.11

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



#assign parameter values

K1=1
K2=0.1
K3=5
K4=0.5 
delta_r=0.02
delta_c=0.02
a=5
b=50



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

#declare right-hand-side for model
def dSdt_original(S,t):
    dS=[0,0] #generate a list to store derivatives
    
    r=S[0]
    c=S[1]
    rd=r/2
    cd=c/2
    
    dS[0]=(a + 10*a*K1*rd**2)/(1+K1*rd**2+ K1*K2*rd**3+K3*cd+K3*K4*cd**2)- delta_r*r
    dS[1]=(b + b*K3*cd)/(1+K1*rd**2+ K1*K2*rd**3+K3*cd+K3*K4*cd**2) - delta_c*c
    return dS



#generate  figure 7.11A

S1_0=[250,50]
S2_0=[250,150]
S3_0=[250,250]
S4_0=[150,250]
S5_0=[100,250]
S6_0=[50,250]
S7_0=[25,250]
S8_0=[17,250]
S9_0=[3,2]
S10_0=[6,25]


S1=odeint(dSdt_original, S1_0, times) #run simulation
S2=odeint(dSdt_original, S2_0, times)
S3=odeint(dSdt_original, S3_0, times)
S4=odeint(dSdt_original, S4_0, times)
S5=odeint(dSdt_original, S5_0, times)
S6=odeint(dSdt_original, S6_0, times)
S7=odeint(dSdt_original, S7_0, times)
S8=odeint(dSdt_original, S8_0, times)
S9=odeint(dSdt_original, S9_0, times)
S10=odeint(dSdt_original, S10_0, times)



plt.figure() 
plt.plot(S1[:,0], S1[:,1], linewidth=2)
plt.plot(S1[:,0], S1[:,1], linewidth=2)
plt.plot(S2[:,0], S2[:,1], linewidth=2)
plt.plot(S3[:,0], S3[:,1], linewidth=2)
plt.plot(S4[:,0], S4[:,1], linewidth=2)
plt.plot(S5[:,0], S5[:,1], linewidth=2)
plt.plot(S6[:,0], S6[:,1], linewidth=2)
plt.plot(S7[:,0], S7[:,1], linewidth=2)
plt.plot(S8[:,0], S8[:,1], linewidth=2)
plt.plot(S9[:,0], S9[:,1], linewidth=2)
plt.plot(S10[:,0], S10[:,1], linewidth=2)
plt.xlabel("cro concentration of (nM)")
plt.ylabel("cI oncentration of (nM)")
plt.legend()
plt.title('Figure 7.11A')


# xrange = arange(0, 2, delta)
# yrange = arange(0, 2, delta)
xrange = np.arange(0, 350, 5)
yrange = np.arange(0, 350, 5)
r, c = np.meshgrid(xrange,yrange)


F = (a + 10*a*K1*(r/2)**2)/(1+K1*(r/2)**2+ K1*K2*(r/2)**3+K3*(c/2)+K3*K4*(c/2)**2)- delta_r*(r)
G = (b + b*K3*(c/2))/(1+K1*(r/2)**2+ K1*K2*(r/2)**3+K3*(c/2)+K3*K4*(c/2)**2) - delta_c*(c)

#plot nullclines
plt.contour(r,c, F, [0])
plt.contour(r,c, G, [0])

plt.xlim(0,250)
plt.ylim(0,250)



#generate  figure 7.11B

delta_r=0.2

S1_0=[250,15]
S2_0=[250,65]
S3_0=[250,100]
S4_0=[250,150]
S5_0=[250,200]
S6_0=[250,250]
S7_0=[50,250]
S8_0=[150,250]
S9_0=[15,20]
S10_0=[15,5]


S1=odeint(dSdt_original, S1_0, times) #run simulation
S2=odeint(dSdt_original, S2_0, times)
S3=odeint(dSdt_original, S3_0, times)
S4=odeint(dSdt_original, S4_0, times)
S5=odeint(dSdt_original, S5_0, times)
S6=odeint(dSdt_original, S6_0, times)
S7=odeint(dSdt_original, S7_0, times)
S8=odeint(dSdt_original, S8_0, times)
S9=odeint(dSdt_original, S9_0, times)
S10=odeint(dSdt_original, S10_0, times)



plt.figure() 
plt.plot(S1[:,0], S1[:,1], linewidth=2)
plt.plot(S1[:,0], S1[:,1], linewidth=2)
plt.plot(S2[:,0], S2[:,1], linewidth=2)
plt.plot(S3[:,0], S3[:,1], linewidth=2)
plt.plot(S4[:,0], S4[:,1], linewidth=2)
plt.plot(S5[:,0], S5[:,1], linewidth=2)
plt.plot(S6[:,0], S6[:,1], linewidth=2)
plt.plot(S7[:,0], S7[:,1], linewidth=2)
plt.plot(S8[:,0], S8[:,1], linewidth=2)
plt.plot(S9[:,0], S9[:,1], linewidth=2)
plt.plot(S10[:,0], S10[:,1], linewidth=2)
plt.xlabel("cro concentration of (nM)")
plt.ylabel("cI oncentration of (nM)")
plt.legend()
plt.title('Figure 7.11B')


# xrange = arange(0, 2, delta)
# yrange = arange(0, 2, delta)
xrange = np.arange(0, 350, 1)
yrange = np.arange(0, 350, 1)
r, c = np.meshgrid(xrange,yrange)


F = (a + 10*a*K1*(r/2)**2)/(1+K1*(r/2)**2+ K1*K2*(r/2)**3+K3*(c/2)+K3*K4*(c/2)**2)- delta_r*(r)
G = (b + b*K3*(c/2))/(1+K1*(r/2)**2+ K1*K2*(r/2)**3+K3*(c/2)+K3*K4*(c/2)**2) - delta_c*(c)

#plot nullclines
plt.contour(r,c, F, [0])
plt.contour(r,c, G, [0])

plt.xlim(0,250)
plt.ylim(0,250)