function f = stageone(x,pars)
% return fit of Izi MSN neuron to Moyer et al basic data
% Moyer data - injection current and corresponding output
inj = [0.22 0.225 0.23 0.2350 0.2400 0.2450 0.2500 0.2550 0.2600 0.2650 0.2700 0.2750 0.2800 0.2850 0.2900 0.2950 0.3000];
injScaled = round(inj .* 1e3);
unmod = [0.0000 0.0000 0.0000 2.0000 4.0000 4.0000 6.0000 6.0000 8.0000 8.0000 10.0000 10.0000 12.0000 12.0000 14.0000 14.0000 16.0000];
% pars = [k,a, b, c, vr, vt, vpeak, T, dt, f_start, f_end, f_time,Istore];
k = pars(1); a = pars(2); b = pars(3); c = pars(4); vr = pars(5); vpeak = pars(6);
T = pars(7); dt = pars(8); f_start=pars(9); f_end=pars(10); f_time=pars(11);Istore = pars(12);
% x = [C, vt, d]
C = x(1); vt =x(2); d = x(3);
% run simulation
t = 0:dt:T;
n = length(t); % number of time points
nInj = numel(inj);
for loop = 1:nInj
v = vr*ones(1,n); u=0*v;
for i = 1:n-1
%--- unmodulated
v(i+1) = v(i) + dt*(k*(v(i)-vr)*(v(i)-vt)-u(i) + injScaled(loop))/C;
u(i+1) = u(i) + dt*a*(b*(v(i)-vr)-u(i));
% spikes?
if v(i+1)>=vpeak
v(i)=vpeak; v(i+1)=c; u(i+1)=u(i+1)+d;
end
end
% time to first spike
temp = find(v == vpeak); isis = diff(temp)*dt;
if temp tfs(loop) = temp(1) * dt; else tfs(loop) = nan; end % time in ms
fI(loop) = sum(v(f_start:f_end) == vpeak) ./ f_time;
if isis fI1st(loop) = 1000./isis(1); else fI1st(loop) = 0; end
end
%-------------- compute fit
% %% RMSE of fI
% fIerror = sqrt(sum((fI-unmod).^2)/nInj); % RMSE
norm = unmod; norm(norm == 0) = 1;
fIerror = sum(abs(fI-unmod)./norm)/nInj; % mean summed relative error
f = fIerror;