# 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, 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



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:

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.