unstable ode45 solution

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

Answer 1

[Kenny Kim]

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);

Answer 2

[Jan]

Your function to be integrated is not differentiable. Matlab's ODE integrators are not designed to handles this. See http://www.mathworks.com/matlabcentral/answers/59582#answer_72047 . Restricting the stepsize is not a reliable solution but only to split the integration to intervals with smooth values.

[Tracey Rochester]

Hi, your answer and link were too advanced for me Jan. Can you simplify please? Or provide guidance on how to achieve it, rather than how not to? Many thanks.