function [pmax1,pmax2,pmin1,pmin2,state,fd]=BGCT_subfun2(delay2,v_sr,v_p1xi,open1,open2,KK,v_p3xi,open3,open4,AA,v_es,v_p3e)


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

dt=0.00005; % Intergation step
time=13;    % Simulation time

% open1=1;    % SNr-TRN pathway, 1 open and 0 close
% open2=1;    % SNr-SRN pathway, 1 open and 0 close
%KK=1;       % Scale factor

% 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);

%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 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);
end


%The following codes are used for bifurcation analysis
TD=4; % starting time used for analysis 
xcut=xxs(1,TD/dt+1:end);
peaksmax=findpeaks(xcut);
peaksmax=sort(peaksmax);
peaksmin=findpeaks(-xcut);
peaksmin=sort(-peaksmin);

if (length(peaksmax)>=1)
    pmax1=peaksmax(end);
    pmax2=peaksmax(1);
    index=find(pmax1-peaksmax>=0.1);
    if (length(index)>=1)
        pmax2=peaksmax(index(end));
    end
else
    pmax1=xcut(1);
    pmax2=xcut(2);
end


if (length(peaksmin)>=1)
    pmin1=peaksmin(end);
    pmin2=peaksmin(1);
else
    pmin1=xcut(1);
    pmin2=xcut(2);
end


%The following codes are used for frequency analysis
Fs=1/dt;
N_n=length(xxs(1,3/dt+1:time/dt));

mean_xx=xxs(1,3/dt+1:time/dt)-mean(xxs(1,3/dt+1:time/dt));
X1=fft(mean_xx);

f1=(0:length(X1)/2-1)*Fs./(length(X1));
Power=X1.*conj(X1)/N_n;

index=find(Power(1:N_n/2)==max(Power(1:N_n/2)));

fd=f1(index); %computing the donamite frequency

%The following codes are used for state analysis
state=0;
if (abs(pmin1-pmin2)>=1 && fd>=2 && fd<=4)
    state=1;   % 
end
if (abs(pmin1-pmin2)>=1 && (fd<2 || fd>4))
    state=2;   %  
end
if (abs(pmin1-pmin2)<=1 && abs(pmax1-pmin2)<=1 && fd==0)
    state=3;   % low firing
end

if (pmin1>248) %straution
    state=4;
end