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.

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