{1- Lytic Phage Model - Batch Culture with Time Delays}
{Three populations of bacteria, S, R1, and R12, and two populations of phage, P1 and P2}
{S is sensitive to all phage, R1 sensitive to P2, and R12 resistant to all phages}
{Model of host range mutations in a single phage}
METHOD EULER
STARTTIME = 0
STOPTIME=80
DT =0.0005 {step size}
DTOUT =0.1 {Output}

Sx=DELAY(S,x) {Sensitive at t-x}
R1x=DELAY(R1,x) {R1 at t-x}
MS1x =DELAY(MS1,x) {S infected with 1 at t-x} 
MS2x=DELAY(MS2,x) {S infected with 2 at t-x}
MR1x= DELAY(MR1,x) {R1 infected with 2 at t-x}
P1x= DELAY(P1,x) {P1 phage at t-x}
P2x= DELAY(P2, x) {P2 phage at t-x}

d/dt (G) =  - (v*S+ v1*R1 + v12*R12)*(G/(G+k))*e
d/dt (S) = v*(G/(G+k))*S - d1*S*P1 - d2*S*P2- NR1
d/dt (R1) = v1*(G/(G+k))*R1 - d2*R1*P2 - NR12 + NR1
d/dt (R12) = v12*(G/(G+k))*R12  + NR12

d/dt (MS1) = d1*S*P1 - d1*Sx*P1x
d/dt (MS2) = d2*S*P2 - d2*Sx*P2x
d/dt (MR1)=  d2*R1*P2 -d2*R1x*P2x

d/dt (P1) = MS1x*b1  -NP12 - dm1*P1
d/dt (P2) = MS2x*b2 +MR1x*b2  -dm2*P2 +NP12

init S =1e6 {Sensitive bacteria}
init P1 =1e5 {Phage1}
init P2 = 0 {phage 2}
init R1=0 {Resistant to P1}
init R12=0 {Resistant to P1 and P2}
init G=1000 {resource}
init MS1=0 {S1 infected with P1}
init MS2=0 {S2 Infected with P2}
init MR1=0 {R1 infected with P2 }

v = 1 {Max growth S }
v1 =.9 {Max growth R1}
v12=0.7 {Max growth R2 }
d1=1e-8 {Adsorption of P1}
d2=1e-8 {Adsorption of P2}
b1=100 {Burst size P1}
b2=100 {burst size P2}


e=5E-7 {Conversion Efficiency}
x=0.1{Latent period}
k=0.25 {Concentration for 1/2 max growth}
dm1=0.1 {rate of mortality of the phage P1}
dm2=0.1 {rate of mortality of phage P2}
vol=10 {Vol of the culture}

mutSR1 = 1e-8 {mutation rate from S to R1}
mutR1R12 = 1e-8  {mutation rate from R1 to R12}
mutP1P2 = 1E-8 {mutation rate from P1 to P2}

ISR1= mutSR1*S*vol*DT  
IR1R12=mutR1R12*R1*vol*DT
IP1P2=mutP1P2*P1*vol*DT

xr1 = RANDOM (0,1) 
xr2 = RANDOM (0,1)
xr3 = RANDOM (0,1)
NR1= IF xr1 < ISR1 THEN PULSE(1/vol,TIME,21) ELSE 0
NR12= IF xr2 < IR1R12 THEN PULSE(1/vol,TIME,21) ELSE 0
NP12= IF xr3 < IP1P2 THEN PULSE(1/vol,TIME,21) ELSE 0