Solve ODE with a time dependent parameter

5 views (last 30 days)
I am working through Comp Neruo Book in my spare time.
I have been manually coding solutions for the ODE's as it has provided me personally with more flexibility, but I am trying to learn how the numerical ODE suite works.
If I have a time varing parameter (Iapp) that changes when t is divisible by 0.1 seconds, how would that be modelled. I would also like to return the paramater as an array for when it. I don't see how that would would be implemented with the "Event" .
The equations below are the hodgkin-huxley model and y is
y = [V, m, h, n]
though this detail is not relevant.
The Iapp value that is returned makes no sense to me and it has no effect on the calculation anyway...
What is the correct approach to modelling this in the numerical ode suite?
tspan = 0:0.0001:0.35;
y0 = [0;0;0;0];
Iapp = 0;
[t,y,Iapp] = ode45(@(t,y) HHmodelD(t,y,Iapp),tspan,y0);
title("Membrane Potential")
function [dydt,Iapp] = HHmodelD(t,y,Iapp)
GL = 30e-9; GnaK = 12e-6; Gk_Max = 3.6e-6; Ena = 45e-3;
Ek = -82e-3; EL = -60e-3; Cm = 100e-12;
% Where I am modifying the Iapp Value
if(mod(t,0.1) == 0)
Iapp(end+1) = Iapp(end) + 0.22e-9;
Iapp(end+1) = Iapp(end);
% Note Iapp(end) in dydt1
dydt1 = (GL/Cm)*(EL - y(1)) +(Gk_Max/Cm)*(y(2)^3) * y(3)*(Ena - y(1)) + (Gk_Max/Cm)*(y(4)^4) *(Ek - y(1)) + Iapp(end)/Cm;
dydt2 = alphaM(y(1))*(1-y(2))-betaM(y(1))*y(2);
dydt3 = alphaH(y(1))*(1-y(3))-betaH(y(1))*y(3);
dydt4 = alphaN(y(1))*(1-y(4))-betaN(y(1))*y(4);
dydt = [dydt1;dydt2;dydt3;dydt4];
function aM = alphaM(V)
aM = ((10^5)*(-V-0.045))/(exp(100*(-V-0.045)) - 1);
function aH = alphaH(V)
aH = 70*exp(50*(-V-0.070));
function aN = alphaN(V)
aN = ((10^4)*(-V-0.060))/(exp(100*(-V-0.060)) - 1);
function bM = betaM(V)
bM = (4*10^3)*exp((-V - 0.070)/0.018);
function bH = betaH(V)
bH = (10^3)/(1+exp(100*(-V-0.040)));
function bN = betaN(V)
bN = 125*exp((-V-0.070)/0.08);

Accepted Answer

Star Strider
Star Strider on 9 Feb 2021
Using discrete values for ‘Iapp’ is likely not appropriate because of the adaptive algorithms used in the MATLAB ODE solvers.
See ODE with Time-Dependent Terms for an approach that will likely work.
Star Strider
Star Strider on 9 Feb 2021
Thank you!
The interpolation is the correct approach to using time-dependent terms in the ODE solution.
If you need to stop the ODE integration and then re-start it to avoid discontinuities, see the documentation section on ODE Event Location for an illustration of how to do this. I’m not certain if the Hodgkin-Huxley models require stopping and then re-starting the integration (there are several versions when I last surveyed them a few years ago) however if that is necessary for your model, using Events is the correct approach.

Sign in to comment.

More Answers (2)

Alan Stevens
Alan Stevens on 9 Feb 2021
If you do want to use the discrete values for Iapp, the following works, though you will have to decide if the results are reasonable:
dt = 0.0001;
tspan = 0:dt:0.35;
y0 = [0;0;0;0];
Iapp = zeros(1,numel(tspan));
for i = 2:numel(tspan)
Iapp(i) = Iapp(i-1);
if abs(mod(tspan(i),0.1))<dt/2
Iapp(i) = Iapp(i)+0.22e-9;
plot(tspan ,Iapp,'.')
[t,y] = ode45(@(t,y) HHmodelD(t,y,Iapp,dt),tspan,y0);
title("Membrane Potential")
function dydt = HHmodelD(t,y,Iapp,dt)
GL = 30e-9; GnaK = 12e-6; Gk_Max = 3.6e-6; Ena = 45e-3;
Ek = -82e-3; EL = -60e-3; Cm = 100e-12;
% Note Iapp in dydt1
ix = round(t/dt+1); % pointer into Iapp for use below
aM = ((10^5)*(-y(1)-0.045))/(exp(100*(-y(1)-0.045)) - 1);
aH = 70*exp(50*(-y(1)-0.070));
aN = ((10^4)*(-y(1)-0.060))/(exp(100*(-y(1)-0.060)) - 1);
bM = (4*10^3)*exp((-y(1) - 0.070)/0.018);
bH = (10^3)/(1+exp(100*(-y(1)-0.040)));
bN = 125*exp((-y(1)-0.070)/0.08);
dydt1 = (GL/Cm)*(EL - y(1)) +(Gk_Max/Cm)*(y(2)^3) * y(3)*(Ena - y(1)) + (Gk_Max/Cm)*(y(4)^4) *(Ek - y(1)) + Iapp(ix)/Cm;
dydt2 = aM*(1-y(2))-bM*y(2);
dydt3 = aH*(1-y(3))-bH*y(3);
dydt4 = aN*(1-y(4))-bN*y(4);
dydt = [dydt1;dydt2;dydt3;dydt4];
This results in the following for the membrane potential

Jan on 9 Feb 2021
Edited: Jan on 9 Feb 2021
Remember that Matlab's ODE integrators are designed to operate on smooth functions. The hard jumps of the parameter will confuse the stepsize estimator. If you are lucky, the integrator stops with an error message. With less luck, the integration provides a final value, which is dominated by rounding errors due to a huge number of integration steps.
If the model changes, stop the integration and restart it:
tspan = 0:0.0001:0.35;
t0 = 0;
y0 = [0;0;0;0];
Iapp = 0;
ready = false;
yAll = [];
tAll = [];
while ~ready
% This time span fro t0 to t0+0.1 including the original
% intermediate steps:
if t0 + 0.1 < tspan(end) % Final interval:
tSpanS = [t0, tspan(tspan > t0 & tspan < t0 + 0.1), t0 + 0.1];
tSpanS = [t0, tspan(tspan > t0)];
% Perform the integration for this interval:
[t, y] = ode45(@(t,y) HHmodelD(t,y,Iapp), tSpanS, y0);
% Append the solution:
tAll = cat(1, tAll, t);
yAll = cat(1, yAll, y);
% New initial values are former final values:
y0 = y(end, :);
t0 = t(end);
% Advance the parameter
Iapp = Iapp + 0.1;
% Has the final time been reached?
ready = abs(t0 - tspan(end)) < 10 * eps(t0);
  1 Comment
Cillian Hayde
Cillian Hayde on 9 Feb 2021
Thank you, I will look into this. It seems very important!

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!