I have found a new hobby: simulating evolutionary models. Mostly I do it in the train on my way to and from work. Yes, you can do a lot in 1 hour a day! And it’s quite captivating.
When I tell people about this, they ask me how I make simulations. Here’s how.
Yes, I program them myself. I use Octave, a free and open-source environment similar to Matlab. Octave supports a large part of the Matlab syntax and library. It also has its own extensions, such as C-like operators. I like to use ++ and +=, but today’s example will be Matlab-compatible.
The example simulates the survival of the fittest, and is only 40-odd lines long.
Here we go:
worldCapacity = 100; survivalScores = [1 1.2]; % of the 2 types initNum = worldCapacity; init2Ratio = 0.1; % The ratio of 2's initially meanNumChildren = 1; maxNumChildren = 3;
Our world has a limited capacity — more about that later. There will be two kinds of cells, type 1 and type 2. As you see, type 2 is slightly better at survival.
First, we create our population:
% init population = ; for i = 1 : initNum population(i) = (rand() < init2Ratio) + 1; end
rand() returns a random number between 0 and 1.
rand() < p is true with probability p.
rand(i, j) returns an i x j matrix.
As you see, initially, only 10% of the cells are type 2.
generation = 0; ratios = ; while 1 generation = generation + 1; % show ratios(generation) = length(find(population == 2)) / length(population); plot(1 : generation, ratios); drawnow if kbhit(1) break; end
kbhit() because Octave doesn’t always handle Ctrl-C well.
%reproduce children = ; for cell = population for i = 1 : maxNumChildren if rand() < meanNumChildren / maxNumChildren children(end + 1) = cell; end end end population = [population, children];
As you see, reproduction is asexual. With maxNumChildren large enough, the distribution becomes close to Gaussian, but that’s not today’s point.
%survive scores = survivalScores(population); prob = scores * (worldCapacity / length(population)); survive = rand(1, length(population)) < prob; population = population(survive); end
Our world has limited resources and can host only about
worldCapacity cells. So they have to compete for survival. There is no absolute survival in this model, i.e., if a cell is in the world alone, it won’t die. The reason a cell can die is that someone else is better at getting resources from the environment.
Most of the time I don’t use absolute survival at all. The reason to use competitive survival rather than (only) absolute survival is that with absolute survival, unless you calibrate the probabilities very well, you population will either die out or explode, making your simulation slow and eating up too much memory. The other reason of course is that in the real world, resources are indeed limited.
Competitive survival works like this: Each organism gets a survival score. They are then all multiplied by the same multiplier, such that their sum becomes
worldCapacity. The products are survival probabilities. Of course, if a “probability” is greater than 1, the organism surely survives. This may be a little unnaturalistic and also gives the population a bias towards an underpopulated world.
What next? The next step is to download the whole program and run the simulation!