Effective Biological Simulation Techniques

1
Effective Biological Simulation Techniques

• Kevin Burrage
• Federation Fellow of the ARC kb_at_maths.uq.edu.au
  
• UQ, Brisbane, Australia

2
Contents

• Introduction
• Biology and Noise
• Cellular Dynamics
• Modelling, Simulation and Visualisation
• Simulation of Chemical Reaction systems
• New simulation algorithms
• Numerical results
• Conclusions

3
Biology and Noise
4
Multiscale Chemically Reacting System
5
I. Biological Evidence of noise

• Stochasticity is evident in all biological
  processes the proliferation of both noise and
  noise reduction systems is a hallmark of
  organismal evolution Federoff et al.(2002).
• Transcription in higher eukaryotes occurs with a
  relatively low frequency in biologic time and is
  regulated in a probabilistic manner Hume
  (2000).
• Gene regulation is a noisy business Mcadams
  et al. (1999).
• Stochastic variation in protein levels is
  predominantly generated at the translational
  level (Ozbudak et al, Nature Genet. 2002, 69).

6

• i) Genetic model of Drosophila ii) Noise in
  eukaryotic gene expression Blake (2003)
  stochastic model for yeast. Model
  predicts linear correlation between
  noise strength and translational
  efficiency.
• Von Dassow 2001 Nature

7
II. Stochastic mechanisms

• Gene expression within a cell is a complex
  process chromatin remodelling, transcription,
  export of RNA, translation of mRNA into proteins.
• Physiological activity and cell differentiation
  within a mammalian cell is controlled by more
  than 10,000 protein coding genes.
• Thousand of genes are expressed at very low level
  copy numbers, which extant gene profiling methods
  cannot reliably detect.
• Initiation of gene transcription is a discrete
  process in which individual protein-coding genes
  in an off state can be stochastically switched
  on, resulting in sporadic pulses of mRNA
  production Sano 2001.

8
Noise can be classified as internal
stochastic fluctuations in number of proteins
external fluctuations in environment and
control parameters. Noise also classified as
phenotypic noise leading to
qualitative differences in a cell phenotype (eg.
lysis-lysogen pathway). Fluctuations cannot
always be viewed as simply small perturbations as
they can, in fact, induce different developmental
pathways. stablizable noise leading to
fluctuations in protein concentrations
(robustness properties of biological systems).
9
Cellular Dynamics
10

• Cells are very complex
• A wide variety of ultrastructure
• Many different types of dynamics processes and
  transports vesicles/tubules etc.
• 
  Three D EM image of a
  mammalian insulin
  secreting cell (Marsh)

11
Frame 10
Frame 20
Untreated cells have little tubulation.
Frame 0
Frame 30
Frame 40
Frame 50
Frame 60
Frame 70
Frame 80
12
Frame 0
Frame 10
Frame 20
Frame 30
Frame 40
Frame 50
Frame 60
Frame 70
Frame 90
Frame 80
Massive tubulation event following EGF treatment,
apparent reduction in the volume of the vacuolar
compartment.
13
Modelling, Simulation and Visualisation
14
III. Modelling Issues

• Even small models have 10s of free parameters
  half lives, binding/reaction rates.
• Real values may be unknownbiologically realistic
  ranges may span orders of magnitude and data
  may be at many different scales.
• Models must be robust under the above constraints
  and a wide range of initial conditions.
• Need to match known gene effects to model
  outcomes - e.g. nonlinear optimisation
  Drosophila model.
• Model prediction must be for biologically valid
  data sets.
• In some cases we can test the topology
  (structure) of a model rather than quantative
  tuning.
• Cells as spatial objects.

15
(No Transcript)
16
IV.Simulations - Grid Implementation

• Simulations of genetic regulatory models on a
  compute grid.
• Computationally intensive if modelling genetic
  regulation within thousands of cells.
• Easy to use by Biologists remote access.
• View the simulation outputs from within a
  browser.
• Parallelisation across the simulations.
• Message passing via MPI.
• Web interface via MathML interactive entry of
  input
• Coded in C.
• Statistics computed on the host.
• Mix of Java, Perl and Shell scripts.
• Graphs displayed within GNUplot.

17
(No Transcript)
18
V. Visualising Cellular Processes

• Need to couple simulations with visualisation and
  data
• Virtual Cell http//www.nrcam.uchc.edu/
• JigCell http//gnida.cs.vt.edu/cellcyclepse/
• BioSpice https//community.biospice.org/
• E-cell http//www.e-cell.org/
• GPL M-Cell http//www.mcell.cnl.salk.edu/
• Visual Cell http//www.acmc.uq.edu.au/cpst/1.html
  

19

• Visualisation of stochastic chemical reaction
  simulations in three D space Brownian diffusion
  (left), Vesicle transport
• (right) - Igor Kromin and Radosav Pantelic (UQ)

20
Simulation of Chemical Reaction Systems
21

• VI. The modelling regimes
• Discrete and stochastic - Finest scale for well
  stirred systems. Based on detailed biochemical
  reactions. Exact description via Stochastic
  Simulation Algorithm (SSA) - Gillespie. Large
  computational time.
• Continuous and stochastic - A bridge connecting
  discrete and continuous models. Valid under
  certain conditions. Described by SDEs The
  Chemical Langevin Equation.
• Continuous and deterministic - The Reaction Rate
  equations. Described by ordinary differential
  equations. Not valid if molecular populations of
  some critical reactant species are small, so that
  microscopic fluctuations along with reaction
  channel feedbacks produce macroscopic effects .

22
VII. Stochastic Simulation Algorithm

• Simulates the time evolution of a well stirred
  chemical reacting system by taking proper account
  of inherent randomness.
• Well-stirred mixture
• N molecular species
• Constant temperature, fixed volume
• M reaction channels
• Dynamical state
• where is the number of
  molecules in the system
• Propensity function the probability,
  given , that one reaction
  will occur somewhere inside in the next
  infinitesimal time interval

23

• When that reaction occurs, it changes the state.
  The amount by which changes is given by
• the change in the number of
  molecules produced by one reaction.
• Generate a tentative reaction time for each
  reaction channel according to
• where are M statistically
  independent samplings of U(0,1)
• Take
• Update

24
An example a system with three species of
molecules and four reactions
The four state transfer vectors are
If the current state is (x1,x2,x3)(100,100,100),
After reaction 2, (x1,x2,x3)(98, 101,100)
4, (x1,x2,x3)(98, 100,101)
25
VIII. Chemical Langevin Equation

• If the system possesses a macroscopically
  infinitesimal time scale, so that during any dt
  on that scale all of the reaction channels fire
  many more times than once yet none of the
  propensity functions change appreciably, we can
  approximate the discrete Markov process by a
  continuous Markov process defined by the Chemical
  Langevin (SDE) Equation
• where are M temporally
  uncorrelated, statistically independent normal
  variables with mean 0 and variance 1

26
New simulation algorithms
27
IX. Method improvement

• The SSA is very expensive must compute m
  reaction times and can be small
• Given a subinterval of length , if we can
  determine how many times each reaction channel
  fires in each subinterval, we can forego knowing
  the precise instants at which the firings took
  place. Thus we could leap from one subinterval
  to the next, rather than one reaction to the
  next.
• How long can that subinterval be? Tau-leaping is
  exact for constant propensity functions, thus
• LEAP CONDITION is selected so that no
  propensity function changes appreciably.
• If the reactant population is large can be
  moderate sized and this gives the Explicit Euler
  method in SDE and ODE regimes.

28
(i)Poisson leap method (Gillespie, JCP,
115(2001),1716) Assumption In the time period
, the number of reactions for reaction
channel is Poisson For the given criterion
, choose a stepsize to satisfy the leap
condition Generate Poisson numbers Update the
system
29
(ii)Binomial leap method (Tian and Burrage, JCP,
2004) Assumption In the time period ,
the number of reactions for reaction channel
is Binomial Let Generate the Binomial random
variables. Update the system Avoids the
possible negative numbers that can arise in the
Poisson leap methods.
30
(iii) The midpoint Poisson and Binomial leap
methods Use the midpoint for better
approximation Predict Update (iv) The
implicit tau-leap method (Rathinam et al.,) For
stiff systems (values of propensity functions
widely varying), small stepsize must be used for
fast reactions. Similar to the ODE case, implicit
tau-leap method is derived in order to improve
stability and to use larger stepsize.
31

• (v) A multi-scale approach - Burrage, Tian
  Burrage Progress in Biophysics and Molecular
  Biology, 2004
• Divide the system into
• Fast reactions large numbers of molecules, use
  the Langevin approach
• Intermediate reactions moderate numbers - use
  the tau-leap methods
• Slow reactions small numbers of molecules, use
  SSA.
• We have used this approach to simulate the
  expression and
• activity of LacZ and LacY proteins in E.coli .
• This problem has 22 reactions and takes several
  hours to do one simulation by the SSA.

32
Numerical Results
33

• Test problem 1
• Reversible dimerization of monomer
• into unstable , which can convert to a
• stable by reaction .
• can also be degraded through .
• Test problem 2
• Ecoli problem with 22 reactions LacZ and LacY
  genes

34

• Results
• For Problem with 4 reactions, PRK method speed up
  of factor of 7 over SSA.
• Tau leap inferior to the other three methods.
• For problem 2, The Binomial leap method is a
  factor 250 times faster than SSA. So significant
  improvements as the size of the problem scales.

35
(No Transcript)
36
Conclusions
37

• Poisson leap method gives substantial
  improvements over SSA.
• Binomial approach robust and even more effective
  than Poisson approach.
• Need to exploit a multiscale approach through 3
  regimes.
• How do we deal with spatial issues.?
• Cells are heterogeneous with complex
  ultrastructure so how do we extend these ideas
  beyond the simple Monte Carlo approach of picking
  a molecule, letting it move and seeing if a
  reaction occurs?

38
Co-workers and Thanks

• Tianhai Tian, Christine Beveridge,Pamela Burrage,
  Francis Clark, Jim Hanan, Thomas Huber,John
  Mattick, Lars Nielsen,Mark Ragan (UQ)

39

• References
• Burrage, K. and Burrage, P.M. (2003) Numerical
  methods for stochastic differential equations
  with applications, SIAM Monograph, to appear
• Gillespie, D. (1977) Exact stochastic simulation
  of coupled chemical reactions, J. Phys. Chem. 81,
  2340.
• Gillespie, D.T. (1992), Markov Processes an
  introduction for Physical Scientists, Academic
  Press.
• Gillespie, D. (2001) Approximate accelerated
  stochastic simulation of chemically reacting
  systems, J. Chem. Phys. 115, 1716-1733.
• Hasty, J. et al, Noise-based switches and
  amplifiers for gene expression, PNAS, 97 (2000),
  2075-2080.
• Hume, D. (2000) Probability in transcriptional
  regulation and its implications for leukocyte
  differentiation and inducible gene expression.
  Blood 96 2323-2328.

40

• H.H.McAdams et.al., It is a noisy business!
  Genetic regulation at the nanomolar scale, Trends
  Genet, 15 (1999), 65-69.
• H.H.McAdams et.al., Stochastic mechanisms in gene
  expression, PNAS, 94 (1997), 814-819.
• P.Smolen et.al., Mathematical modelling of gene
  networks, Neuron, 26 (2000), 567-580.
• N.Thattai et.al., Intrinsic noise in gene
  regulatory networks, PNAS, 98 (2001) 8614-8619.
• Elowitz, M.B., Levine, A.J., Siggia, E.D. and
  Swain, P.S. (2002) Stochastic gene expression in
  a single cell, Science 297, 1183-1186.
• Fedoroff, N. and Fontana, W. (2002) Small numbers
  of big molecules, Science 297, 1129-1131.

41

• Kuznetsov, V.A., Knott, G.D. and Bonner, R.F.
  (2002) General statistics of stochastic process
  of gene expression in eukaryotic cells, Genetics
  161, 1321-1332.
• Sano et al. (2001). Random monoallelic expression
  of 3 genes clustered within 60Kb of mouse
  complex genomic DNA, Genome Res 11, 1833-1845.
• Strogatz, S. Exploring complex networks, Nature
  410, March 8, 2001.
• Tian, T. and Burrage, K (2004) Binomial learp
  methods for simulating chemical kinetics,
  submitted to J. Chem Phys.
• Von Dassow, G. et al. (2000), The segment
  polarity network is a robust developmental
  module, Letters to Nature.