Simulating a Bernoulli experiment with a very large numbers of trials. each trial represents a single bacteria mutating or not mutating I am helping a friend of mine in making a simulation in MatLab involving the mutation of bacteria. First of all I would like to apologize in advance for my lack of knowledge of statistics. 
There is a starting amount of $10^8$ bacteria that grow exponentially over approximately $14400$ time steps. Now for each time step every single bacteria has a probability of $p \approx 1.2\cdot10^{-5}$ to mutate into a different type of bacteria. Now the way we designed the program is that we generate $10^8$ random number between $0$ and $1$ for the first time step and check how many are smaller or equal to $p$ and use that as the amount of bacteria that mutated. We then repeat this process the next time step. The problem however is that this is going to take a very long time.
The problem is that it is important that every single bacteria should have a chance to mutate for every time step. Is there any way we can find the number of mutating bacteria, or the number of successful  Bernoulli trials, for each time step in a faster way. Thanks in advance for any help!
 A: What you've described sounds similar to a pure birth process. I assume that you know that such a process has an analytical solution, the negative binomial distribution (see here). The additional aspect of mutation complicates it, however.
$10^8$ and 14,400 together are big numbers. There's no way around it, even implemented well, a simulation like this will take a long time on a desktop machine (GPU techniques may or may not help, depending on how you break up the problem). You are likely simulating in linear time, i.e., your time steps are equally spaced? Because this is an exponential growth process, it turns out that that is very inefficient. You should look into Gillespie's algorithm or exact stochastic simulation algorithms (SSA) and it's variants. This method is used for many things, including the simulation of discrete molecules in chemical reactions. The idea is to increment time randomly according to the probability that an event (growth and/or mutation in your case) will occur. It may also help if you think of your starting population as $10^8$ distinct trajectories. Unless you supply constraints, each initial bacteria is independent of all of the others.
Though it's applied to chemical reactions' you might find Modeling and Simulating Chemical Reactions by Desmond J. Higham helpful. There is Matlab code at the end.
Also, Matlab's SimBiology Toolbox (my academic license includes this toolbox) includes SSA and Tau-Leaping solvers. Though all the examples given seem to be from chemistry, these should work for a problem like yours (likely not as fast as a good hand-coded implementation). You might even contact The MathWorks to see if they have suggestions for using these tools for a problem like yours.
