function f = stagetwo(x,pars)

% return fit of Izi MSN D1 or D2 type neuron (intrinsic) to corresponding Moyer et al 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);
d1int = [0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	2.0000	6.0000	8.0000	10.0000	12.0000	14.0000	16.0000	16.0000	18.0000];
d2int = [2.0000	2.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	16.0000	18.0000	18.0000];


% pars = [k,a,b, c, vr, vt, vpeak, T, dt, f_start, f_end, f_time, Istore, C, d, k];
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, vt, d]
C = pars(13); vt = pars(14); d = pars(15);

% dopamine stuff
DAtype = pars(16);
DA = pars(17); 

% multipliers
if DAtype == 1
    dD1 = d*(1-DA*x(2));
    vrD1 = vr*(1+DA*x(1));
else
    kD2 = k * (1-x(1)*DA);
end

% 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
        switch DAtype
            case 1
                %--- D1
                v(i+1) = v(i) + dt*(k*(v(i)-vrD1)*(v(i)-vt)-u(i) + injScaled(loop))/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
                %--- D2
                v(i+1) = v(i) + dt*(kD2*(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
            otherwise
                error('Unknown receptor type')
        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
% %% SSE of f-I
% switch DAtype
%     case 1
%         fIerror = sum((fI-d1int).^2); % SSE
%     case 2
%         fIerror = sum((fI-d2int).^2); % SSE
% end

%% summed relative error f-I
switch DAtype
    case 1
        norm = d1int; norm(norm == 0) = 1;
        fIerror =  sum(abs(fI-d1int)./norm); % summed relative error
    case 2
        norm = d2int; norm(norm == 0) = 1;
        fIerror = sum(abs(fI-d2int)./norm); % summed relative error
end


f = fIerror;