I am currently modeling Hodgkin-Huxley equations. I give external input to the model at specific times and see whether action potential occurs. The inputs are short square waves (0.3 ms wide). However, my ode45 gave me nan values that depended on the time interval between my inputs. For example, having 100 ms between two inputs gave me two action potentials, but having 70 ms between two inputs gave me nan for second action potential. I noticed that when nan values come up, the dy(1) value for the ode function increases in magnitude that is outside the realistic range of the model (like around 10^17 or more). This might be more of stable/unstable equation problem but nevertheless I would appreciate any help the community gives.
Here's the plots for the two situations I mentioned above:
This has 100ms interval.
This has 70ms interval.
Here's the wrapper code:
Vm0 = -30;
n0 = 0.2;
m0 = 0.1;
h0 = 0.1;
currentamp = 20;
currentamp2 = 2.2*currentamp;
period = 100;
dt = [0 2000];
y0 = [Vm0 n0 m0 h0];
[t,y]=ode45(@(t,y) hodgkinhuxley(t,y,[currentamp currentamp2 period]),dt,y0);
n = y(:,2);
m = y(:,3);
h = y(:,4);
gNa = 1;
gK = 1;
INa = gNa*m.^3.*h;
IK = gK*n.^4;
hold on;
plot(t,y(:,1));
plot(t,currentamp*(t>=100 & t<=100.3)+currentamp2*(t>=100.3+period & t<=100.3+period+0.3));
title('H-H model');
xlabel('Time (ms)');
ylabel('Voltage (mV)/Current (mA)');
xlim([80 300]);
legend('Potential','Input current');
hold off;
Here's the ode function:
function dy = hodgkinhuxley( t, y, input)
Vm = y(1,1);
n = y(2,1);
m = y(3,1);
h = y(4,1);
period = input(3);
if t>=100 && t<=100.3
inputI = input(1);
elseif t>=100.3+period && t<=100.3+period+0.3
inputI = input(2);
else
inputI = 0;
end
gNa = 120;
gK = 36;
gL = 0.3;
ENa = 50;
EK = -77;
EL = -50;
Cm = 1;
alpha_n = 0.01*(Vm+55)/(1-exp(-(Vm+55)/10));
beta_n = 0.125*exp(-(Vm+65)/80);
alpha_m = 0.1*(Vm+40)/(1-exp(-(Vm+40)/10));
beta_m = 4*exp(-(Vm+65)/18);
alpha_h = 0.07*exp(-(Vm+65)/20);
beta_h = 1/(1+exp(-(Vm+35)/10));
dVm = (1/Cm)*(-gNa*m^3*h*(Vm-ENa)...
-gK*n^4*(Vm-EK)-gL*(Vm-EL)+inputI);
dn = alpha_n*(1-n)-beta_n*n;
dm = alpha_m*(1-m)-beta_m*m;
dh = alpha_h*(1-h)-beta_h*h;
dy = [dVm; dn; dm; dh];
end