function result = issteadySB(model)
% issteadySB: determines whether SBmodel follows the 
% quasi-steady state assumption, i.e., all loops closed

[ic] = SBinitialconditions(model);
% [pn,pv] = SBparameters(model);
% 
% ic = ones(size(ic)).*1e-9;
% pv = ones(size(pv)).*1e-3;
% 
% model = SBinitialconditions(model,ic);
% model = SBparameters(model,pn,pv);

time = 1e+3;
options.AbsTol = 1e-45;
options.RelTol = 1e-6;

ic(ic==0) = options.AbsTol;

result = false;
epsilon = 0.001;

while true
   warning('off','all');
   disp(['Simulating for ',num2str(time),' s']);
   o = SBsimulate(model,'ode15s',time,ic,options);
   ic2 = o.statevalues(end,:);
   ic2 = ic2(:);
   delta = abs((ic2-ic)./ic)./time.*3600;
   if all(le(delta,ones(size(ic)).*epsilon))
      result = true;
      break;
   else
      disp('Steady state not reached.');
      disp(['Max delta: ',num2str(max(delta))]);
      ic = ic2;
      time = time * 2;
   end
end