Problem with time span when integrating with ODE45

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!

Answer 1

[William Rose]

@Samuele Bolotta,

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.

[Samuele Bolotta]

Hey William, thanks for your answer! Yes, it's an amazing model indeed.

This is the function:

function dydt = myode(T,Vmnh,TimeInj,IInj,TimeSyn,Gsyn)

%% Constants
ENa=56; % Reversal potential of Sodium
EK=-102; % Reversal potential of Potassium
El=-49; % Leakage reversal potential
EGlu = 0;  % Reversal potential of Glu
EChl = -76;  % Reversal potential of Chl

%% Capacitance
Cm = 0.01;

%% Values of conductances
GL = 0.003; % Leak maximal conductance
GNa = 1.2; % Maximal conductance of Sodium
GK = 0.36; % Maximal conductance of Potassium

%% Conductances times driving force
INa = (GNa * Vmnh(2,1)^3 * Vmnh(4,1)) *(Vmnh(1,:) - ENa);
IK = (GK * Vmnh(3,1)^4) *(Vmnh(1,:) - EK);
Il = GL *(Vmnh(1,:) - El);
IPSC = Gsyn(1,:) * (Vmnh(1) - EChl);
EPSC = Gsyn(2,:) * (Vmnh(1) - EGlu);

%% 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)];

end

And this is how I call it in the program

[T,yF] = ode45(@(T,Vmnh) myode(T,Vmnh,TimeInj,IInj,TimeSyn,Gsyn), tspan, Vmnh);
%% Vmnh are the starting points.

If I only set the starting and the final point of the time span, as you suggested, I get this (image attached).

[William Rose]

I would not use interp1() inside your dydt function. I realize you need something to generate the correct injection current and current from EPSPs an IPSPs. I would do that with a function call, Create functions iInj(t), iE(t),iI(t). Call them inside your dydt() function. I have done the analogous thing in some simulations.

[William Rose]

@Samuele Bolotta,

Is the figure you provided from your program or from somebody else's? I am surpised your program runs at all, given your dydt() structure. I am pretty sure you must combine your integration variables into a single vector. n your dydt(), they are passed as spearate variables. I don;t think that will work.

Separate point: Non-differentiable functions cause problems for ODE solvers, which generally rely on the assumption that the functions involved have continuous first derivatives. Your injection current appears t have sharp corners, i.e. non-differentiable. YOu may need to smooth it out. You can do that by making your functions smoothly ramp up or down over a short interval.

[William Rose]

Maybe I misunderstod the variables in your dydt(). Help me understand the state vector. Is the following true?

state vector = Vmnh(1:4)

Vmnh(1)=V(t)=membrane potential

Vmnh(2)=m(t)

Vmnh(3)=n(t)

Vmnh(4)=h(t)

Why are there not equaitons for , , etc in your dydt() function? I would expect to see

Why is there a minus sign in the total current equation in your code, between the injected current and the other currents? Do they have different sign conventions?

dydt = [((1/Cm)*(IInj-(INa+IK+Il+IPSC+EPSC))); 

If you post the entire program and error messages (if any) when it runs, that will help. Include the functions alpham(V), alphan(V), betam(V), betan(V), etc.

[William Rose]

I had not seen your secnd post whien I made my subsequent posts. I see that you do have equaitons for I_Na, I_K, etc.

It appears from your plots that ode45() gave up after 1.2 msec. The injected cirrent also differs between your two figures. What's up with that?

YOu use

function dydt = myode(T,Vmnh,TimeInj,IInj,TimeSyn,Gsyn)

but I always only pass t and y (where in your case, y is replaced by Vmnh)

function dydt = myode(T,Vmnh)

But this change would require defining TimeInj,IInj,TimeSyn,Gsyn inside function myode(), or declaring them global, in main() and in myode(). Declaring global variables for the ODE function can really slow down the simulation, I have found.

[William Rose]

If you post all the files necessary to run your code, I will try it later.

[William Rose]

@Samuele Bolotta, regarding my earlier quesiton about the sign conventons for currents: I now understand that you are using the "neurophysiological convention" (see here, p.5), according to which positive ion flux out of the cell is considered to be a negative current, if it is through an ion channel or EPSC or IPSC, and it is considered positive if it is injected via the microelectrode.

[Samuele Bolotta]

Hi William, thanks for the feedback. Regarding your last point, yes exactly! And the figure was from my program.

These are the files necessary to run the code. The main file to run is Z1_interpol.m, this should clarify things! Let me know what you think :)

[William Rose]

@Samuele Bolotta, I don't see any files. Can you try again to attach the code files?

[Samuele Bolotta]

Ops sorry! What about now

[William Rose]

@Samuele Bolotta,Thanks! Program seems to run fine. It runs for the full 150 msec with no errors and theoutput yF contains no NaNs. The plot it produces is below. The results are reasonable. For example, the absence of action potentials from 60-100 msec , despite significant current injection, is expected, due to voltage-driven inactivation of Na channels. And the low level of h during this time period shows that this is in fact the reason.

Suggestions:

• Change the plot order to put m and h together in the legend, since they go together physiologically - they both represent Na channel gating (m is activation, h is inactivation by the "ball-and-chain") (see this paper I co-authored, for more on Na channel inactivation.)
• Change the y-axis units on the second plot to nA (at most), to make it seem semi-realistic to physiologists. Amps is too large by 10 orders of magnitude - see previously cited paper.
• Change to a smaller step size, because 0.25 ms does not resolve the Na channel dynamics well. This is evident in the fact that Na channel activation, m, goes from 15% to 85% in a single step, and membrane voltage goes from resting to max in two time steps, at 15 ms.

The sim runs at a constant step size of 0.25 msec. Thank you for teaching me with this example that one may pass all the time points in tspan, instead of just the start and end times. I din't know that was allowed. Then ode45() is no longer an adaptive stepsize routine, but that seems to be OK - just maybe not super-accurate or efficient. I tried replacing the ful array tspan , in the call to ode45(), with just the start and end times, [0 150], as the ode45 documentation calls for. The program actually runs to competion, from t=0 to 150, in just 49 steps, but almost all the values in the output array yF are NaN or Inf. Your method with a fixed stepsize works, so go with it.

Nice job!

[William Rose]

This prgram can be a bit odd. I tried making the synaptic conductances zero by adding a line to Z1_interpol.m:

%% Synaptic conductances
[Gsyn,TimeSyn] = syncon(0.3, 0.15);
Gsyn=zeros(2,length(TimeSyn));      %force Gsyn=0

The program runs to completion without errors, and the outputs look totally normal up to 57.5 msec, but after that time, the output is all Infs and NaNs. There is no physiological reason this should happen. I was able to "fix" this issue by commenting out the interp1 commands for IPSC, EPSC in myode.m and adding adding two lines to set IPSC and EPSC to zero, as shown:

%% Interpolate currents
IInj = interp1(TimeInj,IInj,T,'linear','extrap');
%IPSC = interp1(TimeSyn,IPSC,T,'linear','extrap');
%EPSC = interp1(TimeSyn,EPSC,T,'linear','extrap');
IPSC=0;         %force IPSC=0; no interpolation
EPSC=0;         %force EPSC=0; no interpolation

Then the main program runs fine. See output below.

This indicates that the interp1 commands inside myode.m are causing problems sometimes. I don't know why.

I had a similar result when I tried changing the IInj waveform by modifying injected_current.m. Instead of using random integer currents for specific time windows, I specified the current exactly, and it was still a series of steps. The relevant portion of injected_current.m thus was as follows:

%IInj = [zeros(1,60), ones(1,80) * b, ones(1,90) * c, ones(1,80) * d, ones(1,90) * e, ones(1,80) * f, ones(1,80) * g, ones(1,41) * h] ;

% Use 15 msec of preconditioning current,
% then 1 msec of IInj=4, then 14 msec of zero, then repeat.
% This sequence is similar to one of the original HH experiments. 
IInj = [-0.1*ones(1,60), 4*ones(1,4), zeros(1,56),...
         0.0*ones(1,60), 4*ones(1,4), zeros(1,56),...
         0.2*ones(1,60), 4*ones(1,4), zeros(1,56),...
         0.5*ones(1,60), 4*ones(1,4), zeros(1,56),...
         1.0*ones(1,60), 4*ones(1,4), zeros(1,57)];
TimeInj = 0:0.25:150;

The program runs to completion with no error messages, but it puts out NaNs and Infs after 45 msec. Why does a sequence of random current steps work fine, but not a sequence of determined steps? I fixed this by replacing IInj=nterp1() in myode.m with

%IInj = interp1(TimeInj,IInj,T,'linear','extrap');
IInj=0.02*T;

which causes IInj to ramp up from 0 to 3 over 150 msec. I also added a simimilar line to Z1_interp, so the plot of IInj versus time would be right. Then the program runs fine. See output below.

Again, this indicates that interp1 inside myode() is the problem. I don't know why interp1 inside myode() works sometimes and fails sometimes. But it is definitely the source of trouble.

[William Rose]

@Samuele Bolotta, You might want to try ode23() instead of ode45(). ode45() blew up on me when I tried to use adaptive stepsize, but ode23 worked fine.

Here is a plot of results with a different input pattern, and adaptive sptepsize. All calls to interp1() are gone. The inputs in this case are a 0.2 msec current pulse at t=10, an IPSC at t=25, an EPSC at t=75, and IPSC+EPSC at t=125. Results are physiologically reasonable. Files attached.

Your program is nice! If you want to do the injected current steps to random levels, like you did before, you can do it by adjusting the values in the arrays of step times (ts) and step amplitudes (as) passed to smoothstep().

[Samuele Bolotta]

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!

[William Rose]

@Samuele Bolotta,You're welcome. Good luck.