how can I turn an ODE equation into a discrete probabilistic model?

I take for example the Verhulst equation for the growth of a population.

$$\frac{dP}{dt} = rP(1-P/K)$$

I was thinking to simulate the population dynamic simply as a loop (for the time) and at every time step, for each member of the population evaluate the probability of procreation (in which case I add 1 to the population counter) and the probability that a member has died, in which case I remove 1 from the population counter.

I don't know though how to use r and K to evaluate such probabilities. Can you give me some hints or link me to some good resource that explains how to do this?

In general, which is the strategy to convert an ODE (or system of ODEs) into it probabilistic equivalent???

Thanks a lot in advance!

  • $\begingroup$ What kind of distribution do you expect the model to have? Uniform? Normal? Other? $\endgroup$ – Paul Mar 28 '12 at 14:53
  • $\begingroup$ What are your constraints on $r$ and $K$? $\endgroup$ – Deathbreath Mar 28 '12 at 16:58
  • $\begingroup$ @Paul my idea is that each "animal" has its own reproduction and death rates that are normally distributed and . However initially I can assume that the reproduction rate is constant. $\endgroup$ – lucacerone Mar 28 '12 at 19:38
  • $\begingroup$ @Deathbreath as I told Paul, I'd like to assume the rates are normally distributed $\endgroup$ – lucacerone Mar 28 '12 at 19:44

If I understand you correctly, you're may be looking for some implementation of the Gillespie Algorithm, which takes advantage of the idea that rates and probabilities can be converted into each other. The algorithm works by ordering the system in terms of event times, which seems to be pretty close to what you want to do.

Keeling and Rohani's book, while focused on mathematical epidemiology, has a decent treatment of several different implementations, and the trade-off between accuracy and computation time between the methods. It shouldn't be terribly difficult to adapt a epidemic model to a population growth model - they work in largely the same fashion.

This question (by me) on this site might also be of interest. I know there are also ways to do this in R, and the implementation looks like it might be fairly approachable to code by hand if you're so inclined.

  • 1
    $\begingroup$ The Gillespie (or SSA, or Kinetic Monte Carlo) algorithm is relatively straightforward to code correctly; the description given in Gillespie's original paper is very clear, so when I was told to implement Gillespie's algorithm for a homework assignment, I followed Gillespie's paper rather than the description in my class notes. $\endgroup$ – Geoff Oxberry Mar 29 '12 at 19:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.