No Pugs

they're evil

science

bacterial colony mutation simulation

Below is a quick ruby program I threw together to try and make a very rough simulation of the impact of rare and barely beneficial mutations in bacteria colonies.

The way this code is here, there’s a strong type of bacteria (with the beneficial mutation) and a weak kind of bacteria (without the beneficial mutation)

The way these numbers are currently set, the bacterial colony starts at 1 weak bacteria and grows from there. There’s a 1 in 10 million chance that a weak bacteria will gain the strong mutation or that a strong bacteria will lose it’s strong mutation. The benefit to the bacteria is a 0.1% improvement in survival rate. Anytime the colony grows above capacity, a bunch are killed off (the current survival rate is 50%, so half die every time it grows above 4 million bacteria.)

Feel free to play with the numbers and see the various outcomes.

A quick explanation of the algorithm: I don’t test every single bacteria to see if it survives. Instead, if 49.9% of them are supposed to die, and let’s say that there’s 5 bacteria. 5 * 49.9 = 2.495. So I kill 2 bacteria and then generate a random number between 0 and 1 to test against 0.495 to see if the 3rd bacteria dies or not. Not testing every single bacteria directly makes it less like nature of course, but it also makes the algorithm much faster.

MUTATION_RATE = 1.0 / 10_000_000
CAPACITY = 4_000_000

BENEFIT = 0.001
SURVIVAL_RATE = 0.5

class Colony
  attr_accessor :survival_rate, :count
  def initialize rate, count = 0
    self.survival_rate = rate
    self.count = count
  end

  def double
    self.count *= 2
  end

  def kill
    self.count = self.count.cut(survival_rate)
  end
end

class Experiment
  attr_accessor :colony_w, :colony_s, :milestones, :capacity, :mutation_rate,
    :generations, :min_w, :min_s, :max_w, :max_s, :first_mutation

  def initialize survival_rate, benefit, capacity, mutation_rate
    self.capacity = capacity
    self.mutation_rate = mutation_rate
    self.colony_w = Colony.new(survival_rate, 1)
    self.colony_s = Colony.new(survival_rate + benefit)
    self.milestones = (0..10).map {|i| i * 0.1}
    self.generations = 0
    self.min_w = colony_w.count
    self.min_s = colony_s.count
    self.max_w = 0
    self.max_s = 0
    update_min_max
  end

  def update_min_max
    self.min_w = [min_w, colony_w.count].min
    self.min_s = [min_s, colony_s.count].min
    self.max_w = [max_w, colony_w.count].max
    self.max_s = [max_s, colony_s.count].max
  end

  def total
    colony_w.count + colony_s.count
  end

  def percent_w
    colony_w.count.to_f / total
  end

  def percent_s
    colony_s.count.to_f / total
  end

  def double
    [colony_w, colony_s].each {|c| c.double}
  end

  def above_capacity?
    total > capacity
  end

  def tick
    run
    update_min_max
    check_milestone
    if generations % 10_000 == 0
      print_status
    end
  end

  def check_milestone
    if percent_s > milestones[0]
      stone = milestones.shift
      puts "colony_s has reached #{stone * 100}% of the population at generation #{generations}.
      total population: #{total}
      "
    end
  end

  def print_status
    puts "gen #{generations}: weak: #{colony_w.count} strong: #{colony_s.count}"
  end
  
  def mutate
    cols = [colony_s,colony_w]
    [cols, cols.reverse].each do |cs|
      mutated = cs[0].count.cut(mutation_rate)
      if mutated > 0
        cs[1].count += mutated
        cs[0].count -= mutated
        
        unless first_mutation
          puts "first mutation arises at gen #{generations}"
          print_status
          self.first_mutation = true
        end
      end
    end
  end

  def run
    double
    mutate
    while above_capacity?
      colony_w.kill
      colony_s.kill
    end
    self.generations += 1
  end
end

Integer.class_eval do
  def cut prob
    full = self * prob
    whole = full.to_i
    partial = full - whole

    rand < partial ? whole + 1 : whole
  end
end

e = Experiment.new SURVIVAL_RATE, BENEFIT, CAPACITY, MUTATION_RATE

e.print_status

while e.percent_w > 0.001 && e.generations <= 500_000_000
  e.tick
end

puts "took #{e.generations} generations to get population_strong from min of #{e.min_s}
to #{e.colony_s.count} (#{e.percent_s * 100}%)
and populations_weak from #{e.min_w} to a max of #{e.max_w} down
to #{e.colony_w.count} (#{e.percent_w * 100})%"

Here’s some of the results I get when running it with different values:

Using the current values

gen 0: weak: 1 strong: 0
first mutation arises at gen 23
gen 23: weak: 4194303 strong: 1
colony_s has reached 0.0% of the population at generation 29.
      total population: 2097152

colony_s has reached 10.0% of the population at generation 3911.
      total population: 2329272

colony_s has reached 20.0% of the population at generation 4317.
      total population: 2620602

colony_s has reached 30.0% of the population at generation 4587.
      total population: 2995564

colony_s has reached 40.0% of the population at generation 4808.
      total population: 3494804

colony_s has reached 50.0% of the population at generation 5010.
      total population: 2097183

colony_s has reached 60.0% of the population at generation 5213.
      total population: 2621926

colony_s has reached 70.0% of the population at generation 5434.
      total population: 3495609

colony_s has reached 80.0% of the population at generation 5703.
      total population: 2622864

colony_s has reached 90.0% of the population at generation 6108.
      total population: 2623759

took 7307 generations to get population_strong from min of 0
to 3258961 (99.0003848879678%)
and populations_weak from 1 to a max of 2097152 down
to 32906 (0.999615112032169)%

So with a benefit of 0.001, it took 29 generations for the mutation to surface, 5010 generations for it to match the frequency of the weaker bacteria, and 7307 generations for the strong bacteria to make up over 99% of the colony.

With a 0.0001 benefit (strong bacteria has a %0.01 better chance of survival over the weak bacteria)

gen 0: weak: 1 strong: 0
first mutation arises at gen 20
gen 20: weak: 2097151 strong: 1
colony_s has reached 0.0% of the population at generation 21.
      total population: 2097152

gen 10000: weak: 2095041 strong: 6691
gen 20000: weak: 2093070 strong: 56031
colony_s has reached 10.0% of the population at generation 27047.
      total population: 2324111

gen 30000: weak: 2091153 strong: 420262
colony_s has reached 20.0% of the population at generation 31090.
      total population: 2613776

colony_s has reached 30.0% of the population at generation 33781.
      total population: 2986635

colony_s has reached 40.0% of the population at generation 35989.
      total population: 3484132

colony_s has reached 50.0% of the population at generation 38015.
      total population: 2090424

gen 40000: weak: 1045189 strong: 1554583
colony_s has reached 60.0% of the population at generation 40043.
      total population: 2613202

colony_s has reached 70.0% of the population at generation 42254.
      total population: 3485008

colony_s has reached 80.0% of the population at generation 44952.
      total population: 2615286

colony_s has reached 90.0% of the population at generation 49018.
      total population: 2620706

gen 50000: weak: 262319 strong: 2870166
gen 60000: weak: 33896 strong: 2649526
gen 70000: weak: 5297 strong: 2445672
took 75938 generations to get population_strong from min of 0
to 2004344 (99.9000174446134%)
and populations_weak from 1 to a max of 2097151 down
to 2006 (0.0999825553866474)%

It took 38015 generations for the strong and weak bacteria to be equal in numbers, 75938 generations for the strong bacteria to go from 0 to over 99% of the population. Almost 10 times longer, which is interesting because the benefit was decreased by 10.

With 0 benefit

Note: I changed the code to only print out the status every 50 million generations for this one, and only to run a maximum of 500 million generations

gen 0: weak: 1 strong: 0
first mutation arises at gen 25
gen 25: weak: 4194303 strong: 1
colony_s has reached 0.0% of the population at generation 27.
      total population: 2097151

colony_s has reached 10.0% of the population at generation 1117956.
      total population: 2097329

colony_s has reached 20.0% of the population at generation 2550113.
      total population: 2097589

colony_s has reached 30.0% of the population at generation 4577301.
      total population: 2097995

colony_s has reached 40.0% of the population at generation 8054126.
      total population: 2097269

colony_s has reached 50.0% of the population at generation 35885018.
      total population: 2101332

gen 50000000: weak: 1052661 strong: 1052518
gen 100000000: weak: 1051453 strong: 1052153
gen 150000000: weak: 1051255 strong: 1050806
gen 200000000: weak: 1052545 strong: 1051996
gen 250000000: weak: 1057002 strong: 1055415
gen 300000000: weak: 1053651 strong: 1055580
gen 350000000: weak: 1057720 strong: 1056790
gen 400000000: weak: 1059992 strong: 1059863
gen 450000000: weak: 1060187 strong: 1059305
gen 500000000: weak: 1058117 strong: 1059018
took 500000001 generations to get population_strong from min of 0
to 1059018 (50.0212787564326%)
and populations_weak from 1 to a max of 2097152 down
to 1058117 (49.9787212435674)%

Looks like it works it’s way to 50% very slowly, and then stays there. This is interesting and isn’t what I expected. I had expected that it would stay at whatever frequency it was at when it reached capacity. I guess what is happening is that whichever type there is more of is going to have more bacteria mutate into the other type than it gets back in return.

Published on 03/20/2009 at 04:00PM under , .

Powered by Typo – Thème Frédéric de Villamil | Photo Glenn