MATLAB: Ode23s gives different results when run inside for loop

for loopMATLABode23s

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

Best Answer

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];
Related Question