ode23s gives different results when run inside for loop

I am trying to record the spike times of a neuron model that I have written for varying stimulus frequencies. When I run ode23s inside of a for loop, it is giving me different results than when I run it by itself. If you run "loop" and compare the recorded event times to those when you run "single" you will see what I mean. Any suggestions would be appreciated.

Below is the code, each is a separate file

%%%%single %%%%

% set parameters

clear -global all

global  spiketimes a A

tstart = 0; %ODE time vector

tend = 50;

dt = .01;

tspan = tstart:dt:tend;

V0 = -60.865; %initial condition

[meq,neq,Teq] = gates(V0);

y0 = [V0,neq];

%stimulus parameters

lambda = .1; 

mean = 40;

a = 1/4;

A = mean/lambda;

% create vector of synaptic event times

spiketimes = 1/lambda:1/lambda:tend;

% call ode solver

options = odeset('MaxStep', dt,'Events',@events);

[t,y,TE,YE,IE] = ode23s(@ode,tspan,y0,options);

%%%%loop %%%%

% set parameters

clear -global all

global  spiketimes a A

tstart = 0; %ODE time vector

tend = 50;

dt = .01;

tspan = tstart:dt:tend;

V0 = -60.865; %initial condition

[meq,neq,Teq] = gates(V0);

y0 = [V0,neq];

%stimulus parameters

lambda = .05;

mean = 40;

a = 1/4;

A = mean/lambda;

spikes = cell(1,3);

stimulus = cell(1,3);

numspikes = zeros(1,3);

numstim = zeros(1,3);

for i = 1:3

% create vector of synaptic event times

      spiketimes = [1/(lambda + (i-1)*.05):1/(lambda + (i-1)*.05):tend];

      stimulus{i} = spiketimes;

      numstim(i) = length(spiketimes);

% call ode solver

options = odeset('MaxStep', dt,'Events',@events);

[t,y,TE,YE,IE] = ode23s(@ode,tspan,y0,options);

spikes{i} = TE;

numspikes(i) = length(TE);

end

%%%%ode %%%%

% define function, y is the column vector containing V,m,n,h

function [dydt] = ode(t,y)

% set parameters

C = 1; %capacitance

gNa = 20; %conductance

gK= 10;

gL = 8;

VNa = 60; %Nernst potentials

VK = -90;

VL = -78;

%allocate variables

dydt = zeros(2,1);

V = y(1);

n = y(2);

[minf,ninf,Tn] = gates(V);

I = regularcurrent(t);

% write differential equations

dydt(1) = 1/C*(-gNa*minf*(V - VNa) - gK*n*(V - VK) - gL*(V - VL) + I);

dydt(2) = 1/Tn*(ninf - n);

end

%%%%gates %%%%

% define function

function [minf,ninf,Tn] = gates(V)

% allocate variables

Vm = -20; %gating variable parameters

Vn = -45;

km = 15;

kn = 5;

Vmaxn = -30; %time constant parameters

Cbase = 1;

Camp = 30;

sigma = 3;

minf = zeros(1,length(V));

ninf = zeros(1,length(V));

Tn = zeros(1,length(V));

vm = Vm*ones(1,length(V));

vn = Vn*ones(1,length(V));

% do calculations

minf = 1./(1+exp((vm - V)/km));

ninf = 1./(1+exp((vn - V)/kn));

Tn = Cbase + Camp.*exp(-(Vmaxn - V).^2/sigma^2);

%%%%events %%%%

%spike listener

function [value,isterminal,direction] = events(t,y)

value = y(1) + 30;

isterminal = 0;

direction = 1;

end

%%%%regularcurrent %%%%

function [I] = regularcurrent(t)

global spiketimes a A

I = 0;

for i = 1:length(spiketimes)

I = I + A*(1/(a*sqrt(pi))*exp( - (t - spiketimes(i))^2/a^2));

end

Answer 1

[A Jenkins]

I'm not sure I know exactly what your output should look like, but I do see a couple differences between your "single" and "loop" models.

1) In "single" your global is called stimtimes, but in "loop" it is called spiketimes. The regularcurrent() function uses the spiketimes varaible only. [Tip: globals are scary]

2) In "loop" you are trying to create the spiketimes vector as follows:

spiketimes = [1/(lambda + (i-1)*.05):1/(lambda + (i-1)*.05),tend];

but I think you mean to use the colon (:) instead of the comma (,) to create a vector similar to your "single" case

spiketimes = [1/(lambda + (i-1)*.05):1/(lambda + (i-1)*.05):tend];

[Ryan]

Thanks for your response. The differences that you've found are typos that I apparently didn't catch from when I was trying to simplify my code. You are correct that I meant to use the colon, and that the global variable should be the same in both cases. I've left the office for today so I can't try it again, but I'm pretty sure I will see the same problem though.

The output that I'm looking at is the vectors contained in the cell "spikes" from "loop" compared to the vector TE from "single". They contain the times that the neuron voltage crosses -30. For the first iteration of the loop, they are identical, but the next two are not. The differences get worse and worse the more iterations of the loop that I run.

[Ryan]

I edited the above post to reflect the errors that you caught

[A Jenkins]

Well your spiketimes is a function of i. So I would expect your answer to change everytime through the loop since your functions use spiketimes, which uses i, which changes everytime through the loop. Do you want TE to match spiketimes?

It seems to be pretty close (I've rounded so it is easier to read):

% added to the end of the for loop:
fprintf('\nRun %d\n', i)
disp(['spiketimes: ' num2str(round(spiketimes))])
disp(['TE:         ' num2str(round(TE'))])

Results:

Run 1
spiketimes: 20  40
TE:         20  40

Run 2
spiketimes: 10  20  30  40  50
TE:         10  20  30  40  50

Run 3
spiketimes: 7  13  20  27  33  40  47
TE:         6  13  20  26  33  40  46

Run 4
spiketimes: 5  10  15  20  25  30  35  40  45  50
TE:         5  10  15  20  25  30  35  40  45  50

Run 5
spiketimes: 4   8  12  16  20  24  28  32  36  40  44  48
TE:         4   8  12  16  20  24  28  32  36  40  44  48

Perhaps the spiketimes calculation is still wrong?

[Ryan]

A Jenkins,

Thanks again for taking a look. Sorry I wasn't very clear about what I'm looking for. I'll try to explain here.

When I run "loop", I save the spiketimes vectors to the "stimulus" cell, and the TE vectors to the "spikes" cell for each loop iteration. Ultimately I want to compare the vectors in "spikes" to the corresponding TE vectors that I get when I run "single". Each iteration of "loop" should be effectively increasing lambda by .05, so for i going from 1 to 3 I'm expecting that the "spikes" vectors be equivalent to the TE vector from "single" for lambda values .05, .1, and .15 respectively. Likewise the "stimulus" vectors from "loop" should be equivalent to "spiketimes" from "single" for the same respective values. The "spiketimes" vectors are indeed matching up with the "stimulus" vectors. But the "spikes" vectors are not matching up with the respective TE vectors (they are close but they should be exact).

Spiketimes acts as an input function for the ODE. As far as I can tell the exact same calculation should be taking place whether I do it in a loop, or if I manually change the value of lambda in "single" each time.

Let me know if anything remains unclear.

Thanks

[A Jenkins]

Ok, thanks for the detailed explanation. Your difference is because you have different values of A in each case.

• In the single case, you set lambda, then calculate A, so it updates each time.
• In the loop case, you set A and leave it, then change the effective lambda in your loop. So A is only right for the first case, then diverges.

So you should just need to add

A = mean/(lambda + (i-1)*.05);

to your loop.

[Ryan]

Ah thanks I'm glad it was something simple. I appreciate your help.