{Hospital Model Program - Extension of that in Bergstrom, Li and Lipsitch PNAS 2004}
{Three antibiotics no cross or double resistance}

METHOD EULER
STARTTIME = 0
STOPTIME=400
DT = 0.01
DTOUT=1 {Print Intervals}
init S =.70 {Initial density of Susceptible}
init X=.20 {Initial density uninfected}
init R1=.03333 {Initial density hosts with bacteria resistant to 1and susceptible to 2 and 3}
init R2=.03333 {Initial density of hosts with bacteria resistant to 2 and susceptible to 1 and 3}
init R3=.03333{Initial density of hosts with bacteria resistant to 3 and susceptible to 1 and 2}
init TT=0
u = 0.1  {rate of patient turnover}
m = 0.7 {fraction of entry of the S type}
m1= 0.03333 {fraction of entry of the R1 type}
m2 = 0.03333 {fraction of entry of the R2 type}
m3 =0.03333 {fraction of entry of the R3 type}
tsum =0.5 {Sum of use of drugs 1, 2 and 3} 	
b =1 {Colonization rate constant}
cycle= 80 {Cycle time}
c1 = 0 {Fitness cost for 1 in the absence of treatment with 1}
c2=0 {Fitness cost for 2 in the absence of treatment with 2}
c3=0 {Fitness cost for 3 in the absence of treatment with 3}
g = 0.03 {rate of loss of colonization of susceptible in the absence of treatment}
si=0.25 {Rate of secondary colonization relative to primary}
a =0 {Compliancy for cycling, a=1 all receive current drug}
	 {a=0, equal probability for each drug, (0 < a < 1) fraction}
             {current drug and (1-a) random with equal for both}
t1=tsum*(1-a)*0.3333+(XT1*a*tsum)*DT
t2=tsum*(1-a)*0.3333+(XT2*a*tsum)*DT
t3=tsum*(1-a)*0.3333+(XT3*a*tsum)*DT
d/dt (S) = u*(m-S)-(t1+t2+t3+g)*S +b*S*X+si*b*(c1*R1+c2*R2+c3*R3 )*S
d/dt(R1) =  u*(m1-R1)-(t2+t3+g)*R1+b*(1-c1)*R1*X-si*b*(c1*S+(c1-c2)*R2)*R1-si*b*(c1*S+(c1-c3)*R3)*R1
d/dt(R2)  = u*(m2-R2)-(t1+t3+g)*R2+b*(1-c2)*R2*X-si*b*(c2*S+(c2-c1)*R1)*R2-si*b*(c2*S+(c2-c3)*R3)*R2
d/dt (R3) = u*(m3-R3)-(t1+t2+g)*R3+b*(1-c3)*R3*X-si*b*(c2*S+(c3-c1)*R1)*R3-si*b*(c3*S+(c3-c2)*R2)*R3
d/dt (X) =u*(1-m-m1-m2-m3-X)+(t1+t2+t3+g)*S + (t2+t3+g)*R1+(t1+t3+g)*R2 +(t1+t2+g)*R3  -b*X*(S+(1-c1)*R1+(1-c2)*R2+(1-c3)*R3)
d/dt (TT) = 1-GT*3
NT=S+R1+R2+R3+X
PR=(R1+R2+R3)/(S+R1+R2+R3+X)
GT= IF TT >3*cycle THEN PULSE(cycle, TIME, 21) ELSE 0
XT1= IF TT < cycle THEN PULSE (1, TIME, 21) ELSE 0
XT2= IF TT > cycle and TT< 2*cycle THEN PULSE (1, TIME, 21) ELSE 0
XT3= IF TT > 2*cycle THEN PULSE (1, TIME, 21) ELSE 0