前言

前面的文章介绍了混沌的基础知识以及三个经典的混沌映射，本文将介绍著名的Hodgkin-Huxley模型。
（传送门：基于混沌系统的文本加密算法研究（一）——混沌及混沌加密的基础知识）
（传送门：基于混沌系统的文本加密算法研究（二）——经典混沌映射）

1952年英国生理学家Alan Hodgkin和Andrew Huxley在对乌贼神经轴突电生理特性研究的基础上，提出了Hodgkin-Huxley（HH）模型。本章将介绍经典的HH模型以，并对HH模型的混沌性质进行研究。

一、 Hodgkin-Huxley模型的数学形式

20世纪50年代，Hodgkin和Huxley利用电压钳位技术对乌贼巨轴突的生理特性进行了深入的研究，在此基础上于1952年总结提出了经典的HH模型。该模型随后成为了很多不同生理结构中神经细胞的模型的基础。经典HH模型由四个微分方程组成，其表达式如式(4-1)所示：
{dVdt=1CM[I−gˉNam3h(V−VNa)−gˉKn4(V−VK)−gˉl(V−Vl)]dmdt=αm(1−m)−βmmdhdt=αh(1−h)−βhhdndt=αn(1−n)−βnn(1)\left\{ \begin{array}{l}\frac{dV}{dt}=\frac{1}{C_M}[I-\bar{g}_{Na}m^3h(V-V_{Na})-\bar{g}_Kn^4(V-V_K)-\bar{g}_l(V-V_l)]\\ \frac{dm}{dt}=\alpha_m(1-m)-\beta_mm\\ \frac{dh}{dt}=\alpha_h(1-h)-\beta_hh\\ \frac{dn}{dt}=\alpha_n(1-n)-\beta_nn\\ \end{array} \right. \tag{1}⎩⎪⎪⎨⎪⎪⎧dtdV=CM1[I−gˉNam3h(V−VNa)−gˉKn4(V−VK)−gˉl(V−Vl)]dtdm=αm(1−m)−βmmdtdh=αh(1−h)−βhhdtdn=αn(1−n)−βnn(1)
其中，V为跨膜电压；m、h、n分别表示钠离子激活变量、钠离子抑制变量和钾离子激活变量；I为外刺激电流；VNa、VK、VlV_{Na}、V_{K}、V_{l}VNa、VK、Vl分别表示钠离子平衡电位、钾离子平衡电位和漏离子平衡电位；gˉNa、gˉK、gˉl\bar{g}_{Na}、\bar{g}_{K}、\bar{g}_{l}gˉNa、gˉK、gˉl分别表示钠离子通道最大电导、钾离子通道最大电导和漏离子通道最大电导；αm、βm、αh、βh、αn和βn\alpha_m、\beta_m、\alpha_h、\beta_h、\alpha_n和\beta_nαm、βm、αh、βh、αn和βn的表达式由实验数据推导得出，如式（2）所示。
αm(V)=αˉm(V−Vˉm)1−exp[−(V−Vˉm)/Kαm]\alpha_m(V)=\frac{\bar{\alpha}_m(V-\bar{V}_m)}{1-exp[-(V-\bar{V}_m)/K_{\alpha_m}]}αm(V)=1−exp[−(V−Vˉm)/Kαm]αˉm(V−Vˉm)βm(V)=βˉmexp(−V/Kβm)\beta_m(V)=\bar{\beta}_mexp(-V/K_{\beta_m})βm(V)=βˉmexp(−V/Kβm)αh(V)=αˉhexp(−V/Kαh)(2)\alpha_h(V)=\bar{\alpha}_hexp(-V/K_{\alpha_h}) \tag{2}αh(V)=αˉhexp(−V/Kαh)(2)βh(V)=βˉh1+exp[−(V−Vˉh)/Kβh]\beta_h(V)=\frac{\bar{\beta}_h}{1+exp[-(V-\bar{V}_h)/K_{\beta_h}]}βh(V)=1+exp[−(V−Vˉh)/Kβh]βˉhαn(V)=αˉn(V−Vˉn)1−exp[−(V−Vˉn)/Kαn]\alpha_n(V)=\frac{\bar{\alpha}_n(V-\bar{V}_n)}{1-exp[-(V-\bar{V}_n)/K_{\alpha_n}]}αn(V)=1−exp[−(V−Vˉn)/Kαn]αˉn(V−Vˉn)βn(V)=βˉnexp(−V/Kβn)\beta_n(V)=\bar{\beta}_nexp(-V/K_{\beta_n})βn(V)=βˉnexp(−V/Kβn)公式（1）和（2）是温度为6.3∘6.3^{\circ}6.3∘C时的经典HH模型的表达式，其中各参数的取值如下表所示。

2、Hodgkin-Huxley模型的混沌分析

对HH模型的混沌分析可以从两个方面进行研究。一方面，可以研究外部混沌信号对HH模型的影响，即将混沌系统生成的混沌序列作为电流输入HH模型中；另一方面，可以通过引入周期性电流刺激，对HH模型的混沌性质进行研究。

（1）外部电流为混沌序列

首先对外部电流为混沌序列的情形进行研究。

（2）外部电流为周期性刺激电流

下面对引入周期性电流刺激的情形进行研究。

总结

本文介绍了Hodgkin-Huxley模型，对HH模型的的混沌性质进行了研究。通过将由Lorenz系统生成的混沌信号作为外部刺激电流，对Lorenz系统控制参数与神经纤维的放电模式的关系进行分析。通过实验仿真可知，随着控制参数\alpha的增加，神经纤维的放电模式经历了从静息状态到混沌放电的转化过程；且当Lorenz系统处于混沌状态时，以Lorenz系统生成的混沌信号作为刺激电流，HH模型也将进入混沌状态。

代码

1、外部电流为混沌序列

%% ode45
clear
clc
sigma = 10; alpha = 92;  beta= 8/3;
K0 = [0.00002, 0.05293, 0.31768, 0.59612, 1, 1, 1];
[T1,Y1] = ode45(@(t,x)HHfun(t,x,sigma,alpha,beta),[0,5000],K0);
plot(T1,Y1(:,1),'k')
xlim([4900,5000])
xlabel('time/ms')
ylabel('V')
figure
plot(T1,Y1(:,5),'k')
xlim([4900,5000])
xlabel('time/ms')
ylabel('x')%% 极值点图(ode45)
clear
sigma = 10; alpha = 0:0.1:100;  beta= 8/3;
K0 = [0.00002, 0.05293, 0.31768, 0.59612, 1, 1, 1];
point = [];
for j = 1:length(alpha)[T,Y1] = ode45(@(t,x)HHfun(t,x,sigma,alpha(j),beta),[0,2000],K0);V = Y1(T>1900,1);maxmin1 = findpeaks(V);maxpoint = max(V);maxmin2 = findpeaks(maxpoint - V);for k = 1:length(maxmin1)point = [point maxmin1(k)*i+alpha(j)];endfor k = 1:length(maxmin2)point = [point (maxpoint - maxmin2(k))*i+alpha(j)];end
end
plot(point,'k.')
xlabel('\alpha')
ylabel('极值点/mV')%% 最大Lyapunov指数(ode45)
clear
clc
sigma = 10; alpha = 20:0.01:100;  beta= 8/3;
d0 = 10^(-10);
N0 = 2000;
N = 1000;
Z = [];
for k = 1:length(alpha)Lya = 0;K0 = [0.00002, 0.05293, 0.31768, 0.59612, 1, 1, 1];K1 = [0.00002, 0.05293+10^(-10), 0.31768, 0.59612, 1 , 1, 1];for j = 1:N0+N[T1,Y1] = ode45(@(t,x)HHfun(t,x,sigma,alpha(k),beta),[0,1],K0);[T2,Y2] = ode45(@(t,x)HHfun(t,x,sigma,alpha(k),beta),[0,1],K1);d1 = sqrt((Y1(end,1) - Y2(end,1))^2 + (Y1(end,2) - Y2(end,2))^2 +(Y1(end,3) - Y2(end,3))^2 +(Y1(end,4) - Y2(end,4))^2 +...(Y1(end,5) - Y2(end,5))^2 +(Y1(end,6) - Y2(end,6))^2 + (Y1(end,7) - Y2(end,7))^2);K0 = Y1(end,:);for i = 1:7K1(i) = Y1(end,i) + (d0/d1)*(Y2(end,i) - Y1(end,i));endif j > N0Lya = Lya + log(d1/d0);endendZ = [Z Lya/N];
end
plot(alpha,Z,'k')
hold on
xlabel('\alpha')
ylabel('最大Lyapunov指数')
plot(alpha,zeros(1,length(alpha)),'r-.')%%
function df = HHfun(t,x,sigma,alpha,beta)
df = zeros(7,1);
df(1) = x(5) - (120*x(2)^3*x(4)*(x(1) - 115) + 36*x(3)^4*(x(1) + 12) + 0.3*(x(1) - 10.599));
df(2) = alpham(x(1))*(1-x(2)) - betam(x(1))*x(2);
df(3) = alphan(x(1))*(1-x(3)) - betan(x(1))*x(3);
df(4) = alphah(x(1))*(1-x(4)) - betah(x(1))*x(4);
df(5) = sigma*(x(6) - x(5));
df(6) = alpha*x(5) - x(5)*x(7) - x(6);
df(7) = x(5)*x(6) - beta*x(7);
endfunction fvalue = f(x1,x2,sigma)
fvalue = sigma*(x2 - x1);
endfunction gvalue = g(x1,x2,x3,alpha)
gvalue = alpha*x1 - x1*x3 - x2;
endfunction lvalue = l(x1,x2,x3,beta)
lvalue = x1*x2 - beta*x3;
endfunction fvalue = f1(m,n,h,V,x1)
fvalue = x1 - (120*m^3*h*(V-115) + 36*n^4*(V+12) + 0.3*(V-10.599));
endfunction gvalue = g1(m,V)
gvalue = alpham(V)*(1-m) - betam(V)*m;
endfunction gvalue = g2(n,V)
gvalue = alphan(V)*(1-n) - betan(V)*n;
endfunction gvalue = g3(h,V)
gvalue = alphah(V)*(1-h) - betah(V)*h;
endfunction am = alpham(V)
am = 0.1*(V - 25)/(1 - exp((25 - V)/10));
endfunction bm = betam(V)
bm = 4*exp(-V/18);
endfunction an = alphan(V)
an = 0.01*(V - 10)/(1 - exp((10 - V)/10));
endfunction bn = betan(V)
bn = 0.125*exp(-V/80);
endfunction ah = alphah(V)
ah = 0.07*exp(-V/20);
endfunction bh = betah(V)
bh = 1/(1 + exp((30 - V)/10));
end

2、外部电流为周期性刺激电流

%% HH模型
clear
clc
A = roundn(rand*2, -2) + 4;
f = roundn(rand*20 , -1) + 120;
V0 = roundn(rand*110, -10) - 10;
m0 = roundn(rand*1, -10);
n0 = roundn(rand*0.45, -10) + 0.3;
h0 = roundn(rand*0.5, -10) + 0.1;
x0 = [V0, m0, n0, h0, 0];
[T,Y] = ode45(@(t,x)HHfun(t,x,A,f),[0,500],x0);
V = Y(:,1);
m = Y(:,2);
n = Y(:,3);
h = Y(:,4);
plot(T,V)%% HH模型（f）
clear
clc
x0 = [0.00002, 0.05293, 0.31768, 0.59612, 0];
A = 5; f = 130;
[T,Y] = ode45(@(t,x)HHfun(t,x,A,f),[0,500],x0);
V = Y(:,1);
m = Y(:,2);
n = Y(:,3);
h = Y(:,4);
% 相空间图
figure
plot(V,m) xlabel('V') ylabel('m')
figure
plot(V,n) xlabel('V') ylabel('n')
figure
plot(V,h) xlabel('V') ylabel('h')
% plot(T,V)
% I = A*(cos(2*pi*f*T/1000) + 1)/2 - 30;
% plot(T,V,'k')
% hold on
% plot(T,I,'k')
% set(gca,'ytick',[-30,0,50,100])
% set(gca,'yticklabel',[0,0,50,100])
% xlabel('time/ms')
% ylabel('V/mV')
% grid on
% xlim([-10,20])%% 极值点分布(f)
clear
clc
A = 5; f = 0:0.1:150;
x0 = [0.00002, 0.05293, 0.31768, 0.59612, 0];
point = [];
for j = 1:length(f)[T,Y] = ode45(@(t,x)HHfun(t,x,A,f(j)),[0,500],x0);V = Y(:,1);maxmin1 = findpeaks(V);maxpoint = max(V);maxmin2 = findpeaks(maxpoint - V);clearvars ifor k = 1:length(maxmin1)point = [point maxmin1(k)*i+f(j)];endfor k = 1:length(maxmin2)point = [point (maxpoint - maxmin2(k))*i+f(j)];end
end
plot(point,'k.')
set(gca,'xtick',0:10:150)
xlabel('f/Hz')
ylabel('V/mV')%% 最大Lyapunov指数图(f)
clear
clc
A = 5; f = 120:0.01:145;
t = 500;
N = 2000;
d0 = 1e-10;
Z = [];
for j = 1:length(f)Lya = 0;x0 = [0.00002, 0.05293, 0.31768, 0.59612, 0];x1 = x0;x1(1) = x1(1) + d0;for k = 1:N[T1,Y1] = ode45(@(t,x)HHfun(t,x,A,f(j)),[0,1],x0);[T2,Y2] = ode45(@(t,x)HHfun(t,x,A,f(j)),[0,1],x1);x0 = Y1(end,:); x1 = Y2(end,:);d1 = sqrt((x0(1) - x1(1))^2 + (x0(2) - x1(2))^2 + (x0(3) - x1(3))^2 + (x0(4) - x1(4))^2 + (x0(5) - x1(5))^2);x1 = x0 + (d0/d1)*(x1 - x0);if k > tLya = Lya + log(d1/d0);endendZ = [Z Lya/(k-t)];
end
plot(f,Z,'k')
hold on
plot(f,zeros(1,length(f)),'r-')
xlabel('f/Hz')
ylabel('Lyapunov')
% xlim([122,146])%% HH模型（A）
clear
clc
x0 = [0.00002, 0.05293, 0.31768, 0.59612, 0];
A = 5; f = 130;
[T,Y] = ode45(@(t,x)HHfun(t,x,A,f),[0,1000],x0);
V = Y(:,1);
V = Y(:,1);
m = Y(:,2);
n = Y(:,3);
h = Y(:,4);
% 相空间图
figure
plot(V,m)
xlabel('V')
ylabel('m')
figure
plot(V,n)
xlabel('V')
ylabel('n')
figure
plot(V,h)
xlabel('V')
ylabel('h')
%
% plot(T,V)
% I = A*(cos(2*pi*f*T/1000) + 1)/2 - 30;
% plot(T,V,'k')
% hold on
% plot(T,I,'k')
% % set(gca,'ytick',[-6,-2,-1,0,1,2,3,4])
% % set(gca,'yticklabel',[0,-2,-1,0,1,2,3,4])
% set(gca,'ytick',[-30,0,50,100])
% set(gca,'yticklabel',[0,0,50,100])
% xlabel('time/ms')
% ylabel('V/mV')
% grid on%% 极值点分布（A）
clear
clc
A = 3:0.01:7; f = 130;
x0 = [0.00002, 0.05293, 0.31768, 0.59612, 0];
point = [];
for j = 1:length(A)[T,Y] = ode45(@(t,x)HHfun(t,x,A(j),f),[0,1000],x0);V = Y(:,1);maxmin1 = findpeaks(V);maxpoint = max(V);maxmin2 = findpeaks(maxpoint - V);clearvars ifor k = 1:length(maxmin1)point = [point maxmin1(k)*i+A(j)];endfor k = 1:length(maxmin2)point = [point (maxpoint - maxmin2(k))*i+A(j)];end
end
plot(point,'k.')
xlabel('A/\muA/cm^{2}')
ylabel('V/mV')%% 最大Lyapunov指数图（A）
clear
clc
A = 4:0.001:6; f = 130;
t = 500;
N = 1000;
d0 = 1e-10;
Z = [];
for j = 1:length(A)Lya = 0;x0 = [0.00002, 0.05293, 0.31768, 0.59612, 0];x1 = x0;x1(1) = x1(1) + d0;for k = 1:N[T1,Y1] = ode45(@(t,x)HHfun(t,x,A(j),f),[0,1],x0);[T2,Y2] = ode45(@(t,x)HHfun(t,x,A(j),f),[0,1],x1);x0 = Y1(end,:); x1 = Y2(end,:);d1 = sqrt((x0(1) - x1(1))^2 + (x0(2) - x1(2))^2 + (x0(3) - x1(3))^2 + (x0(4) - x1(4))^2 + (x0(5) - x1(5))^2);x1 = x0 + (d0/d1)*(x1 - x0);if k > tLya = Lya + log(d1/d0);endendZ = [Z Lya/(k-t)];
end
plot(A,Z,'k')
hold on
plot(A,zeros(1,length(A)),'r-')
xlabel('A/\muA/cm^{2}')
ylabel('Lyapunov')
% xlim([122,146])%% 最大Lyapunov指数图（A、f）
clear
clc
A = 4:0.001:6; f = 130;
t = 500;
N = 1000;
d0 = 1e-10;
Z = [];
for j = 1:length(A)Lya = 0;x0 = [0.00002, 0.05293, 0.31768, 0.59612, 0];x1 = x0;x1(1) = x1(1) + d0;for k = 1:N[T1,Y1] = ode45(@(t,x)HHfun(t,x,A(j),f),[0,1],x0);[T2,Y2] = ode45(@(t,x)HHfun(t,x,A(j),f),[0,1],x1);x0 = Y1(end,:); x1 = Y2(end,:);d1 = sqrt((x0(1) - x1(1))^2 + (x0(2) - x1(2))^2 + (x0(3) - x1(3))^2 + (x0(4) - x1(4))^2 + (x0(5) - x1(5))^2);x1 = x0 + (d0/d1)*(x1 - x0);if k > tLya = Lya + log(d1/d0);endendZ = [Z Lya/(k-t)];
end
plot(A,Z,'k')
hold on
plot(A,zeros(1,length(A)),'r-')
xlabel('A/\muA/cm^{2}')
ylabel('Lyapunov')
% xlim([122,146])%%
function df = HHfun(t,x,A,f)
df = zeros(5,1);
df(1) = A*(cos(2*pi*f*x(5)/1000)+1)/2 -(120*x(2)^3*x(4)*(x(1)-115) + 36*x(3)^4*(x(1)+12) + 0.3*(x(1)-10.599));
df(2) = alpham(x(1))*(1-x(2)) - betam(x(1))*x(2);
df(3) = alphan(x(1))*(1-x(3)) - betan(x(1))*x(3);
df(4) = alphah(x(1))*(1-x(4)) - betah(x(1))*x(4);
df(5) = 1;
endfunction fvalue = f1(m,n,h,z,V,A,f)
fvalue = A*(cos(2*pi*f*z/1000)+1)/2 -(120*m^3*h*(V-115) + 36*n^4*(V+12) + 0.3*(V-10.599));
endfunction gvalue = g1(m,V)
gvalue = alpham(V)*(1-m) - betam(V)*m;
endfunction gvalue = g2(n,V)
gvalue = alphan(V)*(1-n) - betan(V)*n;
endfunction gvalue = g3(h,V)
gvalue = alphah(V)*(1-h) - betah(V)*h;
endfunction am = alpham(V)
am = 0.1*(V - 25)/(1 - exp((25 - V)/10));
endfunction bm = betam(V)
bm = 4*exp(-V/18);
endfunction an = alphan(V)
an = 0.01*(V - 10)/(1 - exp((10 - V)/10));
endfunction bn = betan(V)
bn = 0.125*exp(-V/80);
endfunction ah = alphah(V)
ah = 0.07*exp(-V/20);
endfunction bh = betah(V)
bh = 1/(1 + exp((30 - V)/10));
end

基于混沌系统的文本加密算法研究系列相关推荐

基于混沌系统的文本加密算法研究（二）——经典混沌映射
基于混沌系统的文本加密算法研究(二)--经典混沌映射前言一.一维Logistic混沌映射二.二维Henon混沌映射三.三维Lorenz连续混沌映射总结代码 1.Logistic映射 2.H ...
基于混沌映射的文本加密算法研究系列
基于混沌映射的文本加密算法研究(四) 前言一.传统DES密码算法二.典型的文本混沌加密算法 1.Logistic映射 2.Henon映射 3.Lorenz映射 4.Hodgkin-Huxley模型 ...
混沌图像加密matlab,基于复合混沌系统的彩色图像加密算法及Matlab实现
第27卷第3期湖南城市学院学报 (自然科学版) Vol. 27 No.3 2018年5月 Journal of Hunan City University (Natural Scie ...
基于混沌的图像置乱加密算法及matlab的实现,基于混沌的图像置乱加密算法及MATLAB的实现...
基于混沌的图像置乱加密算法及MATLAB的实现提出了一种基于混沌映射的图像置乱加密算法.借助MATLAB6.5软 (本文共3页) 阅读全文>> 数字水印(Digital Watermar ...
基于混沌的一般彩色图像加密算法的Matlab实现
目录基于混沌的一般彩色图像加密算法的Matlab实现 Matlab代码(密文无损压缩图像保存) 加解密效果基于混沌的一般彩色图像加密算法的Matlab实现说明: 混沌系统使用Loren连续混沌系 ...
【图像加密】基于混沌系统进行灰度图像加密附Matlab代码
1 简介 1 基于混沌系统的图像加密解密 Logistic混沌置乱,先不说有多复杂,其实很简单. Logistic函数是源于一个人口统计的动力学系统,其系统方程形式如下: **X(k+1) = u ...
【图像加密】基于混沌系统的图像加密解密matlab源码含GUI
Logistic混沌置乱,先不说有多复杂,其实很简单. Logistic函数是源于一个人口统计的动力学系统,其系统方程形式如下: X(k+1) = u * X(k) * [1 - X(k)],(k=0 ...
r语言 svm 大样本_r语言基于SVM模型的文本分类研究附数据代码
1 Perceptron 与 SVM 概念介绍 1.1 感知机 (Perceptron) 感知机( perceptron ) 1957 年由 Rosenblatt 提出,是神经网络与支持向量机的基础 ...
基于bert模型的文本分类研究：“Predict the Happiness”挑战
1. 前言在2018年10月,Google发布了新的语言表示模型BERT-"Bidirectional Encoder Representations from Transformers& ...

基于混沌系统的文本加密算法研究系列

基于混沌系统的文本加密算法研究（三）

前言