Stochastic modeling of molecular reaction networks
1 / 44

Stochastic modeling of molecular reaction networks - PowerPoint PPT Presentation

  • Uploaded on
  • Presentation posted in: General

Stochastic modeling of molecular reaction networks. Daniel Forger University of Michigan. Let’s begin with a simple genetic network. We can list the basic reaction rates and stochiometry. numsites = total # of sites on a gene, G = # sites bound

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.

Download Presentation

Stochastic modeling of molecular reaction networks

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript

Stochastic modeling of molecular reaction networks

Daniel Forger

University of Michigan

Let’s begin with a simple genetic network

We can list the basic reaction rates and stochiometry

numsites = total # of sites on a gene, G = # sites bound

M = mRNA, Po = unmodified protein, Pt = modified protein

Transcription trans or 0+M

Translation tl*M+Po

Protein Modification conv*Po-Po, +Pt

M degradation degM*M-M

Po degradationdegPo*Po-Po

Pt degradationdegPt*Pt-Pt

Binding to DNAbin(numsites - G)*Pt -Pt, +G

Unbinding to DNAunbin*G-G

We normally track concentrationLet’s track # molecules instead

  • Let M, Po, Pt be # molecules

  • First order rate constants (tl, unbin, conv, degM, degPo and degPt) have units 1/time and stay constant

  • Zero order rate constant (trans) has units conc/time, so multiply it by volume

  • 2nd order rate constant (bin) has units 1/(conc*time), so divide it by volume

numsites = total # of sites on a gene, G = # sites bound

M = mRNA, Po = unmodified protein, Pt = modified protein

V = Volume

Transcription trans*V or 0 +M

Translation tl*M +Po

Protein Modification conv*Po -Po, +Pt

M degradation degM*M -M

Po degradationdegPo*Po -Po

Pt degradationdegPt*Pt -Pt

Binding to DNAbin/V(numsites - G)*Pt -Pt, +G

Unbinding to DNAunbin*G -G

How would you simulate this?

  • Choose which reaction happens next

    • Find next reaction

    • Update species by stochiometry of next reaction

    • Find time to this next reaction

How to find the next reaction

  • Choose randomly based on their reaction rates








bin/V(numsites - G)*Pt

Random #

Now that we know the next reaction modifies the protein

  • Po = Po - 1

  • Pt = Pt + 1

  • How much time has elapsed

    • a0 = sum of reaction rates

    • r0 = random # between 0 and 1

This method goes by many names

  • Computational Biologists typically call this the Gillespie Method

    • Gillespie also has another method

  • Material Scientists typically call this Kinetic Monte Carlo

Myth 1:“Mass Action Formulations do not account for Stochasticity”

Consider a simple model inspired by the circadian clock in Cyanobacteria







  • Here a protein can be in 3 states, A, B or C

  • We start the system with 100 molecules of A

  • Assume all rates are 1, and that reactions occur without randomness (it takes one time unit to go from A to B, etc.)

Mass Action Representation

Matlab simulation

Mass Action represents a limiting case of Stochastics

  • Mass action and stochastic simulations should agree when certain “limits” are obtained

  • Mass action typically represents the expected concentrations of chemical species (more later)

Myth 2:Stochastic and Mass Action Approaches agree only if there are enough molecules

What matters is the number of reactions

  • This is particularly important for reversible reactions

  • By the central limit theorem, fluctuations dissapear like n-1/2

  • There are almost always a very limited number of genes,

    • Ok if fast binding and unbinding

There are several representations in between Mass Action and Gillespie

  • Chemical Langevin Equations

  • Master Equations

  • Fokker-Planck

  • Moment descriptions

We will illustrate this with an exampleKepler and Elston Biophysical Journal 81:3116

Master Equations describe how the probability of being in each state

Sometimes we can solve for the mean and variance

Distribution of molecules often looks Gaussian

Moment Descriptions

  • Gaussian Random Variables are fully characterized by their mean and standard deviation

  • We can write down odes for the mean and standard deviation of each variable

  • However, for bimolecular reactions, we need to know the correlations between variables (potentially N2)

Towards Fokker Planck

  • Let’s divide the master equation by the mean m*.

  • Although this equation described many states, we can smooth the states to make a probability distribution function


If 1/m* is small, we can then derive a simplifed

Version of the Master equations

Chemical Langevin Equations

  • If we don’t want the whole probability distribution, we can sometimes derive a stochastic differential equation to generate a sample

Adalsteinsson et al. BMC Bioinformatics 5:24


  • Transcription Control

  • Lac Operon

  • Oscillations

  • Accounting for diffusion

Rossi et al. Molecular Cell

Ozbudak et al. Nature 427:737

Guantes and Poyatos PLoS Computational Biology 2:e30

Saddle-Node on an

Invariant Circle









SNIC Bifurcation

Invariant Circle

Limit Cycle








Hopf Bifurcation


limit cycle

Noise Induced oscillations

Liu et al. Cell 129:605

3-D Gillespie

  • Login