ODE45 and a variable that assumes multiple values during the timespan

I have tried in different ways to see what happens to voltage V and gating conductances m, n and h when, at time step x, current I switched from 0 to 0.1, and then at time step x + n it gets back to 0.

This code that I'm posting works: I integrate in chunks and then concatenate the values.

However, this is rather inefficient / not very flexible. Because if I wanted to create another subsection where I change the current to 0.3, then I would have to hard-code everything again. Is there a way to do what I'm doing here, but with more flexibility? E.g. the user specifies how many chunks they want and what values I assumes in those chunks, and the function does it for them. Is the "ODE with Time-Dependent Terms" example that is mentioned in the ode45 documentation a possible alternative?

function ODE_Hodgkin_Huxley (varargin)
    
    %% Initial values
    V=-60; % Initial Membrane voltage
    m1=alpham(V)/(alpham(V)+betam(V)); % Initial m-value
    n1=alphan(V)/(alphan(V)+betan(V)); % Initial n-value
    h1=alphah(V)/(alphah(V)+betah(V)); % Initial h-value
    y0=[V;m1;n1;h1];
    
    first_current = 0;
    second_current = 0.1;
    third_current = 0;
    
    %% Time
    t0 = 0;
    tSwitch1 = 10;
    tSwitch2 = 15;
    tEnd = 25;
    
    %% Matlab's ode45 function
    flag1 = true;
    flag2 = false;
    flag3 = false;
    [time1,V1] = ode45(@ODEMAT, [t0, tSwitch1], y0);
    
    %% Take last value as initial input for next iteration
    v2=V1(end,1);
    m2=V1(end,2);
    n2=V1(end,3);
    h2=V1(end,4);
    y02=[v2;m2;n2;h2]; % Values for next iteration
    
    %% Matlab's ode45 function
    flag1 = false;
    flag2 = true;
    flag3 = false;
    [time2,V2] = ode45(@ODEMAT,[tSwitch1, tSwitch2],y02);
    
    %% Take last value as initial input for next iteration
    v3=V2(end,1);
    m3=V2(end,2);
    n3=V2(end,3);
    h3=V2(end,4);
    y03=[v3;m3;n3;h3]; % Values for next iteration
    
    %% Matlab's ode45 function
    flag1 = false;
    flag2 = false;
    flag3 = true;
    [time3,V3] = ode45(@ODEMAT,[tSwitch2,tEnd],y03);
    
    %% Concatenate values from three iterations
    OD = cat(1,V1(:,1),V2(:,1),V3(:,1));
    ODm = cat(1,V1(:,2),V2(:,2),V3(:,2));
    ODn = cat(1,V1(:,3),V2(:,3),V3(:,3));
    ODh = cat(1,V1(:,4),V2(:,4),V3(:,4));
    time = cat(1,time1,time2,time3);
    
    %% Plots
    %% Voltage
    figure
    subplot(3,1,1)
    plot(time,OD);
    legend('ODE45 solver');
    xlabel('Time (ms)');
    ylabel('Voltage (mV)');
    title('Voltage Change for Hodgkin-Huxley Model');
    
    %% Current
    subplot(3,1,2)
    plot([0,10],0, [10,15],0.1, [15,25],0);
    legend('Current injected')
    xlabel('Time (ms)')
    ylabel('Ampere')
    title('Current')
    
    %% Gating variables
    subplot(3,1,3)
    plot(time,[ODm,ODn,ODh]);
    legend('ODm','ODn','ODh');
    xlabel('Time (ms)')
    ylabel('Value')
    title('Gating variables')
    
   
function [dydt] = ODEMAT(t,y)
    
    %% Constants
    ENa=55; % mv Na reversal potential
    EK=-72; % mv K reversal potential
    El=-49; % mv Leakage reversal potential
    
    %% Values of conductances
    gbarl=0.003; % mS/cm^2 Leakage conductance
    gbarNa=1.2; % mS/cm^2 Na conductance
    gbarK=0.36; % mS/cm^2 K conductancence
    if flag1
        I = first_current;
    elseif flag2
        I = second_current;
    elseif flag3
        I = third_current;
    end

Cm = 0.01; % Capacitance

% Values set to equal input values
V = y(1);
m = y(2);
n = y(3);
h = y(4);

gNa = gbarNa*m^3*h;
gK = gbarK*n^4;
gL = gbarl;

INa=gNa*(V-ENa);
IK=gK*(V-EK);
Il=gL*(V-El);

dydt = [((1/Cm)*(I-(INa+IK+Il))); % Normal case
    alpham(V)*(1-m)-betam(V)*m;
    alphan(V)*(1-n)-betan(V)*n;
    alphah(V)*(1-h)-betah(V)*h];
end

end

Thanks!

Answer 1

[Jan]

The example code in the documentation of ODE45 uses INTERP1 to calculate a parameter in the function to be calculated. The Dormand Prince Runge Kutta integrator with step size control is designed to operator on differentiable functions. This means, that documentation suggests a method, which is definitely driving a numerical method outside the specified limits. If I thave tried to submit something like this to my professor in numerical maths, he would have rejected it. What a pitty, that the person, who has written this documentation, has overseen this.

Nevertheless, ODE45 does compute something. Maybe the step sizes get tiny at the changes and this reduces the quality of the result massively.

Instead of a nested function, I'd prefer providing the parameter by an anonymous function. Then you do not need flags. You can write a function, which encapsulates the repeated calling of the integrator:

% [UNTESTED CODE] written in the forum's editor
function [T, Y] = IntegratorWithSteps(Fcn, t, p, y0)
% Fcn: Function handle
% t:   [t0, t1, ..., tn, tend]
% p:   [p1, p2, ..., pn] with p_i is a column vector of parameters
%      for each interval of t

nSteps = numel(t) - 1;
Tc = cell(1, nSteps);
Yc = cell(1, nSteps);

for k = 1:nSteps
    [TT, YY] = ode45(@(t,y) Fcn(t, y, p(:, k)), [t(k), t(k + 1)], y0);

    y0 = YY(end, :);  % Final position is initial value for next interval
    
    % Do not store the final position except for last interval, because
    % it is repeated as initial position in next interval:
    if k < nSteps
        Tc{k} = TT(1:end-1);
        Yc{k} = YY(1:end-1);
    else
        Tc{k} = TT;
        Yc{k} = YY;
    end
end
T = cat(1, Tc{:});
Y = cat(1, Yc{:});
end

Call this as:

[T, Y] = IntegratorWithSteps(ODEMAT, [0, 10, 15, 25], [0, 0.1, 0], y0)

And modify ODEMAT:

function dydt = ODEMAT(t, y, I)
...
    % Omitting:
    % if flag1
    %     I = first_current;
    % elseif flag2
    %     I = second_current;
    % elseif flag3
    %     I = third_current;
    % end
end

Personally I've written my own ODE45 integrator and includes the simulateous calculation of the sensitivity matrix for initial values and parameters. It uses an event finder also, which allows to set a "stage". Switching the stage is like your flags and allows the function to be integrated to use different parameters. If a stage is changed, the integrator restarts its step size control, such that the discontinuities have no side effects. For scientific work the sensitivity is required, because a trajectory is not a reliable result, if it is extremely sensitive to variations. This reveals instable functions like e.g. rolling beads on a Bean machine (Galton Board).