function f = stagefour(x,pars)
% return fit of tuned Izi MSN basic and D1/D2 synapse 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 from Stage 1 added to 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 - fixed values
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 - from Stage 3 search
PSPampa = pars(27);
PSPnmda = PSPampa / r_ampa_nmda;
PSPgaba = PSPampa / r_ampa_gaba;
% set intrinsic DA model parameters - from Stage 2 search
DAtype = pars(28);
DA = pars(29);
switch DAtype
case 1
% last parameters are mKIR and mLCA
vrD1 = vr*(1+DA*pars(30));
dD1 = d*(1-DA*pars(31));
cD1 = x(1);
case 2
% last parameter is alpha
kD2 = k * (1-pars(30)*DA);
cD2 = x(1);
end
% 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));
%%% switch D1 and D2 type
switch DAtype
case 1
v(i+1) = v(i) + dt*(k*(v(i)-vrD1)*(v(i)-vt)-u(i) + ...
(Gampa(i+1) .* (Eampa - v(i)))+ (1+cD1*DA)*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)-vrD1)-u(i));
% spikes?
if v(i+1)>=vpeak
v(i)=vpeak; v(i+1)=c; u(i+1)=u(i+1)+dD1;
end
case 2
v(i+1) = v(i) + dt*(kD2*(v(i)-vr)*(v(i)-vt)-u(i) + ...
(1-cD2*DA)*(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
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);