{2- Simulation - Model of host defense and phage and bacterial control of an acute infection}
{Partial to full resistance}
{Mutation to resistance}
{Time delays for the latent period}
METHOD EULER
STARTTIME = 0
STOPTIME=50
DT = 0.01 {Step size}
DTOUT = 0.1 {Output interval}
MP=DELAY(M,x)  {Density of infected senstive bacteria at time t-x}
PP=DELAY(P,x) {Density of phage present at time t-x}
SP=DELAY(S,x) {Density of uninfected sensitive bacteria present at time t-x}
MRP=DELAY(MR,x) {Density of resistant infected bacteria at time t-x}
NRP=DELAY(NR,x) {Density of unifected resistant bacteria at time t-x}
d/dt (I) = z*I*(S+NR)/(S+NR+t) {change in immunity}
d/dt (R) =  - (vs*S + vr*NR)*(R/(R+k))*e {Change in resource concentration}
d/dt (S) = vs*(R/(R+k))*S - d*S*P -GM - y*I*S {Change in the density of sensitive bacteria}
d/dt (M) = d*S*P - d*SP*PP {Change in the density of phage infected infected bacteria}
d/dt (MR) = d2*NR*P -d2*NRP*PP {Change in the density of infected resistant bacteria}
d/dt (NR) = (vr*R/(R+k))*NR -d2*NR*P +GM -y*I*NR {Change in the density of reistant bacteria}
d/dt (P) =MP*bv+MRP*bv - dp*P {Change in the density of free phage}
init S =1E6 {Initial density of susceptible bacteria}
init P =0 {Initial density of phage}
init NR=1 {Density of resistant bacteria}
init R=1000 {Initial resource concentration}
init M=0  {initial density of phage infected sensitive bacteria}
init MR=0 {Initial density of phage infected resistant bacteria} 
init I =1 {Initial concentration of immune controller}
vs = 1 {Maximum growth rate sensitive bacteria}
vr =.9 {Maximum growth rate resistant bacteria}
e=5E-7 {Conversion efficiency}
x=0.5 {Latent period}
d=1E-8 {Adsorption rate - sensitive}
d2=1e-12 {Adsorption rate - resistant}
k=0.25 {resource at half maximum rate}
bv=100 {Burst size}
mut = 0 {Mutation rate to resistance}
vol=10 {Volume of the habitat}
dp=0.2 {Death rate of the phage}
z=1 {Maximum rate of growth of the immune response}
t=1e4 {cell density were growth rate of immune response is half maxium}
y= 1E-6 {Immune kill parameter}
bmut=mut*S*vol*dt {probability of a mutant during time interval}
rm = RANDOM (0,1)  {Random number generated}
GM= IF rm < bmut THEN PULSE(1/vol,TIME,21) ELSE 0  {Mutants added}