MATLAB Answers

ODE45 and a variable that assumes multiple values during the timespan

1 view (last 30 days)
Samuele Bolotta
Samuele Bolotta on 9 Mar 2021
Edited: Jan on 9 Mar 2021
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!

Accepted Answer

Jan
Jan on 9 Mar 2021
Edited: Jan on 9 Mar 2021
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).

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by