# Solve ODE with a time dependent parameter

5 views (last 30 days)

Show older comments

Hi,

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?

clear

clc

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

plot(t,y(:,1))

title("Membrane Potential")

ylabel("Volts")

xlabel("seconds")

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;

else

Iapp(end+1) = Iapp(end);

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

end

function aH = alphaH(V)

aH = 70*exp(50*(-V-0.070));

end

function aN = alphaN(V)

aN = ((10^4)*(-V-0.060))/(exp(100*(-V-0.060)) - 1);

end

function bM = betaM(V)

bM = (4*10^3)*exp((-V - 0.070)/0.018);

end

function bH = betaH(V)

bH = (10^3)/(1+exp(100*(-V-0.040)));

end

function bN = betaN(V)

bN = 125*exp((-V-0.070)/0.08);

end

end

##### 0 Comments

### Accepted Answer

Star Strider
on 9 Feb 2021

##### 2 Comments

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.

### More Answers (2)

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;

end

end

plot(tspan ,Iapp,'.')

xlabel('t'),ylabel('Iapp')

[t,y] = ode45(@(t,y) HHmodelD(t,y,Iapp,dt),tspan,y0);

figure

plot(t,y(:,1))

title("Membrane Potential")

ylabel("Volts")

xlabel("seconds")

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

end

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

else

tSpanS = [t0, tspan(tspan > t0)];

end

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

end

### See Also

### Products

### Community Treasure Hunt

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

Start Hunting!