{Episodic Selection Model}
{Model used in Johnsen, P.J, D. Dubnau, and B. R. Levin (2009) Episodic selection and the maintenance of competence} 
{and natural transformation in Bacillus subtilis  GENETICS  181: 1521-1533}
{Model of antibiotic decay, PD/PK and Transformation}
{Hill function for antibiotic concentration dependent growth}
{Stochastic transformation}
{No resource-concentration dependent competence}
{Random entry of antibiotics}
{Serial Transfer Model -  competition with population not capable of generating competent cells}

METHOD EULER
STARTTIME = 0
STOPTIME=500
DT = 0.001
DTOUT =1
{Parameters}
f =0.01 {Maximum rate at which S->C - competence formation}
fc=0.2 {Maximum rate at which C->S}
k=0.25 {Resource concentration at half maximium growth rate}
e=5e-7 {Conversion efficiency}
mic = 1 {MIC}
vmaxs =1.0 {Maximum growth rate sensitive - competent produce}
vmaxc = 0.001 {Maximum growth rate competent}
vmaxn=1.0 {Maximum growth rate sensitive non-competent}
vmaxct=0.95 {Maximum growth rate resistant transformants}
fmins= -2 {Minimum growth rate sensitive - competent producer}
fminc=-0.01 {maximum kill rate competent}
fminn=-2 {Minimum growth rate non-competent producers}
k=1 {Hill Coefficient}
km = 0.25 {Monod constant}
X=1e-15 {rate constant of transformation}
dd =0.0 {DNA decay rate}
da=0.5 {Antibiotic decay}
d=0.01 {dilution factor}
l=1 {Delay}
e=5E-7
apulse=0.2 {Antibiotic added at pulse}
{Variables}
init S= 1E5{Initial density of sensitive competent producing bacteria 1}
init C =0 {Initial density of competent bacteria}
init N=1e4 {Initial density of non-competent}
init CT=0 {Initial resistant transformed Competent bacteria}
init DNAR=1e4 {Concentration of DNA with transforming resistance genes}
init R = 500 {Initial resource concentration}
init A =0.0 {Initial antbiotic concentration}
init TT=0
init TRANS=0
d/dt(TRANS) = TRX
psi= R/(R+km)
vms= vmaxs*psi
vmc= vmaxc*psi
vmct=vmaxct*psi
fms=fmins*psi
fmc=fminc*psi
fmn=fminn*psi
vmn=vmaxn*psi
d/dt (R) = - e*(S*vmaxs+C*vmaxc+CT*vmaxct+ N*vmaxn) *psi + FR -ZR+AR
psis = ((vms-fms)*(A/mic)^k)/((A/mic)^k - fms/vms)
psic = ((vmc-fmc)*(A/mic)^k)/((A/mic)^k - fmc/vmc)
psin = ((vmn-fmn)*(A/mic)^k)/((A/mic)^k - fmn/vmn)
freal=f*(1-R/(km+R))
fcreal=fc*R/(km+R)
d/dt (S) = vms*S - freal*S +fcreal*C - psis*S +GS-HS {Change in the density of sensitive bacteria}
d/dt (C) =vmc*C - fcreal*C+freal*S  - psic*C +GC-HC {Change in density of competent bacteria}
d/dt (N)= vmn*N - psin*N +GN-HN
d/dt (CT)= vmct*CT + Tr  +GCT - HCT {Change in density of transformed cells }
d/dt (DNAR)= - DNAR*dd
d/dt (A) = -da*A + ADD - REM
NT=S+C+CT+N
bm = X*DNAR*C*DT 
rm =RANDOM (0, 1)
Tr = IF rm < bm THEN 1/DT ELSE 0
CS= vms-psis
CC=vmc-psic
CN=vmn-psin

{Serial Transfer stuff}
d/dt (TT)=1-GT
d/dt (TRANS) = TRX
init R = 500
init TT =0
AADD=10
FR = IF TT > 24 THEN 500/DT  ELSE 0
ZR= IF TT > 24 THEN R/DT ELSE 0
GS = IF TT > 24 THEN S*d/DT ELSE 0
HS = IF TT > 24 THEN S/DT ELSE 0
GC = IF TT > 24 THEN C*d/DT ELSE 0
HC = IF TT > 24 THEN C/DT ELSE 0
GN = IF TT > 24 THEN N*d/DT ELSE 0
HN = IF TT > 24 THEN N/DT ELSE 0
GCT = IF TT > 24 THEN CT*d/DT ELSE 0
HCT = IF TT > 24 THEN CT/DT ELSE 0
GT= IF TT >24 THEN 24/DT ELSE 0
TRX= IF TT >24 THEN (TRANS+1)/DT ELSE 0
AR = IF R<0.0001 THEN 0.0001/DT ELSE 0
ba=apulse
ra=RANDOM(0,1)
ADD= IF TT>24 AND ra<ba THEN AADD/DT ELSE 0
REM= IF TT>24 AND ra<ba THEN 2*A/DT ELSE 0