{Model of antibiotic decay, PD/PK and Physiological Tolerance}
{Hill function for antibiotic concentration dependent growth}
{Balaban et al phenotypic switch model}
{Limit to population growth}
METHOD RK4
STARTTIME = 0
STOPTIME=240
DT = 0.02
DTOUT =.1
fmaxs=1  {Maximum growth rate sensitive}
fmins= -10 {Minimum growth rate sensitive}
fmaxt=0.1  {Maximum growth rate tolerant}
fmint=-0.2
kmax = 1e20 {Logistic saturation density}
mic = 1 {MIC same for S and T}
k=1 {Hill Coefficient - same for S and T}
x =0.0 {Rate of conversion of sensitive to tolerant}
y =0.0 {Rate of conversion of tolerant to sensitive}
init A =10 {Initial Antibiotic concentration}
amax =10 {Antibiotic added at dose interval}
init S=10{Initial density of sensitive bacteria}
init T=0 {Initial density of tolerant bacteria}
d =.5 {Antibiotic decay rate}
d/dt (A) = -d*A + ADD  {change in the concentration of the antibiotic}
psis = ((fmaxs-fmins)*(A/mic)^k)/((A/mic)^k - fmins/fmaxs)  {Hill for S }
psit= ((fmaxt-fmint)*(A/mic)^k)/((A/mic)^k - fmint/fmaxt)  {Hill for T}
d/dt (S) =(S*fmaxs - x*S + y*T - psis*S)*(1-N/kmax)  {Change in the density of sensitive bacteria}
d/dt (T) =(T*fmaxt  -  y*T + x*S - psit*T)*(1-N/kmax)  {Change in the density of the tolerant bacteria}
N=S+T
RAT=T/S
dose =8{Dosing interval Lambda}
init TT=0 {Time indicator}
d/dt (TT) = 1-GT*2
ADD = IF TT > dose THEN PULSE(amax*2, TIME, 21) ELSE 0  {Dose added}
GT= IF TT >dose THEN PULSE(dose, TIME, 21) ELSE 0