{Simulation of a Fluctuation Test Experiment}
{Semi-stochastic model - mutation is by a Monte Carlo process and population growth is determinisitic}
{Assumes a total volume of 1 ml}
{Run in the Repeat Batch mode for each experiment}
{Variables:     A- ancestral population, B- mutant population - bacteria per ml}
{R - Concentration of the limiting resource- ug/ml.}
{Parameters: va, vb - maximum growth rates }
{ ka, kb, resource concentraton at va/2, vb/2}
{e - conversion efficiency,  ma - maximum mutation rate to B} 
{In this model, the mutation rate declines with the resource concentration}
{The maximum density of the population, A+B, is equal to the quotient of the inital concentration of the resource R(0) and e.}
{For example if init R =500 and e= 5E-7, the maximum density of the population is 1E9}
{To run this program, use the "Batch Run" routine in the Parameter drop down.} 
{Make 20 or so runs in each batch.   You only need print out the 24 hour B population results (assuming stationary phase is reached by then)}



METHOD EULER

STARTTIME = 0
STOPTIME = 24 
DT = 1E-3  {This should be set at 1E-3 or less}   
DTOUT = 24

{Parameters}
k=0.25 {Resource concentration of half maximum growth rate}
e = 5e-7 {Conversion efficiency}
va=1 { Maximum growth rate population A}
vb=1 { Maximum growth rate population B}
ma=1e-9 {Mutation rate A->B}
{A is the dominant population and we neglect mutation from B->A}

{Variables}
init A =1E2 {Density of parental population}
init B =0  {Density of mutants}
init R = 500 {Initial concentration of the limiting resource}

{Differential equations - population growth and resource consumption}
d/dt (A) = (va*R/(k+R))*A - GM/DT
d/dt (B) = (vb*R/(k+R))*B +GM/DT
d/dt (R) = - (vb*R/(k+R))*B*e -  (va*R/(k+R))*A*e 

{Mutation entry mode}
bm =  A*ma*(R/(R+k))*DT 
rm =RANDOM (0, 1)
GM = IF rm < bm THEN 1 ELSE 0