May 11 - 13, 2008
The dynamics of chemical reaction networks can be modeled either deterministically or stochastically. The deficiency zero theorem for deterministically modeled systems gives conditions under which a unique equilibrium value with strictly positive components exists within each stoichiometric compatibility class (linear subset of Euclidean space in which trajectories are bounded). The conditions of the theorem actually imply the stronger result that there exist concentrations for which the network is "complex balanced." That observation in turn implies that the standard stochastic model for the reaction network has a product form stationary distribution.
Advances in the development of models that can satisfactorily describe biochemical networks are extremely valuable for understanding life processes. In order to get full description of such networks, one has to solve the inverse problem, that is, estimate unknowns (rates and populations of various species) or choose models from a set of hypothesized models using experimental data. In this work we discuss signal processing techniques for resolving the inverse problem of biochemical networks using a stochastic approach based on Bayesian theory. The proposed methods are tested in simple scenarios and the results are promising and suggest application of these methods to more complex networks.
In deterministic as well as stochastic models of chemical kinetics,
stiff systems, i.e. systems with vastly different time scales where the
fast scales are stable, are very common. It is well known that
the implicit Euler method is well suited for stiff deterministic
equations (modeled by ODEs) while the explicit Euler is not. In particular
once the fast transients are over, the implicit Euler allows for the choice of times steps comparable to the slowest time scale of the system.
In stochastic systems (modeled by SDEs) the picture is more complex.
While the implicit Euler has better stability properties over the
explicit Euler, it underestimates the stationary variance. In general one
may not expect any method to work successfully by taking time steps of the
order of the slowest time scale. We explore the idea of interlacing large
implicit Euler steps with a sequence of small explicit Euler steps.
In particular we present our study of a linear test system of SDEs and
demonstrate that such interlacing could effectively deal with
stiffness. We also show uniform convergence of mean and variance.
Local and global stability of biochemical reaction network
Modern biological research provides countless examples of
biochemical interaction networks. For example, at the
level, the nodes of these interaction networks could be
molecules, genes, and gene products. In order to
role played by some of these interactions one often faces
difficulties in trying to interpret the effect of positive
negative feedbacks, nonlinear interactions, and other
signaling between the nodes of the network. We will
connections between reaction network structure and the
complex dynamic behavior. These connections may play an
role in facilitating the understanding of experimental
Reaction-Diffusion models are key components of models in developmental biology.
These reaction diffusion processes can be mathematically modelled using either deterministic partial differential equations or
stochastic simulation algorithms. Here we discuss the stochastic simulations on both linear and non-linear Reaction-diffusion models
The stochastic simulation algorithm, or Gillespie Algorithm, is a tool used to simulate discrete biochemical reaction systems when there are a small to moderate number of molecules. The Gillespie Algorithm has been adapted to handle systems with delays (such as with gene transcription and translation) and approximations (such as tau leaping) have been developed in order to increase computation speed. In this poster we will model discrete biochemical systems via a random
time change representation. We will then demonstrate how natural and
efficient modifications of the stochastic simulation algorithm arise
from such a representation. Also, we will use this representation to
show how to incorporate post-leap checks in the tau-leaping algorithm.
Biochemical systems can often be viewed as discrete-event systems, i.e., as systems that make stochastic state transitions at a strictly increasing sequence of random times. We survey a number of topics pertinent to modeling and simulation of such systems. We first describe several basic models for discrete-event systems, such as generalized semi-Markov processes, stochastic Petri nets, and continuous time Markov chains, and discuss the interplay between the choice of modeling formalism, the compactness of the model representation, and the computational complexity of the resulting simulation algorithm. We then outline a collection of techniques for increasing the efficiency of a simulation, as well as for efficiently estimating the sensitivity of a discrete-event system model with respect to one or more model parameters.
We analyze a heat shock response model developed by Srivastava, Peterson, and Bently (2001) to find appropriate scaling for the abundance of molecules of chemical species and for chemical reaction rates. Extending the method introduced in the paper by Ball, Kurtz, Popovic, and Rempala (2006), we apply the multiple scaling method to the various range of the abundance of chemical species and the chemical reaction rate constants.
In this method, the complicated original system is split into several subsystems with separate fast, intermediate, and slow processes by using the approximated scaled processes with appropriate time change. Then, we can approximate the state of the system from the state of the reduced subsystems. In the slow time scale, we identify behavior of fast processes since they already start to move while the slow species remain mostly constant. In the fast time scale, we identify behavior of the slow processes driven by averaged behavior of fast processes.
Our goal is to find an appropriate scaling that makes the abundance of molecules of each chemical species well-balanced. To make each chemical species balanced, the maximal order of magnitude of the reaction rates of production for each chemical species must be the same as the maximal order of magnitude of the reaction rates of consumption. We also need a balance condition between the maximal order of magnitude of collective reaction rates of inflow and outflow in each atom graph. In case any of balance conditions are not satisfied, we put restriction on the range of time change exponent. In this poster, I will introduce multiscale method briefly and give a simulation result for one set of
exponents meeting our balance conditions.
In the presented work we will describe how to rationalize synthetic biology
using model-driven, molecular-level engineering principles. In the poster
we will focus on the theoretical effort to develop an algorithm for
simulating biomolecular systems across all relevant time and length scales;
from stochastic-discrete to stochastic-continuous and
deterministic-continuous models, we are developing the theoretical
foundation for accurately simulating all biomolecular interactions in
transcription, translation, regulation and induction and how these result in
phenotypic probability distributions at the population level.
The nascent field of synthetic biology offers the promise of engineered gene
networks with novel biological phenotypes. Numerous synthetic gene circuits
have been created in the past decade, including bistable switches,
oscillators, and logic gates. Despite a booming field and although recently
developed designs of regulatable gene networks are ingenious, there are
limitations in routinely engineering synthetic biological systems. Indeed,
there is a need for rationalizing the design of novel regulatable gene
networks that can be used in useful applications.
We are developing multiscale mathematical tools to assist synthetic biology
efforts. Why are multiscale models necessary to assist synthetic biology,
and not simply apply the mathematics developed to model kinetic and
thermodynamic processes in living organisms? Because, although the
principles of thermodynamics, kinetics and transport phenomena apply to
biological systems, these systems differ from industrial-scale chemical
systems in an important, fundamental way: they are occasionally far from the
thermodynamic limit. Indeed, using ordinary differential equations for
simulating the reaction kinetics of these systems can be distinctly false.
The need arises then for stochastic models that account for inherent,
thermal noise, which is manifest as phenotypic distributions at the
population growth/interaction levels.
In the presentation, we will discuss synthetic biological systems, the
development of multiscale models to capture the phenotypic behavior of these
systems and we will present examples of model-driven synthetic
We have investigated how noise propagation in biological reaction networks affects system sensitivities. We have shown that the sensitivities can be enhanced in one region of system parameter values by reducing sensitivities in another region. We have applied this sensitivity compensation effect to enhance the amplification of a concentration detector designed by using an incoherent feed-forward reaction network. We have also shown that metabolic control analysis can be extended for stochastic reaction systems by providing new versions of the summation and connectivity theorems.
In recent years, single-molecule experiments have emerged as a powerful
means to study biophysical/biochemical processes; many new insights are
obtained from this single-molecule perspective. One phenomenon recently
observed in single-molecule biophysics experiments is subdiffusion, which
largely departs from the classical Brownian diffusion theory. In this
talk, by introducing fractional Gaussian noise into the generalized
Langevin equation, we propose a model to describe the subdiffusion. We
will study the analytical properties of the model, compare the model
predictions with experiments, look at its connection with reaction
networks, and explore the implications of the model on enzyme reactions.
Classical stochastic models for chemical reaction networks are
given by continuous time Markov chains. Methods for characterizing these
models will be reviewed focusing primarily on obtaining the models as
solutions of stochastic equations. The relationship between these equations
and standard simulation methods will be described. The primary focus of the
talk will be on employing stochastic analytic methods for these equations to
understand the multiscale nature of complex networks and to exploit the
multiscale properties to simplify the network models.
Engineering intracellular metabolism has a long history, dated back to the late 1950s. The field has enjoyed significant attention in recent years, particularly because of the energy and environmental concerns. Many success stories have highlighted the feasibility of engineering intracellular metabolism to achieve a designed function. In almost all cases, the bottleneck of the problem resides on the biological side. With some exceptions, mathematical modeling has rarely played a predictive role in this development, mainly because of the lack of appropriate model structure and useful parameters. This talk will discuss a few examples of metabolic engineering and examine the role of mathematical modeling in these developments. In particular, we will present an en ensemble modeling technique that aims to narrow the gap been experiment and model prediction.
Actin polymerization and network formation are key processes in cell
motility. Numerous actin binding proteins that control the dynamic
properties of actin networks have been studied and models such as the
dendritic nucleation scheme have been proposed for the functional
integration of at least a minimal set of such regulatory
proteins. In this talk we will describe recent work on
the evolution of the distribution of filament lengths and nucleotide
profiles of actin filaments from both the deterministic and
stochastic viewpoints. For the latter we develop a master
equation for the biochemical processes involved at the individual
filament level and simulate the dynamics by generating numerical
realizations using a Monte Carlo scheme. A new computational algorithm
that is far more efficient than standard methods will also be
Living cells contain such low numbers of active genes, RNAs
and proteins that random fluctuations in concentrations arise
spontaneously. Because many biological processes require reliability and
precision, cells are thought to deal with this problem by using negative
feedback control loops, correcting perturbations by increasing synthesis
at high concentrations and decreasing it at low concentrations. Negative
feedback has been systematically studied in control theory, but focus
has been on macroscopically large systems subject to external
perturbations - like airplanes in random gusts of wind - while other
principles arise in discrete molecular-scale probabilistic processes.
Many biological control loops are also strongly nonlinear, and crucial
aspects are often unknown. Such systems are typically considered
impossible to analyze: how can we make quantitative statements about
strongly nonlinear stochastic systems if we do not even know the nature
of the nonlinearities? By combining exact analytical tools from
functional analysis, information theory and statistical physics, I will
show that seemingly mild constraints, such as short delays or finite
numbers of control molecules, can place fundamental limits on noise
suppression that no arbitrarily elaborate feedback system could ever
overcome - hard physical bounds that can be derived explicitly. Using
approximate methods for families of feedback systems, I will also show
how the limits arise because different types of noise are suppressed
according to incompatible principles, generating frustration trade-offs
where reducing one type of variation inevitably amplifies another.
Finally I discuss partial loopholes in the general laws where
counterintuitive mechanisms can circumvent the trade-offs to some
extent. The general results are illustrated by bacterial plasmids, where
large fluctuations promote extinction and where numerous mechanisms have
evolved to approach the physical limits. I will emphasize intuitive
reasoning and physical observables throughout, and show some preliminary
Several deterministic models will be presented to illustrate
the special types of mechanisms that regulate and control cell metabolism.
Stochastic analysis can be used to investigate such highly nonlinear
deterministic dynamical systems that typically operate far from
equilibrium. Examples will be given that show how specific biological
system properties depend on the natural stochasticity inherent in
diffusion processes and in the making and breaking of chemical bonds. On
many different space and time scales, evolution has selected
stochastic mechanisms for accomplishing specific biological functions.
Finding, characterizing, and explaining these mechanisms requires
cooperative research by biologists and stochastic analysts. The tasks
are daunting because biology is really hard.
Consider the system of reaction rate equations (RRE) describing a chemical network with the reaction rate constants considered to be unknown parameters. The talk shall describe a statistical approach to identifying the "most likely network" from a given set of RRE coefficient estimates. The idea relies on mapping the estimated reaction constants into an appropriate convex region of a network stochiometric space in order to identify the reactions which are most likely to span that region. This approach reduces the original problem to inferring parameters of a certain multinomial distribution which may be solved using the general methods of algebraic statistics.
The reliable genetic engineering of bacterial systems would be
greatly aided by a more quantitative and predictive understanding of
gene expression. In many synthetic biology applications, the
production rate of a protein needs to be precisely tuned to optimize
some performance characteristic of the genetic system. A typical
application might require selecting the optimal expression level of
enzymes in a metabolic pathway. The production rate of a protein can
be precisely tuned by varying the mRNA sequence of the 5' UTR
(untranslated region), including the ribosome binding site (RBS), that
precedes the protein's coding sequence. However, there is currently no
quantitative model that can predict an RNA sequence that yields a
desired production rate of protein.
We have developed and experimentally validated a statistical
thermodynamic model of translation initiation in Escherichia coli. The
model quantifies the effects of the important RNA, RNA-RNA, and
RNA-protein interactions and determines the resulting translation
initiation rate of a specified mRNA transcript. We combine the
quantitative model with Monte Carlo sampling to create a design method
that generates synthetic 5' UTRs. The user inputs a protein coding
sequence and a desired translation initiation rate and the design
method will generate the corresponding RNA sequence. We measure the
accuracy of the design method by generating numerous synthetic 5' UTRs
that drive the expression of the fluorescent protein RFP in a
simplified test system. We then compare the design method's
predictions with flow cytometry data from a test system. The design
method is capable of accurately generating RNA sequences that yield
the desired translation rate, with rates varying across three orders
The presented model and method is a key foundational technology that
will enable us to more rationally engineer synthetic biological
systems. With the emergence of lower cost large-scale DNA synthesis
and the rapid assembly of genetic systems, the development bottleneck
now shifts towards identifying the DNA sequence that ultimately yields
a desired system behavior. These and other quantitative and predictive
models aim to provide that much-needed missing link.
Although viruses are the smallest organisms with the shortest genomes, they have major impacts on human health, causing deadly diseases (e.g., influenza, AIDS, cancer) on a global scale. To reproduce, a virus particle must infect a living cell and divert biosynthetic resources toward production of virus components. A better mechanistic understanding the infection cycle could lead to insights for more effective anti-viral strategies. However, simulating the virus infection cycle is a computationally hard problem because some species fluctuate rapidly while others change gradually in number. We study vesicular stomatitis virus, an experimentally accessible virus for which we have developed a deterministic kinetic model of growth (Lim et al, PLoS Comp Bio, 2006). Stochastic simulation of VSV genome encapsidation is a computationally-intensive chain reaction that produces rapid fluctuations in the nucleocapsid(N) protein that associates with the genome. Analytical results from order statistics enable us to avoid explicit tracking of all intermediate species, with significant reduction in computational burden. This approach can be generalized to a broad class of stochastic polymerization reactions where multi-step chain reactions are modeled as a single reaction with time-delay.
One of the defining characteristics of life is the ability to
keep time, which organisms often achieve by using internal
genetic "clocks" to govern fundamental cellular behavior.
While the gene networks that produce oscillatory expression
signals are typically quite elaborate, certain recurring
network motifs are often found at the core of these biological
clocks. In this lecture, I will describe our recent
experimental and theoretical work on small genetic networks
exhibiting oscillatory behavior. One common motif which may
lead to oscillations is delayed auto-repression. We
constructed synthetic oscillators based on this design
principle, and observed robust and tunable oscillations both in
bacteria and yeast. Since genetic systems typically involve
small number of reacting components, intrinsic fluctuations
play an important role in the dynamics. I will introduce
theoretical and computational tools which can be used for
modeling stochastic dynamics of small genetic networks.
When describing chemical systems in a nonequilibrium steady state, striking differences arise in the predictions between deterministic and stochastic models. Using a classic model first introduced by Schlogl, we study the steady state behaviors of a one dimensional, open system that exhibits bistability. The definitions of thermodynamic quantities, such as flux, chemical potential, and entropy production rate are compared across the two types of models. The stochastic model allows for jumping between the two stable steady states which are separated by an unstable steady state in the deterministic model. The transition rates and exponentially small eigenvalue associated with this jumping are investigated and calculated numerically. Understanding the behavior of the stochastic steady state directly from the master equation is important because even widely used approximations such as Fokker-Planck can be incorrect, as is shown for this example.
This talk will provide an overview of computationally intensive
methods for conducting Bayesian inference for the rate constants of
stochastic kinetic intracellular reaction network models using
single-cell time course data. Inference for the true Markov jump
process is extremely challenging in realistic scenarios, so the true
model will be replaced by a diffusion approximation, known in this
context as the Chemical Langevin Equation (CLE). Inference for the CLE
is also challenging, but the development of effective algorithms is
possible, and turns out to be extremely effective, even in scenarios
where one would expect the diffusion approximation to break down.
Based on joint work with G. Craciun and J. W. Helton.
Dynamical system models of complex biochemical reaction
usually high-dimensional, nonlinear, and contain many
parameters. In some cases the reaction network structure
positive equilibria must be unique for all values of the
the model. In other cases multiple equilibria exist if and
special relationships between these parameters are
describe methods based on homotopy invariance of degree
which allow us
to determine the number of equilibria for complex
networks and how this number depends on parameters in the
The Supra-Chiasmatic Nucleus (SCN) is the central oscillator that keeps time in all mammalian organisms. Central to understanding the circadian rhythms that these oscillators generate are the multitude of interacting genes and proteins, in addition to the complex interactions of coupled SCN cells.
In recent year, quantitative modelling has emerged as an additional tool complimenting experimental techniques in the study of circadian rhythms. In this poster, we discuss a detailed mathematical model for circadian timekeeping within the SCN. Our proposed model consists of a large population of SCN neurons, with each neuron containing a network of biochemical reactions involving the core circadian components. Central to our work is the determination of the model’s unkown parameters, which were obtained from comparing the model’s output to experimental data. From these estimated parameters, additional experimental test of the model are proposed. Our studies highlight the importance of low numbers of molecules of clock proteins, and how this fact addects the accuracy of circadian timekeeping.
When a virus infects a living cell it directs the biosynthetic resources of the cell toward transcription and translation of viral proteins, replication of viral genomes, assembly of virus particles, and release of hundreds to thousands of progeny virus particles to the extracellular environment. For well-characterized viruses one may begin to build kinetic models that predict virus growth behavior based on the dynamics of the underlying molecular processes within the infected cell. Recent experimental study of single BHK cells infected by single particles of vesicular stomatitis virus, a virus that encodes only five genes, reveal a broad distribution of virus productivity. Viral genetic variation and host-cell heterogeneity are unable to fully account for the broad growth behavior. A better understanding of these experiments may follow from stochastic modeling of virus growth.