# This ODE code is a comparable implementation of the BGCT model described in the
# following paper:
#
#     Title:   Bidirectional Control of Absence Seizures by the Basal Ganglia: A
#              Computational Evidence (2014)
#     Authors: Mingming Chen, Daqing Guo*, Tiebin Wang, Wei Jing, Yang Xia, Peng Xu, 
#              Cheng Luo, Pedro A. Valdes-Sosa1, and Dezhong Yao*
#     Journal: PLoS Computational Biology
#     Emails:  twqylsf@gmail.com and dqguo@uestc.edu.cn
# You can use XPPAut (http://www.math.pitt.edu/~bard/xpp/xpp.html) for runing this code.
#

# Intergation step
dt=0.00005 

# Simulation time
time=13    

# SNr-TRN pathway, 1 open and 0 close
p open1=1 

# SNr-SRN pathway, 1 open and 0 close  
p open2=1 

# Scale factor   
p KK=1    

# Maximum firing rate (Table 1 A)
Qmax_i=250 
Qmax_d1=65 
Qmax_d2=65
Qmax_p1=250
Qmax_p2=300
Qmax_xi=500
Qmax_s=250
Qmax_r=250

# Mean firing threshold (Table 1 B)
theta_i=15e-3
theta_d1=19e-3 
theta_d2=19e-3
theta_p1=10e-3 
theta_p2=9e-3
theta_xi=1e-2 
theta_s=15e-3 
theta_r=15e-3

# Coupling strength (Table 1 C)
v_ee=1.0e-3   
v_ei=1.8e-3 
v_re=5.0e-5
v_rs=5.0e-4
v_d1e=1.0e-3
v_d1d1=2.0e-4
v_d1s=1.0e-4
v_d2e=7.0e-4
v_d2d2=3.0e-4
v_d2s=5e-5
v_p1d1=1.0e-4
v_p1p2=3.0e-5
v_p2d2=3.0e-4
v_p2p2=0.75e-4 
v_p2xi=4.5e-4 
v_xip2=4.0e-5 
v_es=1.8e-3
v_se=2.2e-3
v_xie=0.1e-3

# Other parameters (Table 1 D)
gamma_e=100 
alpha=50 
beta=200 
sigma=0.006
v_sn_phi_n=2.0e-3 

# TRN-SRN
p v_sr=8.0e-4 

# STN-SNr
p v_p1xi=3e-4
 
# delay parameter 
p tau=0.05
 
v_sp1=open2*3.5e-5 
v_rp1=KK*open1*3.5e-5

# random initial condition
init x[1]=0.5*0*ran;       x[2]=-1500+6000*ran;     x[3]=0.02*ran;        x[4]=-2*3.5*ran;   x[5]=0.04*ran 
init x[6]=-0.7+1.5*ran;    x[7]=0.001+0.025*ran;    x[8]=-0.4+ran;        x[9]=0.004+0.013*ran
init x[10]=-0.15+0.4*ran;  x[11]=0.0005+0.0035*ran; x[12]=-0.12+0.22*ran; x[13]=-0.001+0.0055*ran
init x[14]=-0.1+0.2*ran;   x[15]=-0.09+0.1*ran;     x[16]=-4+7*ran;       x[17]=0.025*ran;   x[18]=-0.6+2*ran

# sigmoid function for different neural populations

S_i=Qmax_i/(1+exp(-pi/sqrt(3)*(x3-theta_i)/sigma))

S_d1=Qmax_d1/(1+exp(-pi/sqrt(3)*(x5-theta_d1)/sigma))

S_d2=Qmax_d2/(1+exp(-pi/sqrt(3)*(x7-theta_d2)/sigma))

S_p1=Qmax_p1/(1+exp(-pi/sqrt(3)*(x9-theta_p1)/sigma))

S_p2=Qmax_p2/(1+exp(-pi/sqrt(3)*(x11-theta_p2)/sigma))

S_xi=Qmax_xi/(1+exp(-pi/sqrt(3)*(x13-theta_xi)/sigma))

S_s=Qmax_s/(1+exp(-pi/sqrt(3)*(x15-theta_s)/sigma))

S_r=Qmax_r/(1+exp(-pi/sqrt(3)*(x17-theta_r)/sigma))

S_r_lag=Qmax_r/(1+exp(-pi/sqrt(3)*(delay(x17,tau)-theta_r)/sigma))

# --------------------------------cerebral cortex----------------------------------------

x1'=x2
x2'=gamma_e^2*(-x1+S_i)-2*gamma_e*x2

x3'=x4
x4'=alpha*beta*(-x3+v_ee*x1+v_es*S_s-v_ei*S_i)-(alpha+beta)*x4

# ---------------------------------striatum D1--------------------------------------------

x5'=x6
x6'=alpha*beta*(-x5+v_d1e*x1-v_d1d1*S_d1+v_d1s*S_s)-(alpha+beta)*x6

# ---------------------------------striatum D2--------------------------------------------

x7'=x8
x8'=alpha*beta*(-x7+v_d2e*x1-v_d2d2*S_d2+v_d2s*S_s)-(alpha+beta)*x8

# ----------------------------------SNr/GPi-----------------------------------------------

x9'=x10
x10'=alpha*beta*(-x9-v_p1d1*S_d1-v_p1p2*S_p2+v_p1xi*S_xi)-(alpha+beta)*x10

# -----------------------------------GPe--------------------------------------------------

x11'=x12
x12'=alpha*beta*(-x11-v_p2d2*S_d2-v_p2p2*S_p2+v_p2xi*S_xi)-(alpha+beta)*x12

# ------------------------------------STN-------------------------------------------------

x13'=x14
x14'=alpha*beta*(-x13+v_xie*x1-v_xip2*S_p2)-(alpha+beta)*x14

# ------------------------------------SRN--------------------------------------------------

x15'=x16
x16'=alpha*beta*(-x15-v_sp1*S_p1+v_se*x1-v_sr*S_r-v_sr*S_r_lag+v_sn_phi_n)-(alpha+beta)*x16  

# -----------------------------------TRN--------------------------------------------------

x17'=x18
x18'=alpha*beta*(-x17-v_rp1*S_p1+v_rs*S_s+v_re*x1)-(alpha+beta)*x18


@ xlo=0, ylo=-10, xhi=13, yhi=150, yp=x1

@ total=13, dt=0.00005, bounds=10e10, maxstore=10e15, delay=1

done