Today, we’ll go into the nuts and bolts of a general model that can be twisted to fit a wide variety of studies. Its basic principles are very simple, and the real mathematical difficulties start when trying to apply it at different scales, e.g. when trying to understand what happens when the population is huge or when the time of the study is very long. We’ll come around to that, but let’s first start by describing exactly what we’re aiming to do here.

So say hello again to our old friends, the heterotrophic bacteria!

Heterotrophic bacteria

Starting from scratch

Before we can try to understand how entire populations will evolve under climate change, we need to have a grasp at how individuals live their day-to-day lives. Once the behavior of all individuals is properly described at every moment, studying evolution will just be a matter of unfolding events until some sort of ‘endgame’ is reached. Simple, right?

You may think that there is a catch, and indeed, you’ll see some mathemical mumbo jumbo take place for when we’ll want to predict this endgame with just a pen and paper. But when simulation is concerned, the steps I’ve described are pretty much all we need to do in order to see evolution take place. Granted, this process is computationally expensive, so it only works for simple populations in only one location. Don’t expect to model an entire ocean worth of heterotrophic bacteria though! Still, it is a perfectly reasonable tool for precise studies of complex phenomena taking place in a single location.

So, what can happen to a bacterial cell during its short life? Well, in the most general sense, the only interesting things that can happen to it are either it divides, giving birth to a new bacterium, or it dies.

Birth and death process

Of course, these events depend on multiple things: for instance, the division rate will depend on the nutrient concentration, and the mortality rate on competition with the neighbors, or the presence of predators.

As a whole, this type of model is called a birth and death process (the name is self-explanatory I think). Since we’re trying to describe the whole life of the bacteria, we now need some sense of time. When do they give birth, and when do they die?

A first possibility would be to look at time as a series of discrete generations. Each time step corresponds to one generation, and each new generation is generated by looking at the number of bacterial cells that divided or died from the previous generation. But this modeling framework, useful as it may be, won’t be enough for us. Indeed, bacteria divide and die at very different times, and thus generations can very rapidly get mixed. That’s why we want to describe what happens at every moment: we need a model that is continuous in time, and not discrete.

But how do we do that? Well, the solution found to model the times at which the events occur is, in my opinion, very elegant and evocative. For each cell, we define ‘biological clocks’ associated with each event. Here, each individual will have a ‘division clock’ and a ‘death clock’, which are each set up to go off at a random time.

Birth and death clocks

This randomness is weighted with the expected rate at which the events occur though, so we can make predictions about what will happen. For instance, when the environment is rich in dissolved organic matter (the resource for heterotrohic bacteria), we know that the ‘division clock’ will probably go off sooner than in a poor environment, ensuring faster growth of the population1.

We then look at the first clock to go off among all the clocks of all the individual cells, and update the population accordingly. If the first clock to go off is a division clock, we can increase the population by one individual, and if it is a death clock, we decrease the population by one individual.

Here, the division clock rang before the death clock, resulting in an increasing population.

Here, the division clock rang before the death clock, resulting in an increasing population.

We then reset the clocks2 and start over until some sort of equilibrium is reached. If we then look at the size of the bacterial population at each instant, we can expect the evolution to look something like this:

During the studied lapse of time, there were two births and a death in the bacterial population.

During the studied lapse of time, there were two births and a death in the bacterial population.

Let’s now take a step back and look at a concrete simulation:

IBM Population

As you can see, we now know how to simulate an entire population of separate individuals, and track its size across time. The study of such processes is called population dynamics, but you might have noticed that something is missing for our particular interest: darwinian evolution, or natural selection.

Introducing natural selection

As for the length of girafes' necks in our last article, we’ll try to focus on a specific bacterial characteristic (or trait) that is variable, heritable and under selection pressure. There are a number of possible traits, each interesting for a specific study – so for today, we’ll focus on a generic trait that we’ll call $\tau$. We’ll assume that $\tau$ is a number (not necessarily a whole number) representing a unique characteristic: for instance, $\tau$ could be the average neck length of the giraffes of trait $\tau$, but for bacteria it could be the size, the maximum growth rate, etc… What is important is that this trait influences the chances of survival and/or reproduction for the individuals, so that two populations of different traits would have different life histories.

In general, when a bacterial cell divides, both daughter cells inherit the trait of their mother.

Normal division

But in some rare cases, a mutation appears at random. This will slightly alter one of the daughters' trait, which will no longer bare the trait $\tau$ but the slightly different trait $\tau'$. In effect, we can assume that each division has a very small probability of resulting in the apparition of a mutant.

Mutation division

Then, once the mutant appears, we can let our model run its course! Indeed, since we focused on a trait under natural selection, we know that different traits will result in different birth and death rates. In the most simple case, if the birth rate of a population of trait $\tau'$ is higher than the birth rate of population of trait $\tau$ without changing the death rate, we can expect population $\tau'$ to outcompete population $\tau$ resulting in a global shift in population trait, from $\tau$ to $\tau'$.

If you want to see natural selection in action, check out the following simulation:

IBM Evolution

In this simulation, the $x$-axis corresponds to time (as before), but the $y$-axis is no longer the global population size, but the traits that appear in the population. Then, the coloring of each dot represents the size of the population bearing the trait $y$ at time $x$, with blue representing low population densities and red high population densities.

So for instance, at the beginning of the simulation the traits found in the population were spread around the value $\tau_1$, but shifted towards $\tau^* $ by the end. You just witnessed natural selection! We divised a simulation in which trait $\tau$ was an abstract trait governing the division rate and competition between individuals, with the optimal value being $\tau = \tau^*$. We see here that the trait slowly evolves from $\tau=\tau_1$ to its optimal value over time.

Already we have devised a model that allows us to study complex systems and phenomena. Indeed, for instance such a simple system is enough to understand the general mechanism behind speciation for instance, as you can see in the following simulation3:

IBM Branching

Here, we see that over time the population can divide between two trait values – if they stray far enough apart, the two populations will be different enough to be classified in different species!


Now that we have a robust algorithm for simulating natural selection, we should be able to study whatever system we fancy, right? Well, as I told you, this simulation scheme is costly, and not much can be said as such beforehand. But we would like to be able to predict the outcome for different scenarios of simulations, and not just rely on long and costly simulations. To do so, we need to make further assumptions, which will allow us to derive theoretical results that could be useful for studying complex systems.

Therefore, in the next couple of articles I’ll show you how changing scales can allow us to have a glance at the bigger picture of evolution, and how it will be useful to study the microbial loop. Stay tuned!

If you want to go deeper

  • Nicolas Champagnat, Régis Ferrière, Sylvie Méléard, Unifying evolutionary dynamics: From individual stochastic processes to macroscopic models, Theoretical Population Biology, Volume 69, Issue 3, 2006, Pages 297-321, ISSN 0040-5809,

  • Sylvie Méléard, Modèles aléatoires en Ecologie et Evolution, Springer Berlin Heidelberg, 2016,

  1. For all mathematicians out there, the random distribution used for our clocks is the exponential distribution with parameter the rate of the event. This distribution is the only that assures us that the process is Markovian, which is a very desirable property to have. ↩︎

  2. Mathematically speaking, this means ‘doing nothing’ since we chose the distribution to keep the process markovian, but from a simulation point of view this is an important step destined to make sure everything is mathematically sound. ↩︎

  3. You may notice that this simulation seems less ‘precise’ than the previous one: this is actually due to a lower population density, for which random events can upset the population much more than in larger populations. Parametrizing the model in order to be in the configuration needed for a speciation event to occur requires more computing power for the same population size, so in the interest of time (and to save my personnal computer from melting…), I chose to decrease the population size and thus the beauty of the graph. Sorry! ↩︎