Solve ODE with a time dependent parameter

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

Answer 1

[Star Strider]

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.

[Cillian Hayde]

The approach you highlighted seems to work but @Jan ha highlighted a potential issue...

I don't have time to check know but it seems to have worked and my answer is similar to @Alan Stevens. Thats not to say I have modelled this correctly but you have solved my MATLAB question!

clear
clc

tend = 0.35;
tspan = [0 tend];
y0 = [0;0;0;0];

%-----------------------------------------------
% Create Applied Current Vector
Itime = 0:0.001:tend;
Iapp = zeros(1,length(Itime));
Iapp(1) = 0;
for i = 2:1:length(Itime)
    if(mod(Itime(i),0.1) == 0)
        Iapp(i:end) = Iapp(i) + 0.22e-9;
    end
end

%-----------------------------------------------

%Solve
[t,y] = ode45(@(t,y) HHmodelD(t,y,Iapp,Itime),tspan,y0);
   



subplot(2,1,1);
plot(t,y(:,1))
title("Membrane Potential")
ylabel("Volts")
xlabel("seconds")

subplot(2,1,2); 
Iapp1 = interp1(Itime,Iapp,t,'nearest');
plot(t,Iapp1)
title("Applied Current")
ylabel("Current")
xlabel("seconds")

function [dydt] = HHmodelD(t,y,Iapp,Itime)
    
    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

    Iapp = interp1(Itime,Iapp,t,'nearest'); % Interpolate the data set (ITime,Iapp) at time t
%-------------------------------------------------------

    
    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/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

[Star Strider]

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.

Answer 2

[Alan Stevens]

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

[Cillian Hayde]

Thank you for yout help

Answer 3

[Jan]

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

[Cillian Hayde]

Thank you, I will look into this. It seems very important!