Problem with time span when integrating with ODE45

5 views (last 30 days)
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!

Accepted Answer

William Rose
William Rose on 8 Apr 2021
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 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by