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

1 次查看（过去 30 天）
Samuele Bolotta
Samuele Bolotta 2021-3-21
回答： Bjorn Gustavsson 2021-3-22
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 个评论
显示 2更早的评论
darova
darova 2021-3-21
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 2021-3-21
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).

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Bjorn Gustavsson
Bjorn Gustavsson 2021-3-22
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by