function dx = SEImodel4(t,x,par)
% 
% SEI Model 3:
% S' = -beta*c/(K+N)*SI*chi_[0,10]
% E' = beta*c/(K+N)*SI*chi_[0,10] - gamma*E
% I' = gamma*E-d*I
% where
% beta = probability of transmission
% c = limiting contact rate (max)
% K = density at which half maximum contact rate is achieved
% 1/gamma = incubation period
% d = disease-induced death rate
% and
% chi_[0,10] is 1 for t in [0,10] and 0 otherwise.


% set parameter values:
beta = par.beta;
gamma = par.gamma;
d = par.d;
c = par.c;
K = par.K;

% initial pre-infected newts removed after some cohabitation time:
cohab = par.cohab;
Iinit = par.Iinit*heaviside(cohab-t);

% set state variables;
S = x(1);
E = x(2);
I = x(3);
N = S+E+I+Iinit;

% % At the end of the cohabitation period,
% % remove original infected newts
% if t == cohab
%     I = x(3) - I0; 
% end

% initialize dx
dx=zeros(5,1);

% state equations
dx(1) = -beta*c/(K+N)*S*Iinit;
dx(2) = beta*c/(K+N)*S*Iinit - gamma*E;
dx(3) = gamma*E - d*I;
dx(4) = gamma*E; % cumulative number positive
dx(5) = d*I;
end

