# Author : Joel Gutierrez
# Date   : 5/22/07
# Place  : UC Davis
# Class  : Chaos and Nonlinear Dynamics

# Project - Presentation of binocular rivalry model based on the following article.
# Dynamics of Travelling Waves in visual perception
# Hugh R. Wilson, Rondolph Blake, and Sang-Hun Lee

# Right now the system is ignoring the summations! This code will be used to make a 
# movie of the phase portrait showing the nullclines moving in time, therefore
# demonstrating the bistability of a rivaling system.

# Import plotting routines
from pylab import *
from numpy import *

'''
To convert the rendered frames into a video use the command:
ffmpeg -i figs/null_%5d.PNG video/null.mpg


For better quality output:
ffmpeg -i figs/null_%5d.PNG -b 98000 video/null.avi


To view the video with mPlayer:
mplayer video/null.mpg
mplayer video/null.avi

To delete the frames:
rm figs/null_*.PNG

To delete the video:
rm video/null.mpg
'''

# Time
time_step = 0
time_interval = 500
dt = 0.25
time = []

# Control parameter for the model:
tau_E = 20
tau_I = 11
tau_H = 300
E_S = 10
E_T = 10
I_wt = 0.45
H_wt = 0.47

T = []
I_T = []
H_T = []
S = []
I_S = []
H_S = []


#T cell
T.append(20)
I_T.append(10)
H_T.append(15)

#S cell
S.append(10)
I_S.append(8)
H_S.append(3)

# Arrays used to calculate the nullclines. 
Tarray = arange(-1,30,0.01)
Sarray = arange(-1,30,0.01)


# RKThreeD(Firing rate (T) time constant, Inhibition (I) time constant, Adaptation (H) time constant, T, I, H, dT/dt, dI/dt, dh/dt, time step size)
def RKThreeD( tau_E, tau_I, tau_H, p, E_fr, I_fr, H_fr, E_fcn, I_fcn, H_fcn, dt ):
	k1E = dt * E_fcn( tau_E, p, H_fr, E_fr )
	k1I = dt * I_fcn( tau_I,    I_fr, E_fr )
	k1H = dt * H_fcn( tau_H,    H_fr, E_fr )

	k2E = dt * E_fcn( tau_E, p, H_fr + k1E / 2.0, E_fr + k1E / 2.0 )
	k2I = dt * I_fcn( tau_I,    I_fr + k1I / 2.0, E_fr + k1I / 2.0 )
	k2H = dt * H_fcn( tau_H,    H_fr + k1H / 2.0, E_fr + k1H / 2.0 )

	k3E = dt * E_fcn( tau_E, p, H_fr + k2E / 2.0, E_fr + k2E / 2.0 )
	k3I = dt * I_fcn( tau_I,    I_fr + k2I / 2.0, E_fr + k2I / 2.0 )
	k3H = dt * H_fcn( tau_H,    H_fr + k2H / 2.0, E_fr + k2H / 2.0 )

	k4E = dt * E_fcn( tau_E, p, H_fr + k3E, E_fr + k3E )
	k4I = dt * I_fcn( tau_I,    I_fr + k3I, E_fr + k3I )
	k4H = dt * H_fcn( tau_H,    H_fr + k3H, E_fr + k3H )

	E_fr = E_fr  + ( k1E + 2.0 * k2E + 2.0 * k3E + k4E ) / 6.0
	I_fr = I_fr  + ( k1I + 2.0 * k2I + 2.0 * k3I + k4I ) / 6.0
	H_fr = H_fr  + ( k1H + 2.0 * k2H + 2.0 * k3H + k4H ) / 6.0

	return E_fr,I_fr,H_fr

# Differential equations describing the system
def T_dot(tau_T,P_Tplus,H_T,T): 
    return (1.0/tau_T) * (-T + 100*P_Tplus**2/((10+H_T)**2 + P_Tplus**2)) 

def H_Tdot(tau_H,H_T,T):
    return (1.0/tau_H) * (H_wt*T - H_T)

def I_Tdot(tau_I,I_T,T):
    return (1.0/tau_I) * (T - I_T)

def S_dot(tau_S,P_Splus,H_S,S): 
    return (1.0/tau_S) * (-S + 100*P_Splus**2/((10+H_S)**2 + P_Splus**2)) 

def H_Sdot(tau_H,H_S,S):
    return (1.0/tau_H) * (H_wt*S - H_S)

def I_Sdot(tau_I,I_S,S):
    return (1.0/tau_I) * (S - I_S)

# Main program solves the equations at each time step, which describe the state 
# of two rivaling neurons. 
for t in arange(0,time_interval,dt):
	#Store time value 
        time.append(t)

        #Solve for input
        P_T = E_T - I_wt*I_S[time_step]
        P_S = E_S - I_wt*I_T[time_step] 

 	#Ensure input is non-negative
        P_Tplus = max(0,P_T)
        P_Splus = max(0,P_S)


	#Solve the three equations governing the T cell and store values
        T_firing, T_inhib, T_adapt = RKThreeD(tau_E, tau_I, tau_H, P_Tplus, T[time_step],                                                                               I_T[time_step], H_T[time_step], T_dot, I_Tdot, H_Tdot, dt)
	T.append(T_firing)
	I_T.append(T_inhib)
	H_T.append(T_adapt)

	#Solve the three equations governing the S cell and store values
        S_firing, S_inhib, S_adapt = RKThreeD(tau_E, tau_I, tau_H, P_Splus, S[time_step],                                                                               I_S[time_step], H_S[time_step], S_dot, I_Sdot, H_Sdot, dt)
	S.append(S_firing)
	I_S.append(S_inhib)
	H_S.append(S_adapt)

	#Solve for the nullclines and store values
        Tnull = []
        Snull = []

	#Calcuate and save the 100th nullcline
	if time_step%100 == 0:
           for i in Tarray:
               p_t = max(0,E_T - I_wt*i)
               p_s = max(0,E_S - I_wt*i)
               Tnull.append(100*p_t**2/((10 + H_T[time_step])**2 + p_t**2))
               Snull.append(100*p_s**2/((10 + H_S[time_step])**2 + p_s**2))
	   clf
	   figure()	
           title('H_T = ' +str(H_T[time_step])+ ', H_S = ' +str(H_S[time_step]))
           xlabel('Snull')
           ylabel('Tnull')
           plot(Sarray,Tnull,'r',Snull,Tarray,'b')
           axis([-1.0,30,-1.0,30])
           savefig('figs/null_'+str(time_step).zfill(5)+'.PNG')

	#Increment the time index
        time_step = time_step + 1