Using the stochastic simulation algorithm discussed in class:
Q1. Implement the discussed simulation of death at a constant rate, starting with I0 = 10 individuals, with death rate (per individual per unit of time) k = 3d−1 (d = day), and plot the distribution (for an appropriate number of replicate simulations) of the time at which the population is fully extinct.
Q2. Implement the simulation of the emergence of a neutral mutant in an initially small isogenic population growing to saturation (system discussed in class, directly adapted from the work of Luria and Delbrück). We will take r = 1, μ = 5 × 10−7, N = 107 (number of individuals in the final population). Initial conditions will always be 0 mutants and 1 wild-type. The interesting output of the model is the number of mutants in the final population. Plot the distribution of this number of mutants in the final population, and compare with the theoretical distribution proposed by Ma, Sandri and Sarkar. If your simulations are too slow with the suggested parameters, you can change them to lower values of N.