#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
% file synchronized_repressilators.py
% model of synchronized repressilator networks
% from Garcia-Ojalvo et al. (2004) PNAS 101 pp. 10955-10960
% problem 7.8.19

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



#assign parameter values

alpha0=0
alpha=216 
n=2
beta=1
kappa=20
ks0=1
eta=2
ks1=0.01
Q=0.9


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


#assign initial condition.  
#state is [m11 p11 m21 p21 m31 p31 s1 m12 p12 m22 p22 m32 p32 s2]
Sinit=[0, 10, 0, 0, 0, 0, 0, 0, 0, 0, 10, 0, 0, 0]

 
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(14) #generate a list to store derivatives
    
    ma1=S[0]
    pa1=S[1]
    mb1=S[2]
    pb1=S[3]
    mc1=S[4]
    pc1=S[5]
    s1=S[6]
    
    ma2=S[7]
    pa2=S[8]
    mb2=S[9]
    pb2=S[10]
    mc2=S[11]
    pc2=S[12]
    s2=S[13]
    
    se=Q*((1/2)*(s1+s2))
    
    dS[0]= alpha0 + alpha/(1+pc1**n) - ma1
    dS[1]=beta*(ma1-pa1)
    dS[2]= alpha0 + alpha/(1+pa1**n) - mb1
    dS[3]=beta*(mb1-pb1)
    dS[4]= alpha0 + alpha/(1+pb1**n) + kappa*s1/(1+s1) - mc1
    dS[5]=beta*(mc1-pc1)
    dS[6]=-ks0*s1+ks1*pa1-eta*(s1-se)

    dS[7]= alpha0 + alpha/(1+pc2**n) - ma2
    dS[8]=beta*(ma2-pa2)
    dS[9]= alpha0 + alpha/(1+pa2**n) - mb2
    dS[10]=beta*(mb2-pb2)
    dS[11]= alpha0 + alpha/(1+pb2**n) + kappa*s2/(1+s2) - mc2
    dS[12]=beta*(mc2-pc2)
    dS[13]=-ks0*s2+ks1*pa2-eta*(s2-se)
    

    return dS


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

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