In the paper 'Bottlenecks can constrain and channel evolutionary paths' by Jasmine Gamblin, Sylvain Gandon, François Blanquart and Amaury Lambert, we present a model for populations undergoing periodic bottlenecks and show how demographic parameters affect which evolutionary paths are accessible to these populations.
We used simulations to verify that our results derived for large population sizes and exponential growth between bottlenecks remain valid in slightly relaxed settings (i.e. finite population and/or density-dependent growth).
Here is the code that we used to perform these simulations. All functions are in the file sim_pop_functions.py, and the notebook sim_population_bottlenecks.ipynb provides an example of how to use these functions to perform simulations and plot results.
We simulated populations evolving according to a fully stochastic model: individuals are grouped by subpopulations carrying the same genotype, each following a birth-death-mutation process. Double mutants originating from the weakly beneficial single mutant or the strongly beneficial single mutant were grouped separately.
We used the principle of Gillespie's algorithm to simulate the evolution of the system during growth phases: starting from time 0, we draw an exponentially distributed random variable to fix the time increment (which is also the time of the next event). Then we pick randomly the event's identity (birth, death or mutation) proportionally to the corresponding rates. Then we update the population sizes according to this event. This procedure is repeated until the time reaches the length of one cycle. The bottleneck step is then performed by picking for each population a new size given by a binomial variable with parameters (population size, dilution factor).
We resorted to a tau-leaping approximation to speed-up the simulation when population sizes become large. Indeed, when the population size increases events happen more often and time intervals between events become smaller and smaller. Thus we fix a lower limit to the time increment (tau). When rates become too high, we use the fixed time increment tau and simulate Poisson random variables to determine how many of each events happened during that interval. If after updating all population sizes one of them is negative, then we divide the increment by 2 and redraw the Poisson variables.