function dHRF = nda_hrf_f(t, dt, HRF, n, e_, Ts, Tf, Eo, Go, b, Tw, Tmtt, Alpha, Beta)

tindex=round(t/dt+1,0);

fa = HRF(1:2,1);   % Relative cerebral blood flow in the arteriole
v = HRF(3:4,1);     % Relative cerebral blood volume
s = HRF(5:6,1);     % Vasoactive signal
a = HRF(7:8,1);   % Relative pressure to volume relation constant
q = HRF(9:10,1);  % Relative amount deoxy-hemoglobin in balloon/windkessel compart


fw = v.^(Alpha + Beta)./a;
E = (1-Go)*(1-(1-Eo/(1-Go)).^(1./fa));     % Oxygen extraction fraction Zheng 2002 (OEF)
% E = 1-(1-Eo).^(1./fa);     % Oxygen extraction fraction (OEF)
dv = (fa - fw)/Tmtt;

dHRF = zeros(10,1);
dHRF(1:2,1) = s;
dHRF(3:4,1) = (fa - fw)./Tmtt;
dHRF(5:6,1) = (e_*(n(:,tindex)-1) - s/Ts - (fa - 1)./Tf);
dHRF(7:8,1) = (- a + exp(-b*dv))./Tw;
dHRF(9:10,1) = (E/Eo .* fa - fw.*q./v)/Tmtt; 


end