first_fire_probability_finite_difference_method.m

这是基于生物学中实际神经元动作电位产生机制而建立的神经元发火模型。

clear all;
clf;
global MU D TAU0


MU=0.0009;
TAU0=0.1;
D=0.01;
D=0.02;
endt=4.0;

vr=-0.01;
vl=-0.09;

n_v=100;
n_t=100;

x = linspace(vl,vr,n_v);
t = linspace(0,endt,n_t);

delta_v=(vr-vl)./(n_v-1);
delta_t=endt./(n_t-1);

J=10;
J=20;

x(J)
% A surface plot is often a good way to study a solution.
%surf(x,t,u)    
%title('Numerical solution computed with 20 mesh points.')
%xlabel('Distance x')
%ylabel('Time t')

aa=x.^2+MU;
bb=D.*TAU0;
A=0.25*(aa)./delta_v-0.5*D.*TAU0./delta_v./delta_v;
B=1.0/delta_t+D.*TAU0./delta_v./delta_v;
C=-0.25*aa./delta_v-0.5*bb./delta_v./delta_v;
DD=0.5*bb./delta_v./delta_v-0.25*(aa)./delta_v;
E=1.0./delta_t-bb./delta_v./delta_v;
F=0.25.*(aa)./delta_v+0.5*D.*TAU0./delta_v./delta_v;

fp=zeros(n_t,n_v);
fp(:,n_v)=0.0;
fp(1,:)=1.0;
fp(1,n_v)=0.0;

coe_a=zeros(n_v-2,n_v-2);

coe_c=zeros(n_v-2,1);


for k=2:n_v-3;
       coe_a(k,k)=B;
       coe_a(k,k+1)=C(k+1); 
       coe_a(k,k-1)=A(k+1);
end

 coe_a(1,1)=A(2)+B;
 coe_a(1,2)=C(2);
 coe_a(n_v-2,n_v-2)=B;
 coe_a(n_v-2,n_v-3)=A(n_v-1);

 for j=1:n_t-1
  %for j=1:1
     
     for k=2:n_v-2;
           coe_c(k)=DD(k+1).*fp(j,k)+E.*fp(j,k+1)+F(k+1).*fp(j,k+2);
     end
    coe_c(1)=DD(2).*fp(j,1)+E.*fp(j,2)+F(2).*fp(j,3); 
    y=coe_a\coe_c;
    for m=1:n_v-2
        if y(m)<0
            y(m)=0.0;
        elseif y(m)>1.0
            y(m)=1.0;
        end
    end
    fp(j+1,2:n_v-1)=y;
    fp(j+1,1)=fp(j+1,2);
 end
 

plot(t,fp(:,J))
%title('Solution at x = J')
%xlabel('Time t')
%ylabel('firing probability')

FPJ=zeros(n_t,2);
PDFJ=zeros(n_t-1,2);
FPJ(:,1)=t';
FPJ(:,2)=1.0-fp(:,J);

FFP3D=zeros(n_t.*n_v,3);


        

y2=zeros(n_t,1);
y2(1:n_t-1)=fp(2:n_t,J);

y3=fp(:,J)-y2;

PDFJ(:,1)=t(1:n_t-1)';
PDFJ(:,2)=abs(y3(1:n_t-1))./delta_t;
figure
plot(PDFJ(:,1),PDFJ(:,2))

save 'E:\computation program\computation program\NeuronModel\QuadraticNonoinear\FP2.txt' FPJ -ascii 

save 'E:\computation program\computation program\NeuronModel\QuadraticNonoinear\PDF2.txt' PDFJ -ascii 

fi=fopen('P3D.txt','a+');

%save 'E:\computation program\computation program\NeuronModel\QuadraticNonoinear\P3D.txt' FFP3D -ascii 

for j=1:n_t;
    for k=1:n_v;
        fprintf(fi,'%5.6f  %5.6f  %5.6f\n',t(j),x(k),1.0-fp(j,k));
    end
end
status=fclose(fi);

% --------------------------------------------------------------