Problem with time span when integrating with ODE45

25 次查看(过去 30 天)
I am using ODE45 to integrate a system of differential equations and show their dynamics across 150ms. I am struggling with being precise about the time span of integration. More precisely, I have four vectors related to time:
1)
% This is the vector I pass to my ode solver as the time span of
% integration
dt = 0.25;
time = 0:dt:150;
2)
T is the output of the ODE solver and it's correctly a 1x601 vector.
3)
TimeInj, which is the time span of IInj, one of the time-dependent terms that I am interpolating.
TimeInj = 0:0.5:150; % Resulting in a 1x301 vector
4)
TimeSyn, the time span of IPSC and EPSC, the other two time-dependent terms that I am interpolating
TimeSyn = 0:0.25:150; % Resulting in a 1x601 vector
So, each vector has the same time span of 150ms, with one difference: the timestep is 0.5 for TimeInj, resulting in a 1x301 vector and 0.25 for TimeSyn, resulting in a 1x601 vector.
This is the part of my solver where I'm interpolating and integrating:
function dydt = myode(T,Vmnh,TimeInj,IInj,TimeSyn,Gsyn)
..........................
%% Interpolate currents
IInj = interp1(TimeInj,IInj,T);
IPSC = interp1(TimeSyn,IPSC,T);
EPSC = interp1(TimeSyn,EPSC,T);
%% Integration
dydt = [((1/Cm)*(IInj-(INa+IK+Il+IPSC+EPSC)));
alpham(Vmnh(1))*(1-Vmnh(2,1))-betam(Vmnh(1))*Vmnh(2,1);
alphan(Vmnh(1))*(1-Vmnh(3,1))-betan(Vmnh(1))*Vmnh(3,1);
alphah(Vmnh(1))*(1-Vmnh(4,1))-betah(Vmnh(1))*Vmnh(4,1)];
However, as you can see from the image I attached, the timespan for Current IInj and Synaptic Conductances EPSC and IPSC is correctly 150, while the one for the first and four subplot is not correct and has lots of NaN values starting from a certain point. I think I'm doing something wrong with either the time steps or I'm not respecting some rules with the time spans. Thanks!

采纳的回答

William Rose
William Rose 2021-4-8
I think it would be easier to understand your issue if you show your main program that calls ode45(). YOu want to pass to ode45() a 2-elemnt vector containing the start and end time - that's all. You do not pass the step size or a long vector f every time point. The routine will figure out the best stepsize to use in order to achieve "good" accuracy. (You can override its default definition of "good" accuracy but this is rarely necessary in my experience.) I like your model. The HH equaitons are a beautiful accomplishment. It is amazing that H&H figured out what they did, before ion channels were even discovered.
  15 个评论
Samuele Bolotta
Samuele Bolotta 2021-4-11
Hello William, sorry for the late reply but I took some time to process all the feedbacks you gave. Also, I want to thank you for your time.
Using ODE23 instead of ODE45 is indeed a winning choice - and it's interesting that by simply using ODE23 in the code I posted the other day (without changing anything else) instead of ODE45, all the issues you were mentioning disappear. You can set the synaptic conductances to 0 and it also looks like you can keep the interpolation in myode, without the program breaking down! As for your way of injecting current, it is definitely more elegant. Thanks!

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by