On the simulation of biological diffusion processes

Many phenomena of interest in biology can be modeled using diffusion processes satisfying a stochastic differential equation. We consider first the stochastic differential equation dX = mu Xdt + sigma X, dW, where W is a standard Wiener process and representing a population growth process. This is simulated using both a strong Euler scheme involving normal pseudorandom numbers and a weak Euler scheme using Bernoulli pseudorandom numbers. Results are given for the mean of X(1) and its 95% confidence intervals for various numbers of simulations. It is found that there are no significant differences between the results obtained by these two schemes at a particular value of the time step, but the weak scheme takes less computer time than the strong scheme. We also consider the process satisfying dX = (-gamma 1X + gamma 2 (1-X))dt + square root of X(1-X)dW, representing a gene frequency under the influence of random mating and mutation. It is similarly found that the results of simulation by the two schemes are not significantly different. It is concluded that in simulations of many biological diffusion processes it is often advantageous to employ a scheme involving Bernoulli rather than Gaussian random variates not only because it involves fewer machine arithmetic operations but also because problems with large jumps that sometimes occur with extreme values of normal variates are less likely, thus enabling one to employ a larger time step with a concomitant saving in machine time.