How to use interpolation instead of integrating in chunks with ODE45?

I have written this code to integrate voltage V and gating conductances m, n and h in 30ms with injection of current from ms 10 to 20. The code works well and it shows the behaviour I want to see. I know it's suboptimal, but just out of curiosity I'm trying to implement the same code with interpolation instead of integration in chunks. However, I'm failing to do so.

This is the integration with interpolation, which is giving the incorrect result:

    myode = @(T,y0) [((1/Cm)*interp1(t, I, T, 'linear','extrap')-(INa+IK+Il)); % Normal case
        alpham(y0(1,1))*(1-y0(2,1))-betam(y0(1,1))*y0(2,1);
        alphan(y0(1,1))*(1-y0(3,1))-betan(y0(1,1))*y0(3,1);
        alphah(y0(1,1))*(1-y0(4,1))-betah(y0(1,1))*y0(4,1)];
    
    [time,V] = ode45(myode, t, y0, I);

This is the correct one, in chunks:

    [time,V] = ode45(@ODEMAT, [t(chunk), t(chunk+1)], y);   
    
    %%% Calls ODEMAT, in which theres' the integration (extensive code posted below):
    
    dydt = [((1/Cm)*(I(chunk)-(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]; 

This is the code with the integration in chunks:

function Z2_chunks ()
    
    %% 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];
    
    t = [1:30];
    I = [zeros(1,10),ones(1,10),zeros(1,10)];
    
    % Plotting purposes (set I(idx) equal to last value of I)
    idx = numel(t);
    I(idx) = 0.1;
    
    chunks = numel(t) - 1;
    
    for chunk = 1:chunks
        
        if chunk == 1
            V=-60; % Initial Membrane voltage
            m=alpham(V)/(alpham(V)+betam(V)); % Initial m-value
            n=alphan(V)/(alphan(V)+betan(V)); % Initial n-value
            h=alphah(V)/(alphah(V)+betah(V)); % Initial h-value
            y=[V;m;n;h];
        else
            y = V(end, :);  % Final position is initial value for next interval
        end

        [time,V] = ode45(@ODEMAT, [t(chunk), t(chunk+1)], y);
        
        if chunk == 1
            def_time = time;
            def_v = V;
        else
            def_time = [def_time; time];
            def_v = [def_v; V];
        end       
    end
    
    OD = def_v(:,1);
    ODm = def_v(:,2);
    ODn = def_v(:,3);
    ODh = def_v(:,4);
    time = def_time;
    
    %% 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)
    stairs(t,I)
    ylim([0 5*max(I)])
    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
    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(chunk)-(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

This is the code with the integration with interpolation. As you can see, the plots are really different:

function Z1_interpol ()
    
    %% 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];
    
    t = [1:30];
    I = [zeros(1,10),ones(1,10),zeros(1,10)];
    
    % Plotting purposes (set I(idx) equal to last value of I)
    idx = numel(t);
    I(idx) = 0.1;   
    
    V=-60; % Initial Membrane voltage
    m=alpham(V)/(alpham(V)+betam(V)); % Initial m-value
    n=alphan(V)/(alphan(V)+betan(V)); % Initial n-value
    h=alphah(V)/(alphah(V)+betah(V)); % Initial h-value
    y=[V;m;n;h];
    
    %% 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
    Cm = 0.01; % Capacitance
    
    %% Initial values
    V=-60; % Initial Membrane voltage
    m=alpham(V)/(alpham(V)+betam(V)); % Initial m-value
    n=alphan(V)/(alphan(V)+betan(V)); % Initial n-value
    h=alphah(V)/(alphah(V)+betah(V)); % Initial h-value
    y0=[V;m;n;h];
    
    gNa = gbarNa*m^3*h;
    gK = gbarK*n^4;
    gL = gbarl;
    
    INa=gNa*(V-ENa);
    IK=gK*(V-EK);
    Il=gL*(V-El);
    
    myode = @(T,y0) [((1/Cm)*interp1(t, I, T, 'linear','extrap')-(INa+IK+Il)); % Normal case
        alpham(y0(1,1))*(1-y0(2,1))-betam(y0(1,1))*y0(2,1);
        alphan(y0(1,1))*(1-y0(3,1))-betan(y0(1,1))*y0(3,1);
        alphah(y0(1,1))*(1-y0(4,1))-betah(y0(1,1))*y0(4,1)];
    
    [time,V] = ode45(myode, t, y0, I);
    
    OD=V(:,1);
    ODm=V(:,2);
    ODn=V(:,3);
    ODh=V(:,4);
    time = time;
    
    
        %% 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)
    stairs(t,I)
    ylim([0 5*max(I)])
    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')       
    
end

Thanks!

[darova]

integration and interpolation are very different things. Can you explain more what you are trying to do?

[Samuele Bolotta]

Yes. I have a time-dependent variable in ODE45, which is current I. It's not constant - I want to be able to define it as a vector of values that can vary. So that means that this "I(chunk)" can be 0 in a first part of the timespan, but then 0.4 for a while, and then 0.3 etc.:

    dydt = [((1/Cm)*(I(chunk)-(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]; 
    

Is it clearer? Thanks!

[darova]

I still don't understand what is connection between interpolatin and integration. I see you are using interpolation inside ode45. What is the problem

[Samuele Bolotta]

I'm referring to the example in the "ODE with Time-Dependent Terms" section in the ODE45 documentation (https://it.mathworks.com/help/matlab/ref/ode45.html). The problem is that the way I'm using interpolation is definitely wrong, because I'm not getting the same result as the one I obtain with the integration in chunks (which works).

Answer 1

[Bjorn Gustavsson]

Too long code. But this is what I use for cases where I have to integrate ODEs with time-varying data for one parameter:

function drdt = myODE(t,r,data,t4data)

g = -9.82;
D = interp1(t4data,data,t,'pchip');
dydt = [r(3);
       -D*r(3)*abs(r(3));
        r(4);
        g-D*r(4)*abs(r(4))];
end

Perhaps this is the type of solution you're looking for.

HTH