clear
close all

rand('state', sum(100*clock)); %random seed

dt=0.00005; % Intergation step
time=30;    % Simulation time
len1=62;


v_es=0.785*1.0e-3;
v_sr=1.0e-3; 
% v_es=1.799e-3;
delay2=0.06;
v_p3xi=0.101e-3;

% v_sr=1*1.0e-3;
v_p1xi=3.01e-4;
open1=1; open2=0;
KK=1;


open3=0;open4=0;
AA=1;

% Maximum firing rate
Qmax_i=250; 
Qmax_p3=200; 
Qmax_p1=250;
Qmax_p2=300;
Qmax_xi=500;
Qmax_s=250; 
Qmax_r=250; 

% Mean firing threshold
theta_i=15e-3;
theta_p3=9e-3;  
theta_p1=10e-3; 
theta_p2=9e-3;
theta_xi=1e-2; 
theta_s=15e-3; 
theta_r=15e-3;

% Coupling strength
v_sp3=AA*open3*0.031e-3; 
v_rp3=open4*0.031e-3;
% v_p1p3=0.005e-3;
v_xip3=0.032e-3;
v_p3p1=0.029e-3;
% v_p3xi=0.101e-3;
v_p3e=0.8e-3;


v_ee=1.02e-3;   
v_ei=1.781e-3; 
v_re=4.98e-5;
v_rs=4.97e-4;
v_p1p2=3.01e-5; 
%v_p1xi=3e-4; % STN-SNr
v_p2p2=0.748e-4; 
v_p2xi=4.501e-4; 
v_xip2=3.99e-5; 
% v_es=1.799e-3;
v_se=2.201e-3;
v_xie=0.099e-3;
v_sp1=open2*3.51e-5; 
v_rp1=KK*open1*3.51e-5;
delay=0.002;
delay1=0.004; 
% Other parameters
gamma_e=101; 
alpha=51; 
beta=199; 
sigma=0.0058;
v_sn_phi_n=1.99e-3; 

xdelay=zeros(1,ceil(delay/dt)+1);  %  SRN-TRN
xdelay1=zeros(1,ceil(delay1/dt)+1);  % SRN-EPN
xdelay2=zeros(1,ceil(delay2/dt)+1);  % TRN-SRN

xxs=zeros(1,time/dt);   S_p2MFR=zeros(1,time/dt);   S_xiMFR=zeros(1,time/dt); S_p3MFR=zeros(1,time/dt);
S_sMFR=zeros(1,time/dt);  S_rMFR=zeros(1,time/dt);    S_p1MFR=zeros(1,time/dt);
A =[];
B =[];
C =[];
D =[];
E =[];
F=[]; 
G=[];
H=[];
I=[]; 
%random initial condition
x=0.5*[0*rand;-1500+6000*rand;0.02*rand;-2*3.5*rand;0.04*rand;-0.7+1.5*rand;0.004+0.013*rand;-0.15+0.4*rand;0.0005+0.0035*rand;-0.12+0.22*rand;-0.001+0.0055*rand;
    -0.1+0.2*rand;-0.09+0.1*rand;-4+7*rand;0.025*rand;-0.6+2*rand];

xm=zeros(16,1);

for v_p1p3=linspace(0.00,0.10,len1)*1.0e-3 
for t=1:time/dt
    
    xdelay(1:end-1)=xdelay(2:end);
    xdelay(end)=x(13);
    xdelay1(1:end-1)=xdelay1(2:end);
    xdelay1(end)=x(13);
    xdelay2(1:end-1)=xdelay2(2:end);
    xdelay2(end)=x(15);
    
    S_i=Qmax_i/(1+exp(-pi/sqrt(3)*(x(3)-theta_i)/sigma));
    S_p3=Qmax_p3/(1+exp(-pi/sqrt(3)*(x(5)-theta_p3)/sigma));
    S_p1=Qmax_p1/(1+exp(-pi/sqrt(3)*(x(7)-theta_p1)/sigma));
    S_p2=Qmax_p2/(1+exp(-pi/sqrt(3)*(x(9)-theta_p2)/sigma));
    S_xi=Qmax_xi/(1+exp(-pi/sqrt(3)*(x(11)-theta_xi)/sigma));
    S_s=Qmax_s/(1+exp(-pi/sqrt(3)*(x(13)-theta_s)/sigma));
    S_r=Qmax_r/(1+exp(-pi/sqrt(3)*(x(15)-theta_r)/sigma));
    S_s_lag=Qmax_s/(1+exp(-pi/sqrt(3)*(xdelay(1)-theta_s)/sigma));
    S_s_lag1=Qmax_s/(1+exp(-pi/sqrt(3)*(xdelay1(1)-theta_s)/sigma));
    
    S_r_lag=Qmax_r/(1+exp(-pi/sqrt(3)*(xdelay2(1)-theta_r)/sigma));
    
    xm=x;
    
    k1=dt*xm(2);
    k2=dt*(xm(2)+k1/2);
    k3=dt*(xm(2)+k2/2);
    k4=dt*(xm(2)+k3);
    x(1)=xm(1)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*(gamma_e^2*(-xm(1)+S_i)-2*gamma_e*xm(2));
    k2=dt*(gamma_e^2*(-xm(1)+S_i)-2*gamma_e*(xm(2)+k1/2));
    k3=dt*(gamma_e^2*(-xm(1)+S_i)-2*gamma_e*(xm(2)+k2/2));
    k4=dt*(gamma_e^2*(-xm(1)+S_i)-2*gamma_e*(xm(2)+k3));
    x(2)=xm(2)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*xm(4);
    k2=dt*(xm(4)+k1/2);
    k3=dt*(xm(4)+k2/2);
    k4=dt*(xm(4)+k3);
    x(3)=xm(3)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*(alpha*beta*(-xm(3)+v_ee*xm(1)+v_es*S_s_lag1-v_ei*S_i)-(alpha+beta)*xm(4));
    k2=dt*(alpha*beta*(-xm(3)+v_ee*xm(1)+v_es*S_s_lag1-v_ei*S_i)-(alpha+beta)*(xm(4)+k1/2));
    k3=dt*(alpha*beta*(-xm(3)+v_ee*xm(1)+v_es*S_s_lag1-v_ei*S_i)-(alpha+beta)*(xm(4)+k2/2));
    k4=dt*(alpha*beta*(-xm(3)+v_ee*xm(1)+v_es*S_s_lag1-v_ei*S_i)-(alpha+beta)*(xm(4)+k3));
    x(4)=xm(4)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*xm(6);
    k2=dt*(xm(6)+k1/2);
    k3=dt*(xm(6)+k2/2);
    k4=dt*(xm(6)+k3);
    x(5)=xm(5)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*(alpha*beta*(-xm(5)+v_p3e*xm(1)-v_p3p1*S_p1+v_p3xi*S_xi)-(alpha+beta)*xm(6));
    k2=dt*(alpha*beta*(-xm(5)+v_p3e*xm(1)-v_p3p1*S_p1+v_p3xi*S_xi)-(alpha+beta)*(xm(6)+k1/2));
    k3=dt*(alpha*beta*(-xm(5)+v_p3e*xm(1)-v_p3p1*S_p1+v_p3xi*S_xi)-(alpha+beta)*(xm(6)+k2/2));
    k4=dt*(alpha*beta*(-xm(5)+v_p3e*xm(1)-v_p3p1*S_p1+v_p3xi*S_xi)-(alpha+beta)*(xm(6)+k3));
    x(6)=xm(6)+(k1+2*k2+2*k3+k4)/6;
    
    
    k1=dt*xm(8);
    k2=dt*(xm(8)+k1/2);
    k3=dt*(xm(8)+k2/2);
    k4=dt*(xm(8)+k3);
    x(7)=xm(7)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*(alpha*beta*(-xm(7)-v_p1p2*S_p2+v_p1xi*S_xi+v_p1p3*S_p3)-(alpha+beta)*xm(8));
    k2=dt*(alpha*beta*(-xm(7)-v_p1p2*S_p2+v_p1xi*S_xi+v_p1p3*S_p3)-(alpha+beta)*(xm(8)+k1/2));
    k3=dt*(alpha*beta*(-xm(7)-v_p1p2*S_p2+v_p1xi*S_xi+v_p1p3*S_p3)-(alpha+beta)*(xm(8)+k2/2));
    k4=dt*(alpha*beta*(-xm(7)-v_p1p2*S_p2+v_p1xi*S_xi+v_p1p3*S_p3)-(alpha+beta)*(xm(8)+k3));
    x(8)=xm(8)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*xm(10);
    k2=dt*(xm(10)+k1/2);
    k3=dt*(xm(10)+k2/2);
    k4=dt*(xm(10)+k3);
    x(9)=xm(9)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*(alpha*beta*(-xm(9)-v_p2p2*S_p2+v_p2xi*S_xi)-(alpha+beta)*xm(10));
    k2=dt*(alpha*beta*(-xm(9)-v_p2p2*S_p2+v_p2xi*S_xi)-(alpha+beta)*(xm(10)+k1/2));
    k3=dt*(alpha*beta*(-xm(9)-v_p2p2*S_p2+v_p2xi*S_xi)-(alpha+beta)*(xm(10)+k2/2));
    k4=dt*(alpha*beta*(-xm(9)-v_p2p2*S_p2+v_p2xi*S_xi)-(alpha+beta)*(xm(10)+k3));
    x(10)=xm(10)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*xm(12);
    k2=dt*(xm(12)+k1/2);
    k3=dt*(xm(12)+k2/2);
    k4=dt*(xm(12)+k3);
    x(11)=xm(11)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*(alpha*beta*(-xm(11)+v_xie*xm(1)-v_xip2*S_p2+v_xip3*S_p3)-(alpha+beta)*xm(12));
    k2=dt*(alpha*beta*(-xm(11)+v_xie*xm(1)-v_xip2*S_p2+v_xip3*S_p3)-(alpha+beta)*(xm(12)+k1/2));
    k3=dt*(alpha*beta*(-xm(11)+v_xie*xm(1)-v_xip2*S_p2+v_xip3*S_p3)-(alpha+beta)*(xm(12)+k2/2));
    k4=dt*(alpha*beta*(-xm(11)+v_xie*xm(1)-v_xip2*S_p2+v_xip3*S_p3)-(alpha+beta)*(xm(12)+k3));
    x(12)=xm(12)+(k1+2*k2+2*k3+k4)/6;

    k1=dt*xm(14);
    k2=dt*(xm(14)+k1/2);
    k3=dt*(xm(14)+k2/2);
    k4=dt*(xm(14)+k3);
    x(13)=xm(13)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*(alpha*beta*(-xm(13)-v_sp1*S_p1+v_se*xm(1)+v_sp3*S_p3-v_sr*S_r-v_sr*S_r_lag+v_sn_phi_n)-(alpha+beta)*xm(14));
    k2=dt*(alpha*beta*(-xm(13)-v_sp1*S_p1+v_se*xm(1)+v_sp3*S_p3-v_sr*S_r-v_sr*S_r_lag+v_sn_phi_n)-(alpha+beta)*(xm(14)+k1/2));
    k3=dt*(alpha*beta*(-xm(13)-v_sp1*S_p1+v_se*xm(1)+v_sp3*S_p3-v_sr*S_r-v_sr*S_r_lag+v_sn_phi_n)-(alpha+beta)*(xm(14)+k2/2));
    k4=dt*(alpha*beta*(-xm(13)-v_sp1*S_p1+v_se*xm(1)+v_sp3*S_p3-v_sr*S_r-v_sr*S_r_lag+v_sn_phi_n)-(alpha+beta)*(xm(14)+k3));
    x(14)=xm(14)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*xm(16);
    k2=dt*(xm(16)+k1/2);
    k3=dt*(xm(16)+k2/2);
    k4=dt*(xm(16)+k3);
    x(15)=xm(15)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*(alpha*beta*(-xm(15)+v_rs*S_s_lag+v_re*xm(1)-v_rp1*S_p1+v_rp3*S_p3)-(alpha+beta)*xm(16));
    k2=dt*(alpha*beta*(-xm(15)+v_rs*S_s_lag+v_re*xm(1)-v_rp1*S_p1+v_rp3*S_p3)-(alpha+beta)*(xm(16)+k1/2));
    k3=dt*(alpha*beta*(-xm(15)+v_rs*S_s_lag+v_re*xm(1)-v_rp1*S_p1+v_rp3*S_p3)-(alpha+beta)*(xm(16)+k2/2));
    k4=dt*(alpha*beta*(-xm(15)+v_rs*S_s_lag+v_re*xm(1)-v_rp1*S_p1+v_rp3*S_p3)-(alpha+beta)*(xm(16)+k3));
    x(16)=xm(16)+(k1+2*k2+2*k3+k4)/6;
   
    xxs(t)=xm(1);
        S_p2MFR(t)=S_p2;
    S_xiMFR(t)=S_xi;
    S_sMFR(t)=S_s;
    S_rMFR(t)=S_r;
    
    S_p1MFR(t)=S_p1;
     
      S_p3MFR(t)=S_p3;
end

 

phi_E =mean(xxs(5/dt+1:end));
phi_p2=mean(S_p2MFR(5/dt+1:end));
phi_Xi =mean(S_xiMFR(5/dt+1:end));
phi_S =mean(S_sMFR(5/dt+1:end));
phi_R =mean(S_rMFR(5/dt+1:end));

phi_p1=mean(S_p1MFR(5/dt+1:end));
phi_p3=mean(S_p3MFR(5/dt+1:end));

A=[A,phi_E];
B=[B,phi_p2];
C=[C,phi_Xi];
D=[D,phi_S];
E=[E,phi_R];
G=[G,phi_p1];
H=[H,phi_p3];

end
plot(linspace(0.00,0.10,len1),A,'.r')
hold  on 

plot(linspace(0.00,0.10,len1),D,'.k')
hold  on 
plot(linspace(0.00,0.10,len1),E,'.g')
hold  on 
plot(linspace(0.00,0.10,len1),G,'.y')