1 view (last 30 days)

Show older comments

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!

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).

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

Start Hunting!