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

%file symmetric_network.m
%Model of symmetric network from Figure 4.6. This code generates Figures
%4.7, 4.8, 4.9, and 4.19A


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


#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#%%%% Asymmetric parametrization    
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


#assign parameter values

k1=20
k2=20
k3=5
k4=5
n1=1
n2=4


#set time_grid for simulation
t_min=0; t_max=4; dt=0.001
times=np.arange(t_min, t_max+dt, dt) #generate time-grid list

#declare right-hand-side for model
def dSdt_original(S,t):
    dS=[0,0] #generate a list to store derivatives
    dS[0]=k1/(1+np.power(S[1],n2)) - k3*S[0]
    dS[1]=k2/(1+np.power(S[0],n1)) - k4*S[1]
    return dS


#set initial conditions State vector is [s1, s2]
# two trajectories to plot
S1_0=[3,1];
S1=odeint(dSdt_original, S1_0, times) #run simulation
S2_0=[1,3];
S2=odeint(dSdt_original, S2_0, times) #run simulation

#generate fig 4.7A
fig, (ax1, ax2) = plt.subplots(2, 1)
fig.suptitle('Figure 4.7A')

ax1.plot(times, S1[:,0], label="s1", linewidth=2)
ax1.plot(times, S1[:,1], label="s2", linewidth=2)
ax1.legend()
ax1.set_ylabel("concentration")
ax1.set_xlabel("time")
ax2.plot(times, S2[:,0], label="s1", linewidth=2)
ax2.plot(times, S2[:,1], label="s2", linewidth=2)
ax2.legend()
ax2.set_xlabel("time")
ax2.set_ylabel("concentration")



#generate  figure 4.7B

S1_0=[0,0.5]
S2_0=[0.5,0]
S3_0=[1.5,0]
S4_0=[0,1.5]
S5_0=[2.5,0]
S6_0=[3.5,0]
S7_0=[4.5,0.5]
S8_0=[0.5,4.5]
S9_0=[4.5,1.5]
S10_0=[1.5,4.5]
S11_0=[4.5,2.5]
S12_0=[2.5,4.5]
S13_0=[4.5,3.5]
S14_0=[3.5,4.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)
S11=odeint(dSdt_original, S11_0, times)
S12=odeint(dSdt_original, S12_0, times)
S13=odeint(dSdt_original, S13_0, times)
S14=odeint(dSdt_original, S14_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.plot(S11[:,0], S11[:,1], linewidth=2)
plt.plot(S12[:,0], S12[:,1], linewidth=2)
plt.plot(S13[:,0], S13[:,1], linewidth=2)
plt.plot(S14[:,0], S14[:,1], linewidth=2)
plt.xlabel("concentration of $S_1$")
plt.ylabel("concentration of $S_2$")
plt.legend()
plt.title('Figure 4.7B')


#generate nullclines
#s1 nullcline
ns12=np.arange(0, 5.1, 0.1)
ns11=k1/(k3*(1+np.power(ns12,n2)))
#s2 nullcline (k2/2 since divided by (1+1^0))
ns21=np.arange(0, 5.1, 0.1)
ns22=(k2/(1+np.power(ns21,n1)))/k4
#plot nullclines
plt.plot(ns11, ns12, 'k--')
plt.plot(ns21, ns22, 'k--')

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



#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#%%%% Symmetric parametrization    
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#declare parameter values
k1=20
k2=20
k3=5
k4=5
n1=4
n2=4

#set initial conditions State vector is [s1, s2]
# two trajectories to plot
S1_0=[3,1];
S1=odeint(dSdt_original, S1_0, times) #run simulation
S2_0=[1,3];
S2=odeint(dSdt_original, S2_0, times) #run simulation


#generate figure 4.8A
fig, (ax1, ax2) = plt.subplots(2, 1)
fig.suptitle('Figure 4.8A')

ax1.plot(times, S1[:,0], label="s1", linewidth=2)
ax1.plot(times, S1[:,1], label="s2", linewidth=2)
ax1.legend()
ax1.set_ylabel("concentration")
ax1.set_xlabel("time")
ax2.plot(times, S2[:,0], label="s1", linewidth=2)
ax2.plot(times, S2[:,1], label="s2", linewidth=2)
ax2.legend()
ax2.set_xlabel("time")
ax2.set_ylabel("concentration")



#generate  figure 4.8B

S1_0=[0,0]
S2_0=[4.5,4.5]
S3_0=[0, 0.1]
S4_0=[0.1,0]
S5_0=[0.5,0]
S6_0=[0,0.5]
S7_0=[1.5,0]
S8_0=[0,1.5]
S9_0=[4.5,0.5]
S10_0=[0.5,4.5]
S11_0=[4.5,1.5]
S12_0=[1.5,4.5]
S13_0=[4.5,3.5]
S14_0=[3.5,4.5]
S15_0=[4.5,2.5]
S16_0=[2.5,4.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)
S11=odeint(dSdt_original, S11_0, times)
S12=odeint(dSdt_original, S12_0, times)
S13=odeint(dSdt_original, S13_0, times)
S14=odeint(dSdt_original, S14_0, times)
S15=odeint(dSdt_original, S13_0, times)
S16=odeint(dSdt_original, S14_0, times)


plt.figure() 
plt.plot(S1[:,0], S1[:,1], 'k:', linewidth=2)
plt.plot(S2[:,0], S2[:,1], 'k:', 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.plot(S11[:,0], S11[:,1], linewidth=2)
plt.plot(S12[:,0], S12[:,1], linewidth=2)
plt.plot(S13[:,0], S13[:,1], linewidth=2)
plt.plot(S14[:,0], S14[:,1], linewidth=2)
plt.plot(S15[:,0], S13[:,1], linewidth=2)
plt.plot(S16[:,0], S14[:,1], linewidth=2)
plt.xlabel("concentration of $S_1$")
plt.ylabel("concentration of $S_2$")
plt.legend()
plt.title('Figure 4.8B')


#generate nullclines
#s1 nullcline
ns12=np.arange(0, 5.1, 0.1)
ns11=k1/(k3*(1+np.power(ns12,n2)))
#s2 nullcline (k2/2 since divided by (1+1^0))
ns21=np.arange(0, 5.1, 0.1)
ns22=(k2/(1+np.power(ns21,n1)))/k4
#plot nullclines
plt.plot(ns11, ns12, 'k--')
plt.plot(ns21, ns22, 'k--')

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









###generate figure 4.9#######


S1_0=[0,0]
S2_0=[4,4]
S3_0=[0.525,0.5]
S4_0=[0.5,0.525]
S5_0=[4,3.9]
S6_0=[3.9,4]
S7_0=[0.75,0.5]
S8_0=[0.5,0.75]
S9_0=[1.5,2]
S10_0=[2,1.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], 'k:', linewidth=2)
plt.plot(S2[:,0], S2[:,1], 'k:', 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("concentration of $S_1$")
plt.ylabel("concentration of $S_2$")
plt.legend()
plt.title('Figure 4.9')


#generate nullclines
#s1 nullcline
ns12=np.arange(0, 5.1, 0.1)
ns11=k1/(k3*(1+np.power(ns12,n2)))
#s2 nullcline (k2/2 since divided by (1+1^0))
ns21=np.arange(0, 5.1, 0.1)
ns22=(k2/(1+np.power(ns21,n1)))/k4
#plot nullclines
plt.plot(ns11, ns12, 'k--')
plt.plot(ns21, ns22, 'k--')

plt.xlim(0.5, 2)
plt.ylim(0.5, 2)







##generate figure 4.19A

k1=20
k2=20
k3=5
k4=5
n1=2
n2=2



#generate fig 4.19A


fig, axs = plt.subplots(2, 2)
fig.suptitle('Figure 4.19A')


k1=10;
#s1 nullcline
ns12=np.arange(0, 8.1, 0.1)
ns11=k1/(k3*(1+np.power(ns12,n2)))
#s2 nullcline 
ns21=np.arange(0, 8.1, 0.1)
ns22=(k2/(1+np.power(ns21,n1)))/k4
#setup axes
axs[0,0].plot(ns11, ns12, linewidth=1)
axs[0,0].plot(ns21, ns22, linewidth=1)
axs[0,0].set_ylabel("$S_1$ concentration")
axs[0,0].set_xlabel("$S_2$ concentration")

axs[0,0].set_xlim(0,4)
axs[0,0].set_ylim(0,4)


k1=16.2;
#s1 nullcline
ns12=np.arange(0, 8.1, 0.1)
ns11=k1/(k3*(1+np.power(ns12,n2)))
#s2 nullcline 
ns21=np.arange(0, 8.1, 0.1)
ns22=(k2/(1+np.power(ns21,n1)))/k4
#setup axes
axs[0,1].plot(ns11, ns12, linewidth=1)
axs[0,1].plot(ns21, ns22, linewidth=1)
axs[0,1].set_ylabel("$S_1$ concentration")
axs[0,1].set_xlabel("$S_2$ concentration")

axs[0,1].set_xlim(0,4)
axs[0,1].set_ylim(0,4)


k1=20;
#s1 nullcline
ns12=np.arange(0, 8.1, 0.1)
ns11=k1/(k3*(1+np.power(ns12,n2)))
#s2 nullcline 
ns21=np.arange(0, 8.1, 0.1)
ns22=(k2/(1+np.power(ns21,n1)))/k4
#setup axes
axs[1,0].plot(ns11, ns12, linewidth=1)
axs[1,0].plot(ns21, ns22, linewidth=1)
axs[1,0].set_ylabel("$S_1$ concentration")
axs[1,0].set_xlabel("$S_2$ concentration")

axs[1,0].set_xlim(0,4)
axs[1,0].set_ylim(0,4)

k1=35;
#s1 nullcline
ns12=np.arange(0, 8.1, 0.1)
ns11=k1/(k3*(1+np.power(ns12,n2)))
#s2 nullcline 
ns21=np.arange(0, 8.1, 0.1)
ns22=(k2/(1+np.power(ns21,n1)))/k4
#setup axes
axs[1,1].plot(ns11, ns12, linewidth=1)
axs[1,1].plot(ns21, ns22, linewidth=1)
axs[1,1].set_ylabel("$S_1$ concentration")
axs[1,1].set_xlabel("$S_2$ concentration")

axs[1,1].set_xlim(0,4)
axs[1,1].set_ylim(0,4)