% Function implementing the KNDY network model

function y = modelPopulation_v103(t, x, params)

%set fixed model parameterss
d1 = params.d1;
d2 = params.d2;
d3 = params.d3;
p1_basal = params.p1_basal;
p2_basal = params.p2_basal;
maxRate = params.maxRate;
KD = params.KD;
KN = params.KN;
Kr1 =params.Kr1;
Kr2 =params.Kr2;
KI = params.KI;
p3_external =  -log(2/(params.p3_external*60/maxRate*d3 + 1)- 1) ;
y = zeros(size(x));

%% set the parameters being inferred
p1 = 10.^params.fitparams(1);
p2 = 10.^params.fitparams(2);
p3 = 10.^params.fitparams(3);
ee = 10.^(params.fitparams(9)-params.fitparams(3))/maxRate;%
p3Basal = -log((1-10.^(params.fitparams(4)))./(1+10.^(params.fitparams(4))));%;params.p3Basal;%params.fitparams(5);%params.p3Basal;% params.fitparams(5);%params.p3Basal;

%% state variables
Dyn = x(1);
NKB = x(2);
r   = x(3);

%equation for Dyn
y(1) =  p1_basal + p1.*r.^2./(r.^2 + (Kr1).^2) - Dyn*d1;
%equation for NKB
y(2) =  p2.*r.^2./(r.^2 + Kr2.^2)*KD.^2./(Dyn.^2 + KD.^2) + p2_basal - d2*NKB;
%equation for rate of firing

% I = basal + w * r 
% synaptic weight w_i ranges in [0 1] and is a Function of NKB
I = p3_external + p3Basal + p3*( ee + (NKB.^2)./(NKB.^2 + KN.^2)  ).*r;
y(3) =  maxRate * (KI./(exp(-I)+KI-1)-1) - d3*r;

end