Solve ODE with a time dependent parameter
8 次查看(过去 30 天)
显示 更早的评论
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
0 个评论
采纳的回答
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.
2 个评论
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
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/514302/image.jpeg)
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
另请参阅
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!