MATLAB: Unstable ode45 solution

hodgkinhuxleyMATLABnanode45solutionunstable

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:
% initial
Vm0 = -30;
n0 = 0.2;
m0 = 0.1;
h0 = 0.1;
currentamp = 20;
currentamp2 = 2.2*currentamp;
period = 100;
%run HH model ode
dt = [0 2000];
y0 = [Vm0 n0 m0 h0];
[t,y]=ode45(@(t,y) hodgkinhuxley(t,y,[currentamp currentamp2 period]),dt,y0);   
% HH parameters
n = y(:,2);
m = y(:,3);
h = y(:,4);
gNa = 1;
gK = 1;
INa = gNa*m.^3.*h;
IK = gK*n.^4;
% figure;
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)
% parameters from ermentrout ch1
Vm = y(1,1); % mV


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; % mS/cm^2 from handout


gK = 36; % mS/cm^2 from handout
gL = 0.3; % mS/cm^2 from handout
ENa = 50; % mV
EK = -77; % mV
EL = -50; % mV from handout
Cm = 1; % microF/cm^2
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

Best Answer

I got it to work by limiting max step size.
options = odeset('AbsTol',[1e-8 1e-8 1e-8 1e-8],'RelTol',1e-4, 'MaxStep',0.1);
[t,y]=ode45(@(t,y) hodgkinhuxley(t,y,[currentamp currentamp2 period]),dt,y0,options);
Related Question