function yp = omicronodefileBPEAbaseline(t,y,params);

S = y(1);
E = y(2);
Ev = y(3);
Ev2 = y(4);
Ev3 = y(5);
Ev3O = y(6);
I = y(7);
Iv = y(8);
Iv2 = y(9);
Iv3 = y(10);
Iv3O = y(11);
R = y(12);
RO = y(13);
H = y(14);
D = y(15);
CI = y(16);
V = y(17);
CItotal = y(18);

betabar = params(1);
kappa = params(2);
zeta = params(3);
gamma = params(4);
sigma = params(5);
eta = params(6);
etaO = params(7);
seasonalsize = params(8);
seasonalposition = params(9);
fatiguesize = params(10);
fatiguemean = params(11);
fatiguesig = params(12);
betabarv = params(13);
tv = params(14);
Evbar = params(15);
lambda = params(16);
xi = params(17);
betabarv2 = params(18);
tv2 = params(19);
Evbar2 = params(20);
betabarv3 = params(21);
tv3 = params(22);
Evbar3 = params(23);
nuO = params(24);
etaO2 = params(25);
indexparams = params(26);
xiO = params(27);



    psi = seasonalsize*(cos((t+seasonalposition)*2*pi/365)-1)/2;
    kappa = kappa*(1 -normcdf(t,fatiguemean,fatiguesig))+fatiguesize*kappa*normcdf(t,fatiguemean,fatiguesig);
    
    %turn on this line for faster vaccines
    lambda = (t>321)*lambda;
   
    eta = (t<301)*eta + (t>301)*0.5*eta;
  
    etaO = (t<750)*etaO+(t>750)*(t<1095)*etaO/1.2 + (t>1095)*etaO/2;

%new deaths are given by the flow out of the H compartment and that
%determines behavior

beta = betabar*exp(-kappa*zeta*H+psi);
betav = betabarv*exp(-kappa*zeta*H+psi);
betav2 =  betabarv2*exp(-kappa*zeta*H+psi);
betav3 =  betabarv3*exp(-kappa*zeta*H+psi);

%now compute the flows
x = zeros(18,1);

% for S
x(1) = -beta*S*I - betav*S*Iv - betav2*S*Iv2 - betav3*S*(Iv3 + Iv3O) - lambda*S + xi*(R+V) + xiO*RO;
%for E
x(2) = beta*S*I - sigma*E;
%for EV
x(3) = betav*S*Iv - sigma*Ev + (t<(tv+5))*(t>tv-5)*Evbar;
%for EV2
x(4) = betav2*S*Iv2 - sigma*Ev2 + (t<(tv2+2))*(t>tv2)*Evbar2;
%for EV3
x(5) = betav3*S*(Iv3 + Iv3O) - sigma*Ev3 + (t<(tv3+2))*(t>tv3)*Evbar3;
%for EV3O
x(6) = betav3*nuO*(R+V)*(Iv3 + Iv3O) - sigma*Ev3O;
%for I
x(7) = sigma*E - gamma*I;
%for Iv
x(8) = sigma*Ev - gamma*Iv;
%for Iv2
x(9) = sigma*Ev2 - gamma*Iv2;
%for Iv3
x(10) = sigma*Ev3 - gamma*Iv3;
%for Iv3O
x(11) = sigma*Ev3O - gamma*Iv3O;
%for R
x(12) = (1-eta)*gamma*(I+Iv+Iv2) - xi*R - betav3*nuO*R*(Iv3 + Iv3O);
%for RO
x(13) = (1-etaO)*gamma*Iv3 + (1-etaO2)*gamma*Iv3O - xiO*RO; 
%for H
x(14) = eta*gamma*(I+Iv+Iv2) + etaO*gamma*Iv3 + etaO2*gamma*Iv3O - zeta*H;
%for D
x(15) = zeta*H;
%for CI cumulative Omicron new infections
x(16) = betav3*S*(Iv3 + Iv3O);
%for V cumulative vaccinated prior to infection
x(17) = lambda*S - xi*V - betav3*nuO*V*(Iv3 + Iv3O);
%for cumulative first infections
x(18) = beta*S*I+betav*S*Iv+betav2*S*Iv2+betav3*S*(Iv3 + Iv3O)+betav3*nuO*V*(Iv3 + Iv3O) - xi*CItotal;

yp = x;