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

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

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


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

dHRF = zeros(6,1);

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

end