Basal ganglia-corticothalamic (BGCT) network (Chen et al., 2014)

 Download zip file 
Help downloading and running models
Accession:152113
We developed a biophysical model of the basal ganglia-corticothalamic network in this work. "... We demonstrate that the typical absence seizure activities can be controlled and modulated by the direct GABAergic projections from the substantia nigra pars reticulata (SNr) to either the thalamic reticular nucleus (TRN) or the specific relay nuclei (SRN) of thalamus, through different biophysical mechanisms. ... results highlight the bidirectional functional roles of basal ganglia in controlling and modulating absence seizures, and might provide novel insights into the therapeutic treatments of this brain disorder."
References:
1 . Chen M, Guo D, Wang T, Jing W, Xia Y, Xu P, Luo C, Valdes-Sosa PA, Yao D (2014) Bidirectional control of absence seizures by the basal ganglia: a computational evidence. PLoS Comput Biol 10:e1003495 [PubMed]
2 . Chen M, Guo D, Li M, Ma T, Wu S, Ma J, Cui Y, Xia Y, Xu P, Yao D (2015) Critical Roles of the Direct GABAergic Pallido-cortical Pathway in Controlling Absence Seizures. PLoS Comput Biol 11:e1004539 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neural mass;
Brain Region(s)/Organism: Neocortex; Thalamus; Basal ganglia;
Cell Type(s):
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: XPP; MATLAB;
Model Concept(s): Epilepsy;
Implementer(s): Guo, Daqing [dqguo at uestc.edu.cn];
function [pmax1,pmax2,pmin1,pmin2,state,fd]=BGCT_subfun2(delay,v_sr,v_p1xi,open1,open2,KK)
% This MATLAB code is a subfunction version of 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
%     
%     Parameters: "open1" and "open2" are gating parameters for the SNr-TRN
%                 pathway and SNr-SRN pathway, with 1 open and 0 close;
%                 "delay" represents GABA_B delay paramter;
%                 "v_sr" is the coupling strength of the TRN-SRN pathway
%                  v_p1xi is the coupling strength of the STN-SNr pathway

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_d1=65; 
Qmax_d2=65;
Qmax_p1=250;
Qmax_p2=300;
Qmax_xi=500;
Qmax_s=250; 
Qmax_r=250; 

% Mean firing threshold
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
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_p1xi=3e-4; % STN-SNr
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;
v_sp1=open2*3.5e-5; 
v_rp1=KK*open1*3.5e-5;

% Other parameters
gamma_e=100; 
alpha=50; 
beta=200; 
sigma=0.006;
v_sn_phi_n=2.0e-3; 

xdelay=zeros(1,ceil(delay/dt)+1);

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.001+0.025*rand;
    -0.4+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(18,1);


for t=1:time/dt
    
    xdelay(1:end-1)=xdelay(2:end);
    xdelay(end)=x(17);
    
    S_i=Qmax_i/(1+exp(-pi/sqrt(3)*(x(3)-theta_i)/sigma));
    S_d1=Qmax_d1/(1+exp(-pi/sqrt(3)*(x(5)-theta_d1)/sigma));
    S_d2=Qmax_d2/(1+exp(-pi/sqrt(3)*(x(7)-theta_d2)/sigma));
    S_p1=Qmax_p1/(1+exp(-pi/sqrt(3)*(x(9)-theta_p1)/sigma));
    S_p2=Qmax_p2/(1+exp(-pi/sqrt(3)*(x(11)-theta_p2)/sigma));
    S_xi=Qmax_xi/(1+exp(-pi/sqrt(3)*(x(13)-theta_xi)/sigma));
    S_s=Qmax_s/(1+exp(-pi/sqrt(3)*(x(15)-theta_s)/sigma));
    S_r=Qmax_r/(1+exp(-pi/sqrt(3)*(x(17)-theta_r)/sigma));
    S_r_lag=Qmax_r/(1+exp(-pi/sqrt(3)*(xdelay(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-v_ei*S_i)-(alpha+beta)*xm(4));
    k2=dt*(alpha*beta*(-xm(3)+v_ee*xm(1)+v_es*S_s-v_ei*S_i)-(alpha+beta)*(xm(4)+k1/2));
    k3=dt*(alpha*beta*(-xm(3)+v_ee*xm(1)+v_es*S_s-v_ei*S_i)-(alpha+beta)*(xm(4)+k2/2));
    k4=dt*(alpha*beta*(-xm(3)+v_ee*xm(1)+v_es*S_s-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_d1e*xm(1)-v_d1d1*S_d1+v_d1s*S_s)-(alpha+beta)*xm(6));
    k2=dt*(alpha*beta*(-xm(5)+v_d1e*xm(1)-v_d1d1*S_d1+v_d1s*S_s)-(alpha+beta)*(xm(6)+k1/2));
    k3=dt*(alpha*beta*(-xm(5)+v_d1e*xm(1)-v_d1d1*S_d1+v_d1s*S_s)-(alpha+beta)*(xm(6)+k2/2));
    k4=dt*(alpha*beta*(-xm(5)+v_d1e*xm(1)-v_d1d1*S_d1+v_d1s*S_s)-(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*(-x(7)+v_d2e*xm(1)-v_d2d2*S_d2+v_d2s*S_s)-(alpha+beta)*xm(8));
    k2=dt*(alpha*beta*(-x(7)+v_d2e*xm(1)-v_d2d2*S_d2+v_d2s*S_s)-(alpha+beta)*(xm(8)+k1/2));
    k3=dt*(alpha*beta*(-x(7)+v_d2e*xm(1)-v_d2d2*S_d2+v_d2s*S_s)-(alpha+beta)*(xm(8)+k2/2));
    k4=dt*(alpha*beta*(-x(7)+v_d2e*xm(1)-v_d2d2*S_d2+v_d2s*S_s)-(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_p1d1*S_d1-v_p1p2*S_p2+v_p1xi*S_xi)-(alpha+beta)*xm(10));
    k2=dt*(alpha*beta*(-xm(9)-v_p1d1*S_d1-v_p1p2*S_p2+v_p1xi*S_xi)-(alpha+beta)*(xm(10)+k1/2));
    k3=dt*(alpha*beta*(-xm(9)-v_p1d1*S_d1-v_p1p2*S_p2+v_p1xi*S_xi)-(alpha+beta)*(xm(10)+k2/2));
    k4=dt*(alpha*beta*(-xm(9)-v_p1d1*S_d1-v_p1p2*S_p2+v_p1xi*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_p2d2*S_d2-v_p2p2*S_p2+v_p2xi*S_xi)-(alpha+beta)*xm(12));
    k2=dt*(alpha*beta*(-xm(11)-v_p2d2*S_d2-v_p2p2*S_p2+v_p2xi*S_xi)-(alpha+beta)*(xm(12)+k1/2));
    k3=dt*(alpha*beta*(-xm(11)-v_p2d2*S_d2-v_p2p2*S_p2+v_p2xi*S_xi)-(alpha+beta)*(xm(12)+k2/2));
    k4=dt*(alpha*beta*(-xm(11)-v_p2d2*S_d2-v_p2p2*S_p2+v_p2xi*S_xi)-(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_xie*xm(1)-v_xip2*S_p2)-(alpha+beta)*xm(14));
    k2=dt*(alpha*beta*(-xm(13)+v_xie*xm(1)-v_xip2*S_p2)-(alpha+beta)*(xm(14)+k1/2));
    k3=dt*(alpha*beta*(-xm(13)+v_xie*xm(1)-v_xip2*S_p2)-(alpha+beta)*(xm(14)+k2/2));
    k4=dt*(alpha*beta*(-xm(13)+v_xie*xm(1)-v_xip2*S_p2)-(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_sp1*S_p1+v_se*xm(1)-v_sr*S_r-v_sr*S_r_lag+v_sn_phi_n)-(alpha+beta)*xm(16));
    k2=dt*(alpha*beta*(-xm(15)-v_sp1*S_p1+v_se*xm(1)-v_sr*S_r-v_sr*S_r_lag+v_sn_phi_n)-(alpha+beta)*(xm(16)+k1/2));
    k3=dt*(alpha*beta*(-xm(15)-v_sp1*S_p1+v_se*xm(1)-v_sr*S_r-v_sr*S_r_lag+v_sn_phi_n)-(alpha+beta)*(xm(16)+k2/2));
    k4=dt*(alpha*beta*(-xm(15)-v_sp1*S_p1+v_se*xm(1)-v_sr*S_r-v_sr*S_r_lag+v_sn_phi_n)-(alpha+beta)*(xm(16)+k3));
    x(16)=xm(16)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*xm(18);
    k2=dt*(xm(18)+k1/2);
    k3=dt*(xm(18)+k2/2);
    k4=dt*(xm(18)+k3);
    x(17)=xm(17)+(k1+2*k2+2*k3+k4)/6;
    
    k1=dt*(alpha*beta*(-xm(17)+v_rs*S_s+v_re*xm(1)-v_rp1*S_p1)-(alpha+beta)*xm(18));
    k2=dt*(alpha*beta*(-xm(17)+v_rs*S_s+v_re*xm(1)-v_rp1*S_p1)-(alpha+beta)*(xm(18)+k1/2));
    k3=dt*(alpha*beta*(-xm(17)+v_rs*S_s+v_re*xm(1)-v_rp1*S_p1)-(alpha+beta)*(xm(18)+k2/2));
    k4=dt*(alpha*beta*(-xm(17)+v_rs*S_s+v_re*xm(1)-v_rp1*S_p1)-(alpha+beta)*(xm(18)+k3));
    x(18)=xm(18)+(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 state analysis
state=0;
if (abs(pmin1-pmin2)>=1)
    state=1;   % SWD
elseif (abs(pmin1-pmin2)<=1 && abs(pmax1-pmin2)<=1)
    state=2;   % low firing
end

if (pmin1>248) %straution
    state=4;
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


Loading data, please wait...