Solve ODE with a time dependent parameter

10 次查看（过去 30 天）
Cillian Hayde
Cillian Hayde 2021-2-9
评论： Star Strider 2021-2-9
Hi,
I am working through Comp Neruo Book in my spare time.
I have been manually coding solutions for the ODE's as it has provided me personally with more flexibility, but I am trying to learn how the numerical ODE suite works.
If I have a time varing parameter (Iapp) that changes when t is divisible by 0.1 seconds, how would that be modelled. I would also like to return the paramater as an array for when it. I don't see how that would would be implemented with the "Event" .
The equations below are the hodgkin-huxley model and y is
y = [V, m, h, n]
though this detail is not relevant.
The Iapp value that is returned makes no sense to me and it has no effect on the calculation anyway...
What is the correct approach to modelling this in the numerical ode suite?
clear
clc
tspan = 0:0.0001:0.35;
y0 = [0;0;0;0];
Iapp = 0;
[t,y,Iapp] = ode45(@(t,y) HHmodelD(t,y,Iapp),tspan,y0);
plot(t,y(:,1))
title("Membrane Potential")
ylabel("Volts")
xlabel("seconds")
function [dydt,Iapp] = HHmodelD(t,y,Iapp)
GL = 30e-9; GnaK = 12e-6; Gk_Max = 3.6e-6; Ena = 45e-3;
Ek = -82e-3; EL = -60e-3; Cm = 100e-12;
%-------------------------------------------------------
% Where I am modifying the Iapp Value
if(mod(t,0.1) == 0)
Iapp(end+1) = Iapp(end) + 0.22e-9;
else
Iapp(end+1) = Iapp(end);
end
%-------------------------------------------------------
% Note Iapp(end) in dydt1
dydt1 = (GL/Cm)*(EL - y(1)) +(Gk_Max/Cm)*(y(2)^3) * y(3)*(Ena - y(1)) + (Gk_Max/Cm)*(y(4)^4) *(Ek - y(1)) + Iapp(end)/Cm;
dydt2 = alphaM(y(1))*(1-y(2))-betaM(y(1))*y(2);
dydt3 = alphaH(y(1))*(1-y(3))-betaH(y(1))*y(3);
dydt4 = alphaN(y(1))*(1-y(4))-betaN(y(1))*y(4);
dydt = [dydt1;dydt2;dydt3;dydt4];
function aM = alphaM(V)
aM = ((10^5)*(-V-0.045))/(exp(100*(-V-0.045)) - 1);
end
function aH = alphaH(V)
aH = 70*exp(50*(-V-0.070));
end
function aN = alphaN(V)
aN = ((10^4)*(-V-0.060))/(exp(100*(-V-0.060)) - 1);
end
function bM = betaM(V)
bM = (4*10^3)*exp((-V - 0.070)/0.018);
end
function bH = betaH(V)
bH = (10^3)/(1+exp(100*(-V-0.040)));
end
function bN = betaN(V)
bN = 125*exp((-V-0.070)/0.08);
end
end

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

采纳的回答

Star Strider
Star Strider 2021-2-9
Using discrete values for ‘Iapp’ is likely not appropriate because of the adaptive algorithms used in the MATLAB ODE solvers.
See ODE with Time-Dependent Terms for an approach that will likely work.
  2 个评论
Cillian Hayde
Cillian Hayde 2021-2-9
The approach you highlighted seems to work but @Jan ha highlighted a potential issue...
I don't have time to check know but it seems to have worked and my answer is similar to @Alan Stevens. Thats not to say I have modelled this correctly but you have solved my MATLAB question!
clear
clc
tend = 0.35;
tspan = [0 tend];
y0 = [0;0;0;0];
%-----------------------------------------------
% Create Applied Current Vector
Itime = 0:0.001:tend;
Iapp = zeros(1,length(Itime));
Iapp(1) = 0;
for i = 2:1:length(Itime)
if(mod(Itime(i),0.1) == 0)
Iapp(i:end) = Iapp(i) + 0.22e-9;
end
end
%-----------------------------------------------
%Solve
[t,y] = ode45(@(t,y) HHmodelD(t,y,Iapp,Itime),tspan,y0);
subplot(2,1,1);
plot(t,y(:,1))
title("Membrane Potential")
ylabel("Volts")
xlabel("seconds")
subplot(2,1,2);
Iapp1 = interp1(Itime,Iapp,t,'nearest');
plot(t,Iapp1)
title("Applied Current")
ylabel("Current")
xlabel("seconds")
function [dydt] = HHmodelD(t,y,Iapp,Itime)
GL = 30e-9; GnaK = 12e-6; Gk_Max = 3.6e-6; Ena = 45e-3;
Ek = -82e-3; EL = -60e-3; Cm = 100e-12;
%-------------------------------------------------------
% Where I am modifying the Iapp Value
Iapp = interp1(Itime,Iapp,t,'nearest'); % Interpolate the data set (ITime,Iapp) at time t
%-------------------------------------------------------
dydt1 = (GL/Cm)*(EL - y(1)) +(Gk_Max/Cm)*(y(2)^3) * y(3)*(Ena - y(1)) + (Gk_Max/Cm)*(y(4)^4) *(Ek - y(1)) + Iapp/Cm;
dydt2 = alphaM(y(1))*(1-y(2))-betaM(y(1))*y(2);
dydt3 = alphaH(y(1))*(1-y(3))-betaH(y(1))*y(3);
dydt4 = alphaN(y(1))*(1-y(4))-betaN(y(1))*y(4);
dydt = [dydt1;dydt2;dydt3;dydt4];
function aM = alphaM(V)
aM = ((10^5)*(-V-0.045))/(exp(100*(-V-0.045)) - 1);
end
function aH = alphaH(V)
aH = 70*exp(50*(-V-0.070));
end
function aN = alphaN(V)
aN = ((10^4)*(-V-0.060))/(exp(100*(-V-0.060)) - 1);
end
function bM = betaM(V)
bM = (4*10^3)*exp((-V - 0.070)/0.018);
end
function bH = betaH(V)
bH = (10^3)/(1+exp(100*(-V-0.040)));
end
function bN = betaN(V)
bN = 125*exp((-V-0.070)/0.08);
end
end
Star Strider
Star Strider 2021-2-9
Thank you!
The interpolation is the correct approach to using time-dependent terms in the ODE solution.
If you need to stop the ODE integration and then re-start it to avoid discontinuities, see the documentation section on ODE Event Location for an illustration of how to do this. I’m not certain if the Hodgkin-Huxley models require stopping and then re-starting the integration (there are several versions when I last surveyed them a few years ago) however if that is necessary for your model, using Events is the correct approach.

请先登录，再进行评论。

更多回答（2 个）

Alan Stevens
Alan Stevens 2021-2-9
If you do want to use the discrete values for Iapp, the following works, though you will have to decide if the results are reasonable:
dt = 0.0001;
tspan = 0:dt:0.35;
y0 = [0;0;0;0];
Iapp = zeros(1,numel(tspan));
for i = 2:numel(tspan)
Iapp(i) = Iapp(i-1);
if abs(mod(tspan(i),0.1))<dt/2
Iapp(i) = Iapp(i)+0.22e-9;
end
end
plot(tspan ,Iapp,'.')
xlabel('t'),ylabel('Iapp')
[t,y] = ode45(@(t,y) HHmodelD(t,y,Iapp,dt),tspan,y0);
figure
plot(t,y(:,1))
title("Membrane Potential")
ylabel("Volts")
xlabel("seconds")
function dydt = HHmodelD(t,y,Iapp,dt)
GL = 30e-9; GnaK = 12e-6; Gk_Max = 3.6e-6; Ena = 45e-3;
Ek = -82e-3; EL = -60e-3; Cm = 100e-12;
% Note Iapp in dydt1
ix = round(t/dt+1); % pointer into Iapp for use below
aM = ((10^5)*(-y(1)-0.045))/(exp(100*(-y(1)-0.045)) - 1);
aH = 70*exp(50*(-y(1)-0.070));
aN = ((10^4)*(-y(1)-0.060))/(exp(100*(-y(1)-0.060)) - 1);
bM = (4*10^3)*exp((-y(1) - 0.070)/0.018);
bH = (10^3)/(1+exp(100*(-y(1)-0.040)));
bN = 125*exp((-y(1)-0.070)/0.08);
dydt1 = (GL/Cm)*(EL - y(1)) +(Gk_Max/Cm)*(y(2)^3) * y(3)*(Ena - y(1)) + (Gk_Max/Cm)*(y(4)^4) *(Ek - y(1)) + Iapp(ix)/Cm;
dydt2 = aM*(1-y(2))-bM*y(2);
dydt3 = aH*(1-y(3))-bH*y(3);
dydt4 = aN*(1-y(4))-bN*y(4);
dydt = [dydt1;dydt2;dydt3;dydt4];
end
This results in the following for the membrane potential
Jan
Jan 2021-2-9
编辑：Jan 2021-2-9
Remember that Matlab's ODE integrators are designed to operate on smooth functions. The hard jumps of the parameter will confuse the stepsize estimator. If you are lucky, the integrator stops with an error message. With less luck, the integration provides a final value, which is dominated by rounding errors due to a huge number of integration steps.
If the model changes, stop the integration and restart it:
tspan = 0:0.0001:0.35;
t0 = 0;
y0 = [0;0;0;0];
Iapp = 0;
ready = false;
yAll = [];
tAll = [];
while ~ready
% This time span fro t0 to t0+0.1 including the original
% intermediate steps:
if t0 + 0.1 < tspan(end) % Final interval:
tSpanS = [t0, tspan(tspan > t0 & tspan < t0 + 0.1), t0 + 0.1];
else
tSpanS = [t0, tspan(tspan > t0)];
end
% Perform the integration for this interval:
[t, y] = ode45(@(t,y) HHmodelD(t,y,Iapp), tSpanS, y0);
% Append the solution:
tAll = cat(1, tAll, t);
yAll = cat(1, yAll, y);
% New initial values are former final values:
y0 = y(end, :);
t0 = t(end);
% Advance the parameter
Iapp = Iapp + 0.1;
% Has the final time been reached?
ready = abs(t0 - tspan(end)) < 10 * eps(t0);
end

类别

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

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by