Izhikevch模型3_izhikevich_

% ///////////////// text0801_DHH_ 电耦合 + 高斯噪声///////////////////////////// % % Initialize 超声参数 clear; clc; dt=0.01; % Initialize time vector and time step t = 0:dt:100; %Time (ms) I1_inj = 50; %神经元N1感应电流 (uA/cm^2) 56 I2_inj = 30; %神经元N2感应电流 (uA/cm^2); 72 DHH = 0.01; %电耦合度 % r1 = 3; %时滞 R = normrnd(0,1); global num_spike num_spike = 0; % Initialize m, n, and h, which represent K+ channel activation, Na+ % channel activation, and Na+ channel inactivation, respectively. m1 = zeros(1, length(t)); m2 = zeros(1, length(t)); n1 = zeros(1, length(t)); n2 = zeros(1, length(t)); h1 = zeros(1, length(t)); h2 = zeros(1, length(t)); % Initialize alphas and betas, the rate constants alpha_m1 = zeros(1,length(t)); beta_m1 = zeros(1,length(t)); alpha_n1 = zeros(1,length(t)); beta_n1 = zeros(1,length(t)); alpha_h1 = zeros(1,length(t)); beta_h1 = zeros(1,length(t)); alpha_m2 = zeros(1,length(t)); beta_m2 = zeros(1,length(t)); alpha_n2 = zeros(1,length(t)); beta_n2 = zeros(1,length(t)); alpha_h2 = zeros(1,length(t)); beta_h2 = zeros(1,length(t)); % Initialize membrane voltage and conductance vectors V1_m = zeros(1,length(t)); V2_m = zeros(1,length(t)); g_K = zeros(1,length(t)); g_Na = zeros(1, length(t)); % Constants g_K_max = 36; %K+ channel max conductance (mS/cm^2) g_Na_max = 120; %Na+ channel max conductance (mS/cm^2) g_L = 0.3; %Leakage current conductance (mS/cm^2) E_K = -12; %K+ channel Nernst potential (mV) E_Na = 115; %Na+ channel Nernst potential (mV) E_L = 10.6; %Leakage channel Nernst potential (mV) V1_rest = -70; %Resting voltage (mV) V2_rest = -70; %Resting voltage (mV) C1_m = 1.0; %Membrance capacitance (uF/cm^2) C2_m = 1.0; %Membrance capacitance (uF/cm^2) for i = 1:(length(t)-1) % V1_m(1) = -65; % V2_m(1) = -65; % Calculate all alpha and beta values alpha_m1(i) = 0.1*((25 - V1_m(i))/(exp((25-V1_m(i))/10)-1)); beta_m1(i) = 4*exp(-1*V1_m(i)/18); alpha_n1(i) = 0.01*((10-V1_m(i))/(exp((10-V1_m(i))/10)-1)); beta_n1(i) = 0.125*exp(-1*V1_m(i)/80); alpha_h1(i) = 0.07*exp(-1*V1_m(i)/20); beta_h1(i) = 1/(exp((30-V1_m(i))/10)+1); alpha_m2(i) = 0.1*((25 - V2_m(i))/(exp((25-V2_m(i))/10)-1)); beta_m2(i) = 4*exp(-1*V2_m(i)/18); alpha_n2(i) = 0.01*((10-V2_m(i))/(exp((10-V2_m(i))/10)-1)); beta_n2(i) = 0.125*exp(-1*V2_m(i)/80); alpha_h2(i) = 0.07*exp(-1*V2_m(i)/20); beta_h2(i) = 1/(exp((30-V2_m(i))/10)+1); % Initial conditions if(i == 1) m1(i) = alpha_m1(i)/(alpha_m1(i) + beta_m1(i)); n1(i) = alpha_n1(i)/(alpha_n1(i) + beta_n1(i)); h1(i) = alpha_h1(i)/(alpha_h1(i) + beta_h1(i)); m2(i) = alpha_m2(i)/(alpha_m2(i) + beta_m2(i)); n2(i) = alpha_n2(i)/(alpha_n2(i) + beta_n2(i)); h2(i) = alpha_h2(i)/(alpha_h2(i) + beta_h2(i)); end % Calculate conductances g_Na(i) = m1(i)^3*g_Na_max*h1(i); g_K(i) = n1(i)^4*g_K_max; % Calculate currents I1_Na = g_Na(i)*(V1_m(i) - E_Na); I1_K = g_K(i)*(V1_m(i) - E_K); I1_L = g_L*(V1_m(i) - E_L); I1_ion1 = I1_inj - I1_K - I1_Na - I1_L + DHH*(V2_m(i)- V1_m(i)); I1_ion = I1_ion1 + normrnd(0,1); %高斯噪声 I2_Na = g_Na(i)*(V2_m(i) - E_Na); I2_K = g_K(i)*(V2_m(i) - E_K); I2_L = g_L*(V2_m(i) - E_L); I2_ion1 = I2_inj - I2_K - I2_Na - I2_L + DHH*(V1_m(i)-V2_m(i)); I2_ion = I2_ion1 + normrnd(0,0.01); %高斯噪声 % Calculate membrane voltage, m, n, and h using Euler's method V1_m(i+1) = V1_m(i) + I1_ion/C1_m*dt; m1(i+1) = m1(i) + (alpha_m1(i)*(1-m1(i))-beta_m1(i)*m1(i))*dt; n1(i+1) = n1(i) + (alpha_n1(i)*(1-n1(i))-beta_n1(i)*n1(i))*dt; h1(i+1) = h1(i) + (alpha_h1(i)*(1-h1(i))-beta_h1(i)*h1(i))*dt; V2_m(i+1) = V2_m(i) + I2_ion/C2_m*dt; m2(i+1) = m2(i) + (alpha_m2(i)*(1-m2(i))-beta_m2(i)*m2(i))*dt; n2(i+1) = n2(i) + (alpha_n2(i)*(1-n2(i))-beta_n2(i)*n1(i))*dt; h2(i+1) = h2(i) + (alpha_h2(i)*(1-h2(i))-beta_h2(i)*h2(i))*dt; %同步差 e1=V1_m-V2_m; e2=m1-m2; e3=h1-h2; e4=n1-n2; e5=1/4*(abs(e1)+abs(e2)+abs(e3)+abs(e4)); end V1_m = V1_m + V1_rest; V2_m = V2_m + V2_rest; figure; % subplot(3,1,1);plot(t,V1_m,'r--',t,V2_m,'b-'); % legend({'神经元N1','神经元N2'}); % xlabel('time (ms)');ylabel('membrane voltage (mV)'); % subplot(3,1,2);plot(t,V1_m,'r--'); % subplot(3,1,3);plot(t,V2_m,'b--'); % title('V1,V2 相平面图'); plot(t,V1_m,'r--',t,V2_m,'b-'); legend({'V1:I1=50uA','V2:I2=30uA'}); xlabel('time (ms)');ylabel('membrane voltage (mV)'); figure; subplot(5,1,1);plot(t,e1 );ylabel('e1');title('DHH=0.01 τ=3时同步差'); subplot(5,1,2);plot(t,e2);ylabel('e2'); subplot(5,1,3);plot(t,e3);ylabel('e3'); subplot(5,1,4);plot(t,e4);ylabel('e4'); subplot(5,1,5);plot(t,e5);ylabel('e5');xlabel('time (ms)'); % legend({'e1','e2','e3','e4'}); figure; plot(V1_m,V2_m); % legend({'神经元N1','神经元N2'}); xlabel('V1_m');ylabel('V2_m'); title('V1,V2 相平面图'); % figure; % plot(t,e1);ylabel('e1');title('DHH=0.01 τ=3时同步差');