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

1 Ansicht (letzte 30 Tage)
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!
  4 Kommentare
darova
darova am 21 Mär. 2021
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
Samuele Bolotta am 21 Mär. 2021
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).

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Bjorn Gustavsson
Bjorn Gustavsson am 22 Mär. 2021
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

Weitere Antworten (0)

Community Treasure Hunt

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

Start Hunting!

Translated by