-1
$\begingroup$

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, as many as I define at the beginning, and depending on the number n of timesteps in which current I changes its value (which is known because the user defines it), I will call ode45 n times, every time using the last values of previous iteration as starting values. However, I am aware that under ODE45 for MATLAB there is a section for time-dependent terms. Someone has suggested that it is not correct because 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 driving a numerical method outside the specified limits. So is this correct? Can I keep my way of approaching the problem? Thanks!

function ODE (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];

    t(1) = 0;
    t(2) = 10;
    I(1) = 0; % Current in chunk 1

    t(3) = 15;
    I(2) = 0.1; % Current in chunk 2

    t(4) = 25;
    I(3) = 0; % Current in chunk3

    t(5) = 30;
    I(4) = 0;

    % 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

    function [def_temp,def_volt] = DE(varargin)

        gL=0.003; % mS/cm^2 Leakage conductance
        Cm = 0.01; % Capacitance
        EL=-49; % mv Leakage reversal potential


        dt = 0.01;
        clear chunk
        for chunk = 1:chunks
            temp = t(chunk):dt:t(chunk+1)-dt;
            volt = 1/gL * (-exp(-temp*(gL/Cm))*(I(chunk) + 60*gL + gL*EL) + I(chunk) + gL*EL); % Exact solution

            if chunk == 1
                def_volt = volt;
                def_temp = temp;
            else
                def_volt = [def_volt, volt];
                def_temp = [def_temp, temp];
            end


        end

    end
end

$\endgroup$
1
  • $\begingroup$ There is also odextend which automates some of those steps. $\endgroup$ – Lutz Lehmann Mar 12 at 10:11
0
$\begingroup$

Let's take it from the beginning. I was unable to run your code as it must be missing some input, but I see that your input to the differential equations is the current $I$ which you control as a user, and it looks something like this:enter image description here

You will notice that this input is discontinuous, and thus non-differentiable. Unfortunately matlab does not have any standard solver that handle discontinuities, and even though to the best of my knowledge this is not directly mentioned in MATLAB's documentation, it is something that is mentioned in this excellent article that outlines and compares the capabilities of MATLAB's ODE solver suits with those of other ODE suits such as Python, Julia, etc.

With that said, as 'annoying' as it is (and trust be I have been in the same place), the only way to tackle the problem in MATLAB is by the way that you have done. This is by splitting your discontinuous domain to its three chunks for which it is continuous; in this case from 0-15, 15-25 and 25-30 seconds. You can then run the ODE solver in a loop for each chunk respectively. So to answer one of your questions, your approach to the problem seems correct.

The second part of your question relates to MATLAB's ODE45 example that contains the interp1 command, and has the following code:

clear ; clc

ft = linspace(0,5,25);
f = ft.^2 - ft - 3;

gt = linspace(1,6,25);
g = 3*sin(gt-0.25);

tspan = [1 5];
ic = 1;
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[t,y] = ode45(@(t,y) myode(t,y,ft,f,gt,g), tspan, ic, opts);

fig = figure ;
fig.Color = 'w' ;

plot(t, y, 'k', ft, f, 'r', gt, g, 'b')

xlim(tspan)
grid on

xlabel('time (s)')
legend('ODE Solution', 'Input Function f(t)', 'Input Function g(t)') 

function dydt = myode(t,y,ft,f,gt,g)
f = interp1(ft,f,t); % Interpolate the data set (ft,f) at time t
g = interp1(gt,g,t); % Interpolate the data set (gt,g) at time t
dydt = -f.*y + g; % Evaluate ODE at time t
end

If you run the code above, you will notice that the input functions $f(t)$ and $g(t)$ which are linearly interpolated within the ODE solver, are both continuous and hence differentiatble. As a result, the ODE performs performs within its specified limits, so no problem in MATLAB's example.

I hope this answered your question :)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.