function f = stagethree(x,pars)
% return fit of tuned Izi MSN basic type neuron to Moyer et al f-f plot data
% approximate r values to get same overall synaptic Hz...
N_ctx = 84; alpha_ctx = 0; r_ctx = [4:0.5:8]; % r_ctx = 8;
N_gaba = 84; alpha_gaba = 0; r_gaba = r_ctx; % r_gaba = 8;
% pars = [k,a,b, c, vr, vpeak, T, dt, f_start, f_end, f_time, Istore, C, vt, d];
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);
% newly tuned parameters added to end of pars list
C = pars(13); vt = pars(14); d = pars(15);
% fitted linear function to Moyer data
B1 = pars(16); B2 = pars(17);
% synaptic parameters
r_ampa_nmda = pars(18); r_ampa_gaba = pars(19);
Egaba = pars(20); Enmda = pars(21); Eampa = pars(22); Mg = pars(23);
ts_ampa = pars(24); ts_nmda =pars(25); ts_gaba = pars(26);
% set conductances
PSPampa = x(1);
PSPnmda = PSPampa / r_ampa_nmda;
PSPgaba = PSPampa / r_ampa_gaba;
% run simulation
t = 0:dt:T;
n = length(t); % number of time points
nHz = numel(r_ctx);
SynExp_ampa = exp(-dt / ts_ampa);
SynExp_nmda = exp(-dt / ts_nmda);
SynExp_gaba = exp(-dt / ts_gaba);
for loop = 1:nHz
v = vr*ones(1,n); u=0*v;
Ggaba = zeros(1,n);
Gampa = zeros(1,n);
Gnmda = zeros(1,n);
% generate the spike trains
Sctx = spkgen([0:dt:T], N_ctx, r_ctx(loop), alpha_ctx);
Sgaba = spkgen([0:dt:T], N_gaba, r_gaba(loop), alpha_gaba);
S = Sctx + Sgaba;
for i = 1:n-1
%%% update synaptic input
Gampa(i+1) = Gampa(i) + (PSPampa .* Sctx(i)./ts_ampa);
Gampa(i+1) = Gampa(i+1) * SynExp_ampa;
Gnmda(i+1) = Gnmda(i) + (PSPnmda .* Sctx(i)./ts_nmda);
Gnmda(i+1) = Gnmda(i+1) * SynExp_nmda;
Ggaba(i+1) = Ggaba(i) + (PSPgaba .* Sgaba(i) ./ ts_gaba);
Ggaba(i+1) = Ggaba(i+1) * SynExp_gaba;
% update NMDA plug
B_nmda = 1 ./ (1 + (Mg/3.57) * exp(-v(i)*0.062));
%%% unmodified
v(i+1) = v(i) + dt*(k*(v(i)-vr)*(v(i)-vt)-u(i) + (Gampa(i+1) .* (Eampa - v(i))) ...
+ B_nmda*(Gnmda(i+1) .* (Enmda - v(i))) + (Ggaba(i+1) .* (Egaba - v(i))) )/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
brate(loop) = sum(S) / (T*1e-3);
ff(loop) = sum(v(f_start:f_end) == vpeak) ./ f_time;
temp = find(v == vpeak); isis = diff(temp)*dt;
if isis ffISI(loop) = 1000./mean(isis); else ffISI(loop) = 0; end
end
%-------------- compute fit
% calculate expected value from fitted function
expff = B1 + B2 .* brate;
expff(expff < 0) = 0; % rectify fit!!
% f = sum((ff-expff).^2); % SSE
% SRE
norm = expff; norm(norm == 0) = 1;
f = sum(abs(ff-expff)./norm);