cpg.m

用matlab实现模拟神经元模型

clear all
%parameter of DSI
% cell passive parameter
CmD=4.17;       %uF/cm2 membrane capacity
RmD=30;         %M Ohs  membrane resistance
ElD=-40;         %mV  resting potential
% threshold parameter
TfD=-48;        %mV finite threshold
ToD=-0;         %mV initial threshold after spking
tauTD=50;       %ms time constance of threshold
% Ik parameter
gmaxkD=0.25;    %mS max conductance of k
taukD=20;       %ms time constance of k conductance
EkD=-80;        %mV equilibrium potential of c
% Ic parameter
gmaxcD=0.00223;  %mS max conductance of c
taucoD=20;       %ms onset time constance of c conductance
taucdD=5000;     %ms decay time constance of c conductance
EcD=-80;        %mV equilibrium potential of c
% parameter of synapses
gmaxRD=0.012; tauoRD=25; taudRD=15000; EsRD=10;
gmaxDD=0.001; tauoDD=100; taudDD=500; tauDD=2000; EsDD=10;
gmaxVD=0.6; tauVD=150; EsVD=-80;
gmaxCD1=0.045; tauCD=400; EsCD1=10;
gmaxCD2=0.0015; tauoCD=400; taudCD=5000; EsCD2=-80;

%parameter of VSI
% cell passive parameter
RmV=15;         %M Ohs  membrane resistance
CmV=6.13;       %uF/cm2 membrane capacity
ElV=-55;         %mV  resting potential
% threshold parameter
TfV=-43;        %mV finite threshold
ToV=-30;         %mV initial threshold after spking
tauTV=10;       %ms time constance of threshold
% Ik parameter
gmaxkV=0.08;    %mS max conductance of k
taukV=5;       %ms time constance of k conductance
EkV=-80;        %mV equilibrium potential of c
% Ic parameter
gmaxcV=0.0075;  %mS max conductance of c
taucoV=5;       %ms onset time constance of c conductance
taucdV=400;     %ms decay time constance of c conductance
EcV=-80;        %mV equilibrium potential of c
% Ia parameter
gmaxa=1.287;
Bm=40; Cm=-5.5;taum=10;
Bh=68; Ch=6.7;tauh=760;
% parameter of synapses
gmaxVV=0.000025; tauoVV=750; taudVV=3000; EsVV=10;
gmaxDV=0.003; tauoDV=80; taudDV=800; EsDV=-80;
gmaxCV=0.003; tauoCV=100; taudCV=1500; EsCV=10;

%parameter of C2
% cell passive parameter
RmC=12;         %M Ohs  membrane resistance
CmC=4.17;       %uF/cm2 membrane capacity
ElC=-50;         %mV  resting potential
% threshold parameter
TfC=-30;        %mV finite threshold
ToC=-10;         %mV initial threshold after spking
tauTC=100;       %ms time constance of threshold
% Ik parameter
gmaxkC=0.033;    %mS max conductance of k
taukC=40;       %ms time constance of k conductance
EkC=-80;        %mV equilibrium potential of c
% Ic parameter
gmaxcC=0.00083;  %mS max conductance of c
taucoC=40;       %ms onset time constance of c conductance
taucdC=10000;     %ms decay time constance of c conductance
EcC=-80;        %mV equilibrium potential of c
% parameter of synapses
gmaxDC1=0.02; tauoDC1=200; taudDC1=800; EsDC=10;
gmaxDC2=0.00001; tauoDC2=200; taudDC2=2000;
gmaxVC1=0.07; tauVC=400; EsVC=-80;
gmaxVC2=0.000008; tauoVC=400; taudVC=2000;

% number of spike
ND=2; NV=1; NC=1; NR=1;
% fire tick
tfD=[0 0]; tfV=[0]; tfC=[0]; tfR=[0];
%initial state
VD=ElD; VV=ElV; VC=ElC;
TD=ToD; TV=TfV; TC=TfC;
m=0; h=1;

% interval of time
th=1;       %ms
% constant external input
IeD=-0.2; IeV=0; IeC=0;

for t=1:th:30000
    % stimulus 1
    if(t==1000) NR=NR+1; tfR=[tfR t]; end
    if(t==1100) NR=NR+1; tfR=[tfR t]; end
    if(t==1200) NR=NR+1; tfR=[tfR t]; end
    if(t==1300) NR=NR+1; tfR=[tfR t]; end
    if(t==1400) NR=NR+1; tfR=[tfR t]; end
    if(t==1500) NR=NR+1; tfR=[tfR t]; end
    if(t==1600) NR=NR+1; tfR=[tfR t]; end
    if(t==1700) NR=NR+1; tfR=[tfR t]; end
    if(t==1800) NR=NR+1; tfR=[tfR t]; end
    if(t==1900) NR=NR+1; tfR=[tfR t]; end
  
    % stimulus 2
    %if(t==1000) IeC=4; end
    %if(t==20000) IeC=0; end
    
    %DSI
    % synapses from RAMP
    gRD=0;
    if(NR>1)
        for n=2:1:NR
            gRD=gRD+taudRD*(exp(-(t-tfR(n))/taudRD)-exp(-(t-tfR(n))/tauoRD))/(taudRD-tauoRD);
        end
    end   
    IsRD=(VD-EsRD)*gmaxRD*gRD;
    % synapses from DSI
    gDD=0;
    if(ND>1)
        for n=2:1:ND
            gDD=gDD+taudDD*(exp(-(t-tfD(n))/taudDD)-exp(-(t-tfD(n))/tauoDD))/(taudDD-tauoDD)...
                +(t-tfD(n))*exp(-(t-tfD(n))/tauDD)/tauDD;
        end
    end
    IsDD=(VD-EsDD)*gmaxDD*gDD;
    % synapses from VSI
    gVD=0;
    if(NV>1)
        for n=2:1:NV
            gVD=gVD+(t-tfV(n))*exp(-(t-tfV(n))/tauVD)/tauVD;
        end
    end   
    IsVD=(VD-EsVD)*gmaxVD*gVD;
    % synapses from C2
    gCD1=0;
    gCD2=0;
    if(NC>1)
        for n=2:1:NC
            gCD1=gCD1+(t-tfC(n))*exp(-(t-tfC(n))/tauCD)/tauCD;
            gCD2=gCD2+taudCD*(exp(-(t-tfC(n))/taudCD)-exp(-(t-tfC(n))/tauoCD))/(taudCD-tauoCD);
        end
    end
    IsCD=(VD-EsCD1)*gmaxCD1*gCD1+(VD-EsCD2)*gmaxCD2*gCD2;
    % Ik and Ic
    gkD=0;
    gcD=0;
    if(ND>1)
        for n=2:1:ND
            gkD=gkD+(t-tfD(n))*exp(-(t-tfD(n))/taukD)/taukD;
            gcD=gcD+taucdD*(exp(-(t-tfD(n))/taucdD)-exp(-(t-tfD(n))/taucoD))/(taucdD-taucoD);
        end
    end
    IkD=(VD-EkD)*gmaxkD*gkD;
    IcD=(VD-EcD)*gmaxcD*gcD;
    % voltage
    VD=VD+th*(IeD-IsRD-IsDD-IsVD-IsCD-IkD-IcD-(VD-ElD)/RmD)/CmD;
    
    %VSI
    % synapses from VSI
    gVV=0;
    if(NV>1)
        for n=2:1:NV
            gVV=gVV+taudVV*(exp(-(t-tfV(n))/taudVV)-exp(-(t-tfV(n))/tauoVV))/(taudVV-tauoVV);
        end
    end   
    IsVV=(VV-EsVV)*gmaxVV*gVV;
    % synapses from DSI
    gDV=0;
    if(ND>1)
        for n=2:1:ND
            gDV=gDV+taudDV*(exp(-(t-tfD(n))/taudDV)-exp(-(t-tfD(n))/tauoDV))/(taudDV-tauoDV);
        end
    end
    IsDV=(VV-EsDV)*gmaxDV*gDV;
    % synapses from C2
    gCV=0;
    if(NC>1)
        for n=2:1:NC
            gCV=gCV+taudCV*(exp(-(t-tfC(n))/taudCV)-exp(-(t-tfC(n))/tauoCV))/(taudCV-tauoCV);
        end
    end
    IsCV=(VV-EsCV)*gmaxCV*gCV;
    % Ik and Ic
    gkV=0;
    gcV=0;
    if(NV>1)
        for n=2:1:NV
            gkV=gkV+(t-tfV(n))*exp(-(t-tfV(n))/taukV)/taukV;
            gcV=gcV+taucdV*(exp(-(t-tfV(n))/taucdV)-exp(-(t-tfV(n))/taucoV))/(taucdV-taucoV);
        end
    end
    IkV=(VV-EkV)*gmaxkV*gkV;
    IcV=(VV-EcV)*gmaxcV*gcV;
    % Ia
    mf=1/(1+exp((VV+Bm)/Cm));
    hf=1/(1+exp((VV+Bh)/Ch));
    m=m+th*(mf-m)/taum;
    h=h+th*(hf-h)/tauh;
    Ia=(VV-ElV)*gmaxa*m*h;
    % voltage
    VV=VV+th*(IeV-IsVV-IsDV-IsCV-IkV-IcV-Ia-(VV-ElV)/RmV)/CmV;

    % C2
    % synapses from DSI
    gDC1=0;
    gDC2=0;
    if(ND>1)
        for n=2:1:ND
            gDC1=gDC1+taudDC1*(exp(-(t-tfD(n))/taudDC1)-exp(-(t-tfD(n))/tauoDC1))/(taudDC1-tauoDC1);
            gDC2=gDC2+taudDC2*(exp(-(t-tfD(n))/taudDC2)-exp(-(t-tfD(n))/tauoDC2))/(taudDC2-tauoDC2);;
        end
    end
    IsDC=(VC-EsDC)*(gmaxDC1*gDC1+gmaxDC2*gDC2);
    % synapses from VSI
    gVC1=0;
    gVC2=0;
    if(NV>1)
        for n=2:1:NV
            gVC1=gVC1+(t-tfV(n))*exp(-(t-tfV(n))/tauVC)/tauVC;
            gVC2=gVC2+taudVC*(exp(-(t-tfV(n))/taudVC)-exp(-(t-tfV(n))/tauoVC))/(taudVC-tauoVC);
        end
    end   
    IsVC=(VC-EsVC)*(gmaxVC1*gVC1+gmaxVC2*gVC2);
    % Ik and Ic
    gkC=0;
    gcC=0;
    if(NC>1)
        for n=2:1:NC
            gkC=gkC+(t-tfC(n))*exp(-(t-tfC(n))/taukC)/taukC;
            gcC=gcC+taucdC*(exp(-(t-tfC(n))/taucdC)-exp(-(t-tfC(n))/taucoC))/(taucdC-taucoC);
        end
    end
    IkC=(VC-EkC)*gmaxkC*gkC;
    IcC=(VC-EcC)*gmaxcC*gcC;
    % voltage
    VC=VC+th*(IeC-IsDC-IsVC-IkC-IcC-(VC-ElC)/RmC)/CmC;
    
    % fire or not
    % DSI
    if(ND>1)
        TD=TfD+(TD-TfD)*exp(-(t-tfD(ND))/tauTD);
    else
        TD=TD;
    end
    if( VD>=TD)
        plot_VD(t)=0;
        ND=ND+1;
        tfD=[tfD t];
        VD=ElD;
        TD=ToD;
    else
        plot_VD(t)=VD;
    end
    % VSI
    if(NV>1)
        TV=TfV+(TV-TfV)*exp(-(t-tfV(NV))/tauTV);
    else
        TV=TV;
    end
    if( VV>=TV)
        plot_VV(t)=0;
        NV=NV+1;
        tfV=[tfV t];
        VV=ElV;
        TV=ToV;
    else
        plot_VV(t)=VV;
    end
    % C2
    if(NC>1)
        TC=TfC+(TC-TfC)*exp(-(t-tfC(NC))/tauTC);
    else
        TC=TC;
    end
    if( VC>=TC)
        plot_VC(t)=0;
        NC=NC+1;
        tfC=[tfC t];
        VC=ElC;
        TC=ToC;
    else
        plot_VC(t)=VC;
    end
    
end
subplot(3,1,1);
plot(plot_VD,'r');
subplot(3,1,2)
plot(plot_VC,'g');
subplot(3,1,3);
plot(plot_VV,'b');