{Simulation of a Fluctuation Test Experiment}
{Run in the Repeat Batch mode for each experiment 20 or so replicas}
{In this model, the mutation rate declines with the resource}
{The final cell density in this simulation is the} {quotient of the initial resoruce concentration}
{init R, and the convestion efficiency e.}
{For example, if init R=500 and e=5E-7, the final}
{density of each culture would be 500/5E-7= 1E9.}
{To calculate the mutation rate, you need the final}  {density of these cultures and the}
{number of mutants in each culture.} 
METHOD EULER
STARTTIME = 0
STOPTIME = 24
DT = 1E-4   {Step size}
DTOUT =24  {Output time}
ka =0.25 {resource concentration at 1/2 the maximum} {rate of the ancestral population}
kb =0.25 {resource concentration at 1/2 the maximum}
{rate of growth of the mutant population}
ea = 5e-7 {Conversion efficiency ancestral population}
eb=5e-7 {Conversion efficiency mutant population}
va=1 {Maximum growth rate of ancestral cell population}
vb=1 {Maxumum growth rate of mutant population}
ma=1e-8 {Mutation rate A to B}
d/dt (B) = (vb*R/(kb+R))*B +GM  {Rate of growth of the} {mutant}
init B =0 {Initial mutant population size}
d/dt (A) = (va*R/(ka+R))*A - GM {Rate of growth of the} {ancestral population}
init A =1E3  {Initial size of the ancestral population}
d/dt (R) = - (vb*R/(kb+R))*B*eb -  (va*R/(ka+R))*A*ea  {Rate of change in concentration of the resource}
init R = 500 {initial resource concentration}
bm =  A*ma*(R/(R+ka))*DT {probabilty of a mutant B} {during interval DT}
rm =RANDOM (0, 1)
GM = IF rm < bm THEN PULSE (1,TIME,21) ELSE 0  {mutant addition condition}