DTIC ADA458856: Stochastic Analysis of Gene Regulatory Networks using Finite State Projections and Singular Perturbation

0.41 MB·English

by  Defense Technical Information Center

Stochastic analysis of gene regulatory networks using finite state projections
and singular perturbation
Brian Munsky, Slaven Peles and Mustafa Khammash
Department of Mechanical Engineering
University of California
Santa Barbara, CA 93106-5070
Abstract—Considerable recent experimental evidence sug- This consideration has given rise to a relatively new branch
geststhatsignificantstochastic fluctuationsarepresentingene of computational research in chemistry: stochastic modeling
regulatory networks. The investigation of stochastic properties
of chemical reactions at the mesoscopic level.
in genetic systems involves the formulation of a mathematical
ForadiscretepopulationmodelswithN chemicalspecies,
representationofmolecularnoiseanddevisingefficientcompu-
tational algorithms for computing the relevant statistics of the the configuration of the system at any time is specified by
modeledprocesses.However,thecomplexityofgeneregulatory the population of each species: x:=[ ξ1 ξ2 ... ξN ]T.
networks poses serious computational difficulties and makes The configuration space of the system is the N dimensional
any quantitative prediction a difficult task. Monte Carlo based
non-negative integer lattice that contains all possible config-
approaches are typically used in study of complex stochastic
urations. Over time, the reacting system follows a trajectory,
systems, but they often suffer from long computation times,
slow convergence, and offer little analytic insight. The recently whereeachindividualreactionisrepresentedbyajumpfrom
proposed Finite State Projection (FSP) approach provides an the current configuration to another within the configuration
analytical alternative that avoids many of the shortcomings space. With the assumptions that the system is well mixed,
of Monte Carlo methods, but thus far it has only been
has constant volume and is kept at constant temperature,
demonstrated for a certain class of problems. In this paper
the problem becomes autonomous of time, and the waiting
weshowthattheapplicabilityofthefiniteprojectionapproach
canbeenhancedbytakingadvantageoftoolsfromthefieldsof time between jumps can be modeled as an exponentially
moderncontroltheoryanddynamicalsystems.Inparticular,we distributedrandomvariable[1].Theprobabilityofeachcon-
present an approach that utilizes singular perturbation theory figurationcanthenbedescribedwiththelineartimeinvariant
in conjunction with the Finite State Projection approach to
ordinarydifferentialequationknownastheChemicalMaster
improve the computation time and facilitate model reduction.
Equation [2].
Wedemonstratetheeffectivenessoftheresultingslowmanifold
FSPalgorithmonasimpleexamplearisinginthecellularheat Until recently, it was believed that the CME was in-
shock response mechanism. tractableexceptinthesimplestofcircumstances.Asaresult,
most discrete state stochastic models have been studied
I. INTRODUCTION almost exclusively with the simple Monte Carlo simulation
Through evolution living organisms have developed com- techniqueknownasthestochasticsimulationalgorithmSSA
plexrobustcontrolmechanismswithwhichtheycanregulate [3].Hererandomnumbersaregeneratedforeveryindividual
intracellular processes and adapt to changing environments. reactioneventinordertodetermine(i)whenthenextreaction
These mechanisms are triggered by small changes in sen- will occur, and (ii) which reaction it will be. However, for
sitive control parameters, yet their functions are unaffected most systems, huge numbers of individual reactions may
by relatively large variations in the system that may even occur, and the SSA is often too computationally expensive.
include damage of whole parts of the mechanism itself. Two types of recent approximations have been made to
Understanding control processes that take place inside a cell improve efficiency of the SSA: time scale based system
is key to understanding fundamental biological functions of partitioningroutinesandτ leapingroutines.Inthefirsttype,
living organisms. However, despite dramatic new develop- efficiency is gained by separating the system dynamics into
ments in experimental study of intracellular processes, the slow and fast parts [4]–[8]. If one is interested in short time
quantitative study of these systems is significantly hindered scales,the slowprocesses areconsidered tobe constant, and
by computational complexity. onlythedynamicsofthefastpartitionareconsidered.Ifone
Formanychemicalprocessesreactantsarepresentinlarge is interested in long term behaviors of the system, then the
numbersandthereactionscanbemodeledatthemacroscopic fast dynamics are averaged in order to concentrate upon the
level ordinary or stochastic differential equations (ODEs slow dynamics. In the second type of approximation to the
or SDEs). Inside a cell, however, matters can be entirely SSA,someprobabilisticaspectsofthesystemareassumedto
different. Here some chemicals such as proteins or RNA be constant over short discrete time intervals, and under this
molecules may be present in only a few copies per cell. assumption one can develop a faster Monte Carlo algorithm
Since these trace chemicals may effect important changes that leaps from one time step to the next [9]–[14]. Both
in the cell’s behavior, they cannot be ignored, yet they approximationtypes,systempartitioningandτ leaping,have
obviously cannot be treated as continuous concentrations. been very successful in increasing the scope of problems to
Report Documentation Page Form Approved
OMB No. 0704-0188
Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and
maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information,
including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington
VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it

does not display a currently valid OMB control number.
1. REPORT DATE 3. DATES COVERED
2006 2. REPORT TYPE 00-00-2006 to 00-00-2006
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER
Stochastic analysis of gene regulatory networks using finite state
5b. GRANT NUMBER
projections and singular perturbation
5c. PROGRAM ELEMENT NUMBER
6. AUTHOR(S) 5d. PROJECT NUMBER
5e. TASK NUMBER
5f. WORK UNIT NUMBER
7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION
Department of Electrical and Computer Engineering,University of REPORT NUMBER
California,Santa Barbara,CA,93106
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S)
11. SPONSOR/MONITOR’S REPORT
NUMBER(S)
12. DISTRIBUTION/AVAILABILITY STATEMENT
Approved for public release; distribution unlimited
13. SUPPLEMENTARY NOTES
14. ABSTRACT
15. SUBJECT TERMS
16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF
ABSTRACT OF PAGES RESPONSIBLE PERSON

a. REPORT b. ABSTRACT c. THIS PAGE 7
unclassified unclassified unclassified
Standard Form 298 (Rev. 8-98)
Prescribed by ANSI Std Z39-18
which the SSA may be reasonably applied. willtransitiontox duringthetimestep,dt.Thisprobability
i
Moving in an entirely different direction, we recently can be written as:
introducedtheFiniteStateProjection(FSP)methodasanew M !
X
approachfordirectlyapproximatingthesolutiontotheCME p (t+dt)=p (t) 1− a (x )dt
i i µ i
[15].TheFSPprovidesasystematicmethodtoperformbulk µ=1
system reductions on the possibly infinite dimensional ODE M
X
inordertofindafinitedimensionalanalyticalapproximation + p(x −ν ,t)a (x −ν )dt. (1)
i µ µ i µ
to the CME. Unlike Monte Carlo algorithms, the FSP does µ=1
not rely on random number generation and provides precise
From Eqn 1 it is easy to derive the differential equation
lower and upper bounds on its approximation accuracy.
known as the Chemical Master Equation, or CME [3]:
While the earliest implementations of the FSP have shown
M M
much promise in dealing with a handful of biologically X X
p˙ (t)=−p (t) a (x )+ p(x −ν ,t)a (x −ν ).
inspired models [15], [16], the method is still limited in i i µ i i µ µ i µ
µ=1 µ=1
problems where accuracy requirements demand considera-
(2)
tion of an exorbitantly large state space. In these cases,
model reduction techniques can be utilized to make the By combining all possible reactions that begin or end with
problemdimensionmoremanageable.Forexample,drawing the configuration point, xi, the time derivative of the proba-
uponconceptsofreachabilityandcontrolabilityfrommodern bility density of xi can be written in vector form as:
control theory, one can significantly reduce the order of the  −PM a (x ) T  p(x ,t) 
FSP [17]. Additionally, since the SSA benefits so strongly µ=1 µ i i
fromsystempartitioningmethods,itisnaturaltoaskwhether  a1(xi−ν1)   p(xi−ν1),t) 
tphaepesra,mweeorshsoimwilhaorwtoothlsemFaSyPimcapnrorveeapuphoungethbeeFnSePfi.tsInfrtohmis p˙i(t)= a2(xi...−ν2)   p(xi−...ν2),t) . (3)
time-scale separation based system partitioning approaches aM(xi−νM) p(xi−νM),t)
and the application of linear perturbation theory [18]. This
Based upon our enumeration of X, we can combine the
times-scale separation is exploited to develop the Slow-
equations for every p (t) into ordinary differential equation
Manifold FSP algorithm, which we present here. i
(ODE):
In the next section we provide some background on the P˙(t)=AP(t), (4)
basic FSP method. In section III we illustrate how linear
perturbation theory can be used to reduce the order of the where P(t) is the probability distribution at time t. The
FSP problem and we formulate the Slow-Manifold FSP generator matrix A is uniquely defined by the reaction stoi-
algorithm. In section IV we demonstrate this algorithm on chiometries and propensities and the choice of the enumera-
an example from the field of systems biology. Finally, in tionofX.Ahasthepropertiesthatitisindependentoft;all
section V, we conclude with remarks on the advantages of of its diagonal elements are non-positive; all its off-diagonal
these approaches over the original FSP and outline a few elements are non-negative; and all its columns sum to zero.
directions for future work on the topic. ThesolutiontoEqn4beginningatt=0andendingatt=t
f
is the expression: P(t )=Φ(0,t )·P(0). In the case where
f f
II. BACKGROUND the configuration set is finite dimensional, the operator,
We consider a well-mixed, fixed-temperature, and fixed Φ(0,t ), is the exponential of At , and one can compute
f f
volume system of N distinct reacting chemical species. the solution: P(t ) = exp(At )P(0). However, when A is
f f
At any point in time, we use an integer vector x := infinite dimensional or extremely large, the corresponding
[ ξ ξ ... ξ ]T to describe the discrete populations analytic solution is unclear or vastly difficult to compute.
1 2 N
oftheN molecularspecies.Eachxisauniqueconfiguration In these cases, one may devise a systematic means of
of the system, and the nonnegative set NN is the set of approximating the full system using finite dimensional sub-
all possible configurations. We will a priori fix a sequence systems. This truncation approach lies behind the rationale
x ,x ,...ofelementsinNN anddefineX:=[x ,x ,... ]T for the Finite State Projection (FSP) method [15].
1 2 1 2
astheorderedconfigurationset.Beginningatanyx ,suppose We must introduce some convenient notation. Let J =
i
that there can be a maximum of M possible reactions, {j ,j ,j ,...}denoteanindexset.ForanytwoindexsetsI
1 2 3
where each reaction leads to a different configuration: x → and J, let J ⊆I denote that I contains every element from
i
x = x +ν . For each µth reaction, ν is the associated J, and let JSI and JTI be the union and intersection,
j i µ µ
stoichiometry(directionaltransitionontheconfigurationset), respectively, of J and I. Let J0 denote the complement of
and a (x ) is the associated propensity function. Define the set J, where unless stated otherwise the complement is
µ i
p (t)=p(x ,t) as the probability that the system will have taken with respect to the set of all non-negative integers. If
i i
theithconfigurationattimet.Atanincrementallysmalltime X is an enumerated set {x ,x ,x ,...}, then X denotes
1 2 3 J
later,p (t+dt)isequaltothesumof(i)theprobabilitythat the subset {x ,x ,x ,...}. Furthermore, let v denote
i j1 j2 j3 J
the system begins in x at t and remains there until t+dt, the subvector of v whose elements are chosen according to
i
and (ii) the probability that the system is in x 6=x at t and J, andlet A denote thesubmatrix of A such thatthe rows
j i IJ
have been chosenaccording to I and the columnshave been Step 0 Choose an initial finite set of states, X , for the FSP.
Jo
chosen according to J. For example, if I and J are defined Initialize a counter, i=0.
as {3,1,2} and {1,3}, respectively, then: Step 1 Use propensity functions and stoichiometry to form A .
Ji
 a b c   g k  Compute ΓJi =1T exp(AJitf)PJi(0).
Step 2 If Γ ≥1−ε, Stop.
 d e f  = a c . Ji
exp(A t )P (0) approximates P (t ) to within ε.
g h k d f Ji f Ji Ji f
IJ Step 3 Add more states to find X .
Ji+1
For convenience let AJ denote the principle submatrix of A, Increment i and return to Step 1.
inwhichbothrowsandcolumnshavebeenchosenaccording
toJ.Unlessotherwisestated,allvectorinequalitiesaretobe
For the basic FSP algorithm, if we wish to find a solution
interpretedcomponentwise.Withthisnotationwecanrestate
thatisaccuratetowithinεatatimet ,wemustfindafinite
the following theorem from [15]: f
Theorem 1 If A ∈ Rn×n has no negative off-diagonal setofconfigurationssuchthattheprobabilityofeverleaving
that set during the time interval [0,t ] is less than ε. For
elements, then for any index set, J, f
many problems, including the examples shown in [15], [16],
[expA] ≥exp[A ]≥0. (5) this set of configurations may be small enough that we can
J J
easily compute a single matrix exponential to approximate
the solution to the CME. However, in other situations that
Ref. [15] provides a detailed proof of this theorem. The
may not be the case. In the next section we apply some
matrix A in Eqn 4 has no negative off-diagonal terms and
concepts of linear perturbation theory to reduce the order
therefore satisfies the assumptions of Thm 1. Consider two
of the FSP problem and thereby significantly improve the
finite configuration subsets X and X , where J ⊇ J .
J1 J2 2 1 efficiencyandexpandtheapplicabilityoftheFSPalgorithm.
Sincethefullprobabilitydensityvector,P(0)isnonnegative,
Thm 1 assures that: III. THESLOW-MANIFOLDFSPALGORITHM
[exp(A t )] P (0)≥exp(A t )P (0)≥0. In Step 1 of the FSP algorithm above we are faced with
J2 f J1 J1 J1 f J1 solving an N dimensional LTI system:
where P is the probability distribution of the elements
of X indJe1xed by J . This result guarantees that as one P˙(t)=AP(t), (9)
1
increases the finite configuration subset, the approximate
where A is made up from the propensity functions of the
solution increases monotonically for each element in the
various reactions. In gene networks we can often identify
configuration subset.
m clusters of configurations within which transitions occur
In addition to being non-negative, P(t) sums to exactly
quite frequently, while transitions between the clusters are
one. These properties and the nonnegativity of the off-
relatively rare. Let H be a generator matrix that is made
diagonalelementsofAallowonetostatethefollowingresult
up from the frequent transitions within these clusters, and
[15].
let G be the generator made up of the remaining transitions
Theorem 2 Consider a Markov process in which the
from one cluster to another such that A = H+G. With
probabilitydistributionevolvesaccordingtothelinearODE,
some permutation of the configuration space, H has a block
P˙(t) = AP(t), where A has no negative off-diagonal
diagonal structure, H = diag{H ,H ,...,H }, where
1 2 m
entries. If for some finite index set J, ε>0, and t ≥0,
f each block H is itself a generator corresponding to a single
i
1Texp(A t )P (0)≥1−ε, (6) cluster of fast interconnected configurations.
J f J
Each H is itself a generator and therefore has an eigen-
i
then
value equal to zero, which corresponds to a left eigenvector
exp(A t )P (0)≤P (t ), and (7) ui = 1T and a positive right eigenvector, vi, where vi
J f J J f
is scaled such that u v = 1Tv = 1. We use these
||P (t )−exp(A t )P (0)|| ≤ε. (8) i i i
J f J f J 1 eigenvectors to define.
While Theorem 1 guarantees that as we add points to    
u 0 ... v 0 ...
1 1
the finite configuration subset, the approximate solution
U= 0 u2 ... , and V= 0 v2 ... .
mofonhootwonicclaolslye itnhcereaapspesro,xTimheaotiroenmi2s ptorovthideestruaecesortliufitcioante.  ... ... ...   ... ... ... 
Together the two theorems lead us to the FSP algorithm
Let S = [ V S ] be a square matrix in which each
presented in [15]: R
column is an right eigenvector of H that has been scaled
to have an absolute sum of one. The inverse of S is given
The Finite State Projection Algorithm (cid:20) (cid:21)
U
byS−1 = suchthatwehavethefollowingsimilarity
Inputs Propensity functions and stoichiometry for all reactions. SL
Initial probability density vector, P(0). transformation for H:
(cid:20) (cid:21)
Final time of interest, t . 0 0
f S−1HS= , Λ=diag(λ ,...,λ ).
Total amount of acceptable error, ε>0. 0 Λ m+1 N
where the first m diagonal elements correspond to the zero build matrices U and V.
eigenvalues of the H blocks, and we label non-zero eigen- Estimate ε=||G V|| /|λ |.
i Jk 1 m+1
values of H so that 0 > Re{λ } ≥ Re{λ },... ≥ Compute transient error γ =|S P(0)| exp(λ t ).
m+1 m+2 L 1 m+1 f
Re{λ }. If we apply the coordinate transformation Step 3 Find PFSP(t )=Vexp(UG Vt )UP (0)
N Jk f Jk f Jk
(cid:20) y1(t) (cid:21)=S−1P(t)=(cid:20) UP(t) (cid:21), Step 4 If Γand≥co1m−puδt,eSΓtoJkp.=1TPFJkSP(tf).
y2(t) SLP(t) PFSJPk(t ) approximates P (t ) within δ+γ+O(ε).
Jk f Jk f
Eqn 9 becomes: Step 5 Add more states to find X .
Jk+1
(cid:20) (cid:21) (cid:20) (cid:21)(cid:20) (cid:21) Increment k and return to Step 1.
y˙ (t) UGV UGS y (t)
1 = R 1 .
y˙ (t) S GV Λ+S GS y (t)
2 L L R 2
(10)
Above we have used the non-traditional error estimate
If we can make a clear distinction between frequent and
notation δ+γ+O(ε) to mean the following. If δ is largest,
rare reactions, then we are assured that the magnitude of
then the dominant error is most likely the result of the
GV is small with respect to Λ. For convenience we define
projection, and the slow manifold truncation error can be
ε = 1 ||GV|| and Q = 1GV, and we can make
|Re{λm+1}| 1 ε ignored. If γ is largest then the time tf is too short for
coordinatesubstitution(y1,y2)=(y,εz)togetthestandard thetransientdynamicstosufficientlydiminishandadditional
singular perturbation model, eigenvectors must be included in the truncation. Finally, if
(cid:20) y˙(t) (cid:21) (cid:20) εUQ εUGS (cid:21)(cid:20) y(t) (cid:21) ε is larger than δ and γ, then there is insufficient separation
= R .
εz˙(t) εS Q ε(Λ+S GS ) z(t) between the slow and fast dynamics and an alternative
L L R
(11) reduction scheme may be required.
Wecanapplylinearperturbationtheorytoshowthat|z(t)|≤ In the following section we will illustrate this algorithm
γ(t)+O(ε), where γ = |z(0)|exp(λ t) is the transient on model of the heat shock response in E. coli.
m+1
error of the fast manifold. For the slow manifold, we have
IV. EXAMPLE
y(t)=exp(UGVt)y(0)+O(ε),
In order to deal with frequently changing situations,
for all time t≥01. living organisms require mechanisms to deal with harsh
The reverse transformation then yields environments.Onesuchmechanismisthecellularheatshock
response. At elevated temperatures proteins inside the cell
(cid:20) (cid:21) (cid:20) (cid:21)
y (t) y(t)
P=S 1 =S begintomisfoldandlosetheirgeometrydependentfunction-
y (t) εz(t)
2 ality. In an attempt to stop this misfolding, the cell produces
=Vy(t)+O(ε) heat shock proteins including molecular chaperones and
=Vexp(UGVt)y(0)+O(ε) proteases.Byhelpingtorefolddeformedproteins,molecular
chaperones can restore the original functionality of those
=Vexp(UGVt)UP(0)+O(ε). (12)
proteins.Proteases,ontheotherhand,helpdegradedamaged
It is important to note that due to truncation of Eqn 10, proteins before those damaged proteins impair the function
only contributions of the first m eigenvectors of H affect of the cell. The core of the heat shock response mechanism
the approximate solution. Therefore, instead of calculating in E. coli is the synthesis of the σ -RNAP complex [19].
32
full eigensystems for each block H , it suffices to find only Here we study a simplified model for σ -RNAP synthesis
i 32
eigenvectors associated with zero eigenvalues. using the slow manifold finite state projection algorithm.
Applying this model reduction to the original FSP Atnormalphysiologicaltemperaturesσ proteinisfound
32
algorithm yields the following algorithm which we name almost exclusively in a complex σ -DnaK. However, as the
32
the Slow-Manifold FSP algorithm: temperature increases this complex becomes less stable and
the probability of finding free σ inside the cell increases.
32
The Slow-Manifold FSP Algorithm Once free, σ can combine with RNA polymerase to form
32
a stable σ -RNAP complex, which initiates production of
Inputs Propensities and stoichiometries for all reactions. 32
heat shock proteins. This simple regulatory mechanism is
Initial probability density vector, P(0).
summarized by a set of three reactions,
Final time of interest, t .
f
Target FSP error, δ >0. s (cid:10)s →s , (13)
1 2 3
Step 0 Choose initial set of states, X , for the FSP.
Jo
Initialize a counter, k =0. where s1, s2 and s3 correspond to the σ32-DnaK complex,
Step 1 Use fast reactions connecting states within XJk to the σ32 heat shock regulator and the σ32-RNAP complex,
form H =diag{H ,...,H }. respectively. This model of the heat shock subsystem has
Jk 1 m
Use remaining reactions to form G . been analyzed before using various computational methods
Jk
Step 2 Find eigenvalues and vectors of each H and including Monte Carlo implementations [7], [20].
i
In the biological system, the relative rates of the reactions
1Foradetailedderivationoftheseerrors,seetheappendixof[18]. are such that the reaction from s2 to s1 is by far the fastest,
Step0:Wespecifytheproblemparametersasgivenabove.
7 We choose a target 1-norm error of δ = 0.001 for the
distribution at time t , and we choose an initial FSP set
6 f
X that includes all configurations such that s ≤200 and
5 J0 3
s ≤11.
s 4 2
3 Step 1: We form the fast generator H . For the config-
3 Jk
urations in X , reaction 1 has propensities ranging from
2 J0
18,000 to 20,000, reaction 2 has propensities ranging up
1 to 110,000, and reaction 3 has a propensity that ranges no
0 higher than 22. Thus there is a clear separation between the
0 1 2 3 4 5 s6 7 8 9 10 11 12 fast transitions (reactions 1 and 2) and the slow transitions
2
(reaction 3). Furthermore, since s changes only during the
3
(a)
slow reaction, each level of s = i defines a cluster of fast
3
interconnected configurations; these are shown as horizontal
rows in Figure 1. We will use the index set I to denote the
i
set of configurations in the ith cluster. The generator H is
H i
3 formed from all of the reactions that begin and end within
H2 X .Inthiscasethereare12pointsineachith clustercorre-
Ii
H spondingtos =i,ands ={0,1,2,...,11};ifwewereto
1 3 2
H change the number of points in our projection, the size and
0
number of clusters would change accordingly. The full fast
0 1 2 3 4 5 6 7 8 9 10 11
generator is defined H = diag(H ,H ,...,H ). For
(b) J I0 I1 I200
the reaction rates given above, the first nonzero eigenvalue
Fig.1. (a)Twodimensionalintegerlatticerepresentingpossibleconfigu- of H can be computed to be λ ≈−4×104.
J m+1
rationsofthetoyheatshockmodel.Heres2ands3arepopulationsoffree The rare transition generator G is then formed from
σ32 molecules and σ32-RNAP compounds, respectively, while s1 is the J
populationofσ32-DnaKcompounds.Reactionss1 (cid:10)s2,arerepresented the remaining transitions coming from two components: (i)
bwyithbiddiiargecotnioalnaalrrhoworsi.zoTnhtealtoatrarlownusmabnedrorefaσct3i2onisscso2ns→tants,3soisthreepcrheesmenitceadl rXeactio,nan3dc(oiri)reaspllontrdainnsgititoonsatatrkainnsgititohne fsryosmtemsetfroXmIithtoe
state of the system is uniquely defined by s2 and s3 alone. (b) The same Ii+1
lattice after applying the finite state projection. Unlikely states have been set XJ to its complement XJ0.
aggregatedintoasinglesinkstate.Eachhorizontalrowofconfigurationsis Step 2: We build the matrices U and V from the left and
separatedfromtherestbytheslowreaction3andtheisusedtoformthe
right eigenvectors corresponding to the zero eigenvalues of
fastblockgeneratorHi
H . We compute the norm ||G V|| ≈2.0, and use this to
J J 1
estimate how well the time scale is separated between the
frequent and rare events:
and σ molecules infrequently escape from DnaK long
32
enough to form the σ32-RNAP complex. The purpose of ε= ||GJV||1 = 2 =0.5×10−4.
this mechanism is to strike a balance between fixing the |λ | 4×104
m+1
damage produced by heat and saving the cell’s resources,
We also compute the transient error due to the initial con-
as a significant portion of cell energy is consumed when
ditions, γ = |S P(0)| exp(λ t ) ≈ exp(−1.2×107),
producing heat shock proteins. The optimal response to the L 1 m+1 f
which is far smaller that δ, and can therefore be neglected.
heat shock is not massive, but measured production of heat
Step3:Weapplyperturbationtheoryandapproximatethe
shock proteins, which leaves sufficient resources for other
FSP solution as
cellular functions. We use the following set of parameters
values for the reaction rates: PFSP(t )≈Vexp(UG Vt )UP (0),
Jk f Jk f Jk
c =10, c =4×104, c =2, which has error γ +O(ε) compared to the unreduced FSP
1 2 3
s (0)=2000, s (0)=s (0)=0. solution. In order to determine whether this FSP solution is
1 2 3
sufficient, we compute Γ = 1TPFSP(t ) = 7.3×10−8.
Jk Jk f
With only the reactions above, the total number of σ – It is not.
32
free or in compounds–is constant, so that s +s +s = Step 4: Since the Γ < 1 − δ, we need to add more
1 2 3 J0
const. With this constraint the reachable states of this three configurations and return to Step 1. If we define X to
J1
species problem can be represented on a two dimensional include all configurations such that s ≤ 250 and s ≤
3 2
latticeasshowninFigure1a.Forourinitialconditionsthere 11, we find Γ = 0.034, which is little better. If we
J1
are 2,001,000 reachable states, and the full chemical master further expand the configuration set to include up to 300
equation is too large to be tackled directly, and we wish to or 350 molecules of s , we compute Γ = 0.921 and
3 J2
find an approximate solution using the slow manifold FSP Γ = 0.99997 respectively. The FSP solution on the set
J3
algorithm at time t =300s. To do this we follow the steps X exceeds the precision of our stopping criteria, and we
f J3
of the algorithm as follows: know that PFSP(t ) approximates the true solution within
J3 f
TABLEI
0.03
(TOP)ACOMPARISONOFTHECOMPUTATIONALEFFORTAND
0.025 ACCURACYFORTHEFSPANDTHESLOW-MANIFOLDFSP
y
sit ALGORITHMSFOREACHEXPANDINGCONFIGURATIONSET
en 0.02 XJ1 ⊂XJ2 ⊂XJ3 .(BOTTOM)ACOMPARISONOFTHETOTAL
D COMPUTATIONALEFFORTANDACCURACYOFTHEFSP,THE
bility 0.015 SLOW-MANIFOLDFSP,THESSAANDASLOWSCALESSAALGORITHM
a FORTHESOLUTIONOFTHECMEARISINGINTHETOYHEATSHOCK
b 0.01
o MODELATtf =300s.ALLCOMPUTATIONSHAVEBEENPERFORMEDIN
r
P MATLAB7.2ONA2.0MHZPOWERPCDUALG5.
0.005
Configurationset ||Error|| Comp.Time(s)
1
2000 220 240 260Pop280ula3t0i0on 32o0f s340 360 380 400 FSP Slow-ManifoldFSP
3 XJ1:s3≤250,s2≤11 0.97 253 0.39
XJ2:s3≤300,s2≤11 0.08 439 0.54
Fig.2. Probabilitydistributionfors3 calculatedattf =300s.Theslow XJ3:s3≤350,s2≤11 2×10−5 647 0.67
manifoldFSPsolution(dots)approximateswelltheFSPsolution(smooth
Method #Simulations Time(s) ||Error||
solid lines). The jagged grey line represents the corresponding slow scale 1
SSAsolutionafter104 realizations.SeealsoTableI. FSP –a 1472 ≤2×10−5
SSA 103 >20000 ≈0.25
Slow-ManifoldFSP – 1.88 ≈6.6×10−4
δ = 0.001. Furthermore, since ε = 0.0005 and γ (cid:28) δ, the SlowScaleSSA 103 82 ≈0.24
SlowScaleSSA 104 826 ≈0.066
error introduced by the time scale separation is of the same
SlowScaleSSA 105 8130 ≈0.027
order.
Forthisproblemweareinterestedindetermininghowthe aTheFSPandSlow-ManifoldFSPalgorithmsarerunonlyonce
population of σ -RNAP compounds grows in time if the withaspecifiedallowabletotalerrorof10−3.
32
temperatureisconstantandslightlyabovenormalphysiolog-
icallevel.Thisnumberisproportionaltothenumberofheat
shock proteins produced in the cell. With the original FSP bulk system reductions could bring models closer into the
algorithm,wehavecomputedtheprobabilitydistributionfor fold of solvable problems. Here we have shown that the
s3 at tf = 300s, and Figure 2 illustrates this probability FiniteStateProjectionmethodcanbefurtherenhancedwhen
distribution. The FSP solution is shown with solid lines, and solving the chemical master equation for systems involving
the dots represent the distribution of s3 as computed using multiple time scales. In this work we have presented a
our current slow manifold FSP algorithm. The difference Slow Manifold version of the FSP algorithm. Based upon
between the two results is indistinguishable. However, if we singular perturbation theory and our previous work in [18],
are to compare the relative computational effort of the two this algorithm provides a powerful computational tool for
algorithms, we find a significant improvement. Table I(top) studyingintracelularprocessesandgeneregulatorynetworks.
shows the computational time and accuracy for each step By casting the stochastic chemical kinetics problem in a
of the FSP and the slow manifold FSP algorithm. From familiar linear system context, the FSP provides valuable
the table we find that the time scale reduction reduces the insight into the bewildering complexity that intracellular
computational effort at each step by a factor of almost processes exhibit. Model reductions not only speed up
1000 compared to the original FSP algorithm. For further computations, but they also enable us to distinguish which
comparison I(bottom) shows the total computational effort aspects of a system are important and which are not. It is
oftheFSPandtheslowmanifoldFSP,aswellastwoMonte plausible that the main features of intracellular dynamics
Carlo methods: the original SSA and an SSA approximation can be captured in a relatively small subset of the state
very similar to the Slow Scale SSA [7] in which we have space, as the results obtained by FSP reduction on the heat
appliedthesametimescalereductionapproach.The1-norm shock model suggest. Depending on the observation time of
difference between the original FSP and the slow manifold interest, some reactions can be neglected, while some will
FSP was found to be 6.6×10−4, which is indeed on the contribute only through their averages. Preliminary success
same order as ε=5×10−4 and is smaller than the required with our approach gives us a hope that relatively simple
errortoleranceδ =1×10−3.However,bymakingthissmall models for intracellular processes can be discovered in the
sacrifice in accuracy, the slow manifold FSP converges in a process of reducing the FSP solution.
fraction of the time of any of the other methods.
Of course, one can easily envision that additional model
reductionsmaybepossibletoevenfurtherenhancethepower
V. CONCLUSION
of the Finite State Projection. Indeed other reductions based
Until recently, it was thought that the chemical master upon control theory [17] are already becoming apparent.
equation could not be solved analytically except for the Also, in our computations we have used off the shelf nu-
most trivial systems. Previous work on the Finite State merical routines for eigensystem calculations and matrix ex-
Projection demonstrated that for many biological systems, ponentiation. Further improvements in computational speed
can be achieved if these routines are optimized for matrices
which define master equations and their special properties.
We intend to investigate these possibilities in the future.
ACKNOWLEDGMENT
This material is based upon work supported by the
National Science Foundation under Grant NSF-ITR CCF-
0326576 and the Institute for Collaborative Biotechnologies
through Grant DAAD19-03-D-0004 from the U.S. Army
Research Office.
REFERENCES
[1] vanKampen,StochasticProcessesinPhysicsandChemistry,3rded.
Epsevier,2001.
[2] D.McQuarrie,“Stochasticapproachtochemicalkenetics,”J.Applied
Probability,vol.4,pp.413–478,1967.
[3] D. T. Gillespie, “Exact stochastic simulation of coupled chemical
reactions,”J.Phys.Chem.,vol.81,no.25,pp.2340–2360,May1977.
[4] E.HaseltineandJ.Rawlings,“Approximatesimulationofcoupledfast
andslowreactionsforstochasticchemicalkinetics,”J.Chem.Phys.,
vol.117,no.15,pp.6959–6969,Jul.2002.
[5] C. V. Rao and A. P. Arkin, “Stochastic chemical kinetics and the
quasi-steady-stateassumption:Applicationtothegillespiealgorithm,”
J.Chem.Phys.,vol.118,no.11,pp.4999–5010,Mar.2003.
[6] Y. Cao, D. T. Gillespie, and L. R. Petzold, “Multiscale stochastic
simulation algorithm with stochastic partial equilibrium assumption
forchemicallyreactingsystems,”J.Comp.Phys.,vol.206,no.2,pp.
395–411,July2005.
[7] Y. Cao, D. Gillespie, and L. Petzold, “The slow-scale stochastic
simulationalgorithm,”J.Chem.Phys.,vol.122,no.014116,Jan.2005.
[8] H.SalisandY.Kaznessis,“Accuratehybridstochasticsimulationofa
systemofcoupledchemicalorbiologicalreactions,”J.Chem.Phys.,
vol.112,no.054103,2005.
[9] D. T. Gillespie, “Approximate accelerated stochastic simulation of
chemically reacting systems,” J. Chem. Phys., vol. 115, no. 4, pp.
1716–1733,Jul.2001.
[10] D. T. Gillespie and L. R. Petzold, “Improved leap-size selection for
accelerated stochastic simulation,” J. Chem. Phys., vol. 119, no. 16,
pp.8229–8234,Oct.2003.
[11] M. Rathinam, L. R. Petzold, Y. Cao, and D. T. Gillespie, “Stiffness
in stochastic chemically reacting systems: The implicit tau-leaping
method,” J. Chem. Phys., vol. 119, no. 24, pp. 12784–12794, Dec.
2003.
[12] T. Tian and K. Burrage, “Binomial leap methods for simulating
stochastic chemical kinetics,” J. Chem. Phys., vol. 121, no. 21, pp.
10356–10364,Dec.2004.
[13] Y. Cao, D. T. Gillespie, and L. R. Petzold, “Avoiding negative
populationsinexplicitpoissontau-leaping,”J.Chem.Phys.,vol.123,
no.054104,2005.
[14] Y.CaoandL.R.Petzold,“Accuracylimitationsandthemeasurement
oferrorsinthestochasticsimulationofchemicallyreactingsystems,”
J.Comp.Phys.,vol.212,no.1,pp.6–24,Feb.2006.
[15] B.MunskyandM.Khammash,“Thefinitestateprojectionalgorithm
forthesolutionofthechemicalmasterequation,”J.Chem.Phys.,vol.
124,no.044104,2006.
[16] B. Munsky, A. Hernday, D. Low, and M. Khammash, “Stochastic
modeling of the pap-pili epigenetic switch,” Proc. FOSBE, pp. 145–
148,August2005.
[17] B. Munsky and M. Khammash, “A reduced model solution for the
chemical master equation arising in stochastic analyses of biological
networks,”Proc.45thIEEECDC,Dec.2006.
[18] S. Peles, B. Munsky, and M. Khammash, “Model reduction using
multiple time scales in stochastic gene regulatory networks,” no.
CCDC-06-0828,2006.
[19] H. El Samad, H. Kurata, J. Doyle, C. Gross, and K. M., “Surviving
heatshock:Controlstrategiesforrobustnessandperformance,”PNAS,
vol.102,no.8,p.27362741,2005.
[20] H.ElSamad,M.Khammash,L.Petzold,andD.Gillespie,“Stochastic
modelingofgeneregulatorynetworks,”Int.J.RobustNonlin.,vol.15,
pp.691–711,2005.