Approximation Techniques for Stochastic Analysis of Biological Systems
Thakur Neupane, Zhen Zhang, Curtis Madsen, Hao Zheng, Chris J. Myers
AApproximation Techniques for StochasticAnalysis of Biological Systems
Thakur Neupane, Zhen Zhang ∗ , Curtis Madsen, Hao Zheng, and Chris J.Myers Abstract
There has been an increasing demand for formal methods in thedesign process of safety-critical synthetic genetic circuits. Probabilistic modelchecking techniques have demonstrated significant potential in analyzing theintrinsic probabilistic behaviors of complex genetic circuit designs. However,its inability to scale limits its applicability in practice. This chapter addressesthe scalability problem by presenting a state-space approximation method toremove unlikely states resulting in a reduced, finite state representation of theinfinite-state continuous-time Markov chain that is amenable to probabilisticmodel checking. The proposed method is evaluated on a design of a genetictoggle switch. Comparisons with another state-of-art tool demonstrates bothaccuracy and efficiency of the presented method.
Thakur NeupaneUtah State University, Logan, UT USA, e-mail: [email protected]
Zhen ZhangUtah State University, Logan, UT USA, e-mail: [email protected]
Curtis MadsenBoston University, Boston, MA USA, e-mail: [email protected]
Hao ZhengUniversity of South Florida, Tampa, FL USA, e-mail: [email protected]
Chris J. MyersUniversity of Utah, Salt Lake City, UT USA, e-mail: [email protected] ∗ Corresponding author. 1 a r X i v : . [ c s . ET ] J a n Neupane et al.
Computational biologists typically construct models to better understandand explore the possible behaviors of biological systems [35]. By using formalmethods, such as model checking, to analyze these models, researchers areable to ensure that certain properties hold in biological system designs [19]. Inorder to numerically model check a system, the system’s state space must beenumerated. For systems that are highly concurrent and have infinite states,such as genetic circuits (i.e., the collections of genes within DNA that interactto control the behavior of cells, see Section 4 for more details), enumeratingthe state space can be computationally intractable due to the state spaceexplosion problem. Techniques such as partial order reduction that reducethe number of reachable states in a system have shown some promise intackling this problem [3, 4, 9], but these methods often rely on transitiondependencies based on the disablings (and/or enablings) and commutativityof independent transitions. Most models of genetic circuits do not containtransitions that disable/enable other transitions leading researchers to seekother solutions to this problem.Another way to reduce the state space of a system is to introduce thresholdabstractions to collapse multiple states of the system together [31]. This typeof abstraction works very well in systems where there exist groups of statesin equivalence classes. This is often the case in genetic circuits where firing asingle transition does not have a great effect on the likelihood of firing othertransitions in the system. Although this type of abstraction has previouslybeen successfully applied to genetic circuits, selecting the threshold values iscurrently done in a manual ad hoc manner.This chapter presents an alternative method for deriving a reduced, fi-nite state representation of a genetic circuit’s behavior. This method worksby computing the approximate probability of reaching each state on-the-fly and stops exploring different paths when the cumulative path probabil-ity drops below a predetermined value, and these paths are routed to anabstract absorbing state, which accumulates probability leakage during theMarkovian analysis. The resulting continuous-time Markov chain (CTMC)can be analyzed using probabilistic model checking approaches to determinethe probability that the original genetic circuit satisfies a desired temporallogic property given in continuous stochastic logic (CSL) [2, 25]. This chapterillustrates the utility of this method by applying it to a model of the genetictoggle switch and by comparing the results to a previous approach wherethe thresholds were determined by hand to produce the finite state repre-sentation [29, 31]. Additionally, this method is compared with a state-of-artstochastic hybrid analysis tool on several benchmarks, and comparisons ofresults demonstrate both accuracy and efficiency of our proposed method. pproximation Techniques for Stochastic Analysis of Biological Systems 3
To improve the scalability of probabilistic model checking, bisimulation mini-mization (e.g., [11, 12, 13]) has been extended to the probabilistic setting [22]to achieve up to a logarithmic state space reduction.
Probabilistic abstrac-tion (e.g., [23, 10, 21]) applies coarser state merging to achieve better reduc-tion, while ensuring a simulation relation between the abstract and concreteMarkov models. A transition on the abstract Markov model has a range ofprobabilities, represented by an interval with the maximal and minimal prob-abilities for taking the transition. In particular, [23] presents a theoreticalframework for reducing discrete-time Markov chains (DTMCs) and CTMCsusing a three-valued abstraction and for model checking these abstractions.However, how to partition the state space in this framework is not discussed,nor is the refinement of the abstractions in the case that inconclusive resultsare produced. Although these reduction techniques can be powerful, theymay not be effective in alleviating the exponential state growth caused byconcurrency as they are not designed to tackle concurrency in the first place.Unfortunately, concurrency is inherent in most synthetic biological systems.To address the state explosion problem, some approaches attempt to trun-cate the state space. For instance, [32] presents a method for selectively ex-ploring states involving rare events; however, this technique requires the mod-ification of parameters in the system to help guide this exploration. Otherapproaches attempt to dynamically explore the state space and continuallyadd states until the resulting state space satisfies a desired level of preci-sion [5, 33, 34].A probabilistic counter-example guided abstraction refinement approach isdeveloped in [21, 38]. Predicate abstraction is applied to programs in a proba-bilistic guarded command language, and counter-examples are represented asfinite Markov chains, where additional predicates are extracted by using anSMT solver in the case that such counter-examples are spurious. [26] presentsa compositional verification approach to probabilistic systems using assume-guarantee reasoning. Both component assumptions and guarantees are rep-resented as probabilistic safety properties. Component verification can beexpensive in this approach as it is reduced to a linear programming problem.Furthermore, assumptions are derived manually. Additionally, [38, 26, 21]are all based on probabilistic automata, which support non-determinism butwith discrete-time semantics.In [31], genetic circuit models are converted into CTMCs using operatorsite reduction abstractions relying on quasi steady-state approximations. Toavoid the state explosion problem, the authors employ a state aggregationmethod to collapse states together based on user provided thresholds. Whilethe application of probabilistic model checking to the reduced CTMC canproduce results in a fraction of the time of simulation-based approaches, thismethod is incapable of quantifying the error introduced by this aggregationand relies on user input for good choices of thresholds. While there has been
Neupane et al. work to address the former [1], our method attempts to alleviate the latterby automatically determining a finite state representation by removing statesthat are found to be extremely unlikely during the generation of the CTMCfrom the genetic circuit model.A similar approach to the one presented in this chapter is the sliding win-dow abstraction [20]. This method approximates a solution to the chemicalmaster equation (CME) by dividing the time period of interest into smalltime steps, iteratively constructing a window of an abstract state space thatpreserves the probability mass at the current time step, and then “sliding”the window in each subsequent time step to include newly generated stateswith significant probability while abstracting away those with negligible prob-ability contribution until the last time step has elapsed. This method ef-fectively performs transient CTMC analysis on a manageable approximatedstate space and successively updates a state space approximation by followingthe direction in which probability mass moves as time evolves. The abstractstate space construction is based on a worst-case estimation of lower andupper bounds on the populations of the chemical species.A more recent improvement of the sliding window implementation is the
STochastic Analysis of biochemical Reaction networks (STAR) [28]. It com-putes approximate solutions to population Markov processes using a stochas-tic hybrid model that combines moment-based and state-based representa-tions of probability distributions, and has been optimized to drop unlikelystates and add likely states on-the-fly.Our approach differs from the sliding window method in that it does notrequire many manual factors (e.g., several different initial states to com-pute a state update, a limited window size, etc.) to compute its state space.Additionally, the method presented in this chapter has the potential to op-timize the choice of the termination indicator factor to preserve accuracywhile requiring a manageable state space. Finally, our approach is based ona reaction-based abstraction model, and as a result, is readily applied to ge-netic circuit models while the method in [20] focuses on Markov chains thatare specified by a finite set of transition classes.
The high-level modeling formalism used in this chapter is the stochastic chem-ical kinetic (SCK) model [35].
Definition 1.
A SCK model is a tuple M = (cid:104) S , R , x (cid:105) which is composed of n chemical species S = { S , . . . , S n } , m reaction channels R = { R , . . . , R m } ,and an initial molecule count of each chemical species at the beginning ofanalysis (i.e., x : S n → N ). A reaction R i = (cid:104) α i , v i (cid:105) includes a propensityfunction α i : N n → R + that corresponds to the probability of a reaction, and pproximation Techniques for Stochastic Analysis of Biological Systems 5 the state change vector v i ∈ Z n that corresponds to the change in moleculecount for each species due to reaction R i .A reaction R i can occur in state x ∈ X , if its propensity is greater thanzero (i.e., α i ( x ) > α i essentially determines thelikelihood that R i occurs in the current state. After a reaction R i occurs, thestate is updated as follows: x (cid:48) = x + v i .The execution of reactions in an SCK model creates a state graph as definedbelow: Definition 2.
A SG is a tuple G = (cid:104) X , δ, x (cid:105) where– X is a non-empty set of states,– δ ⊆ X × R × X is the set of state transitions.– x : S n → N is the initial state.Note that |G | represents the state count of G .For most SCK models of real biological networks, they incur an infinitenumber of states. Therefore, the goal of this chapter is to find a finite subset ofthe states that sufficiently represents states that are actually likely to occur.Once a finite state graph is obtained, properties can be verified on this stategraph using probabilistic model checking. Probabilistic model checking is a formal verification method for checkingquantitative properties of probabilistic systems. The models of interest in-clude DTMCs and CTMCs, both of which belong to a class of stochasticprocesses that are used to reason about random phenomena in applicationdomains such as synthetic biology. Both Markov models are essentially atransition system with each transition labeled by a discrete probability forDTMCs or a transition rate for CTMCs. A DTMC is a transition systemwith a discrete probability labeled on each transition [25], which describesthe likelihood of a single step moving from one state to another. A CTMC,on the other hand, is a transition system with a transition rate r ( s, s (cid:48) ) la-beled on the transition emanating from state s to s (cid:48) . This rate determines theprobability of executing this transition within t time units, which is e − r ( s,s (cid:48) ) t .The rate r ( s, s (cid:48) ) uniquely characterizes an exponential distribution to governthe average state residence time of state s , which is r ( s,s (cid:48) ) . CTMCs allowfor modeling of real-time systems, as the delay of a transition can be anyarbitrary real value.Properties to verify using probabilistic model checking are specified using Probabilistic Computation Tree Logic (PCTL) [18] for DTMCs and CSL forCTMCs. PCTL extends
Computation Tree Logic (CTL) [6] by replacing exis-tential and universal path quantifiers with a probability operator, and henceexpresses probabilistic properties for a DTMC. In addition to path probabil-ities, two traditional properties of CTMCs are the transient and steady-state behaviors. Transient analysis reports the probability of being in each state ofthe Markov chain at a particular time instant, and steady-state analysis gives
Neupane et al. the corresponding probability in the long-run. Model checking algorithms forPCTL (e.g. [7, 8, 18]) have identical structure to the model checking algo-rithm for CTL. Model checking CTMC first discretizes the CTMC into an embedded
DTMC, from which many properties of the corresponding CTMCcan be deduced, for example, checking state reachability properties regard-less of how long it takes, and the expected time objectives. For checking statereachability within some time bound, the CTMC is discretized into a uni-formized
DTMC with the iterative numerical method uniformization [16, 17].The uniformized DTMC preserves the state resident time so that its transientbehavior is equal (up to some accuracy) to the corresponding CTMC.In order to perform probabilistic model checking on CTMCs, CSL can beused. CSL properties consist of state formulae (formulae that are either trueor false in a specific state) and path formulae (formulae that are either true orfalse along a specific path). CSL properties are specified using the followinggrammar:
Prop ::= U (T , Ψ, Ψ ) | F (T , Ψ ) | G (T , Ψ ) | St ( Ψ ) Ψ ::= true | Ψ ∧ Ψ | ¬ Ψ | φ (cid:62) φ | φ > φ | φ = φφ ::= v i | c i | φ + φ | φ − φ | φ ∗ φ | φ/φ | Prop
T ::= true | T ∧ T | ¬ T | t (cid:62) c i | t > c i | t = c i where v i is a variable, c i is a constant, and t stands for time in the system. InCSL, Ψ is a state formula that can be either comparisons between numericalexpressions, φ , or other state formula combined using logical connectives.A CSL property, Prop , is a path property over state formula. For example,the
Until property is of the form U ( T, Ψ , Ψ ), and it returns the probabilitythat along paths originating in the current state, Ψ remains true until Ψ becomes true during the time specified by time expression, T . The eventuallyoperator, F , is simply a shorthand for an until property where Ψ is true . Theglobally true operator, G , is another shorthand that specifies that Ψ remainstrue during the time in which T evaluates to true. The steady-state operator, St , returns the probability that when the SCK model reaches a steady statethat it has reached a state where Ψ is true. Finally, CSL formulae, Prop , canbe nested within other formula, creating recursive properties.For example, the CSL property St ( x > ∧ y (cid:62)
10) would return theprobability that in the steady state, the system reaches a state where thevariable x is greater than 5 and the variable y is greater than or equal to 10.Alternatively, the CSL property F ( t > ∧¬ ( t (cid:62) , x > ∧ y (cid:62)
10) wouldreturn the probability that the system follows an execution path originating inthe initial state where the variable x becomes greater than 5 and the variable y becomes greater than or equal to 10 sometime between 100 and 200 timeunits non-inclusive. For a path to satisfy this property, the system does notneed these conditions to hold true for the entire 100 time unit interval; theyjust both need to become true simultaneously at some point within this timeframe. pproximation Techniques for Stochastic Analysis of Biological Systems 7 A genetic circuit is constructed using DNA, and it typically includes, ata minimum, regions that act as promoters , ribosome binding sites (RBS), coding sequences (CDS), and terminators . The promoters are regions where transcription is initiated when an RNA polymerase (RNAP) molecule binds,and then begins to walk along the DNA copying the sequence to form a mes-senger RNA (mRNA) molecule until it reaches the location of the terminator.The terminator causes the RNAP to be released and thus ends transcription.The RBS region when copied to an mRNA results in a region that binds toa ribosome to initiate the translation process. During translation, the CDSregion on the mRNA is used as instructions following the genetic code toselect the amino acids to use to construct a protein . Proteins are a funda-mental component for almost all molecular functions within a cell. Proteinscan also bind to promoters to activate or repress transcription, i.e., increasingor decreasing the associated promoters binding affinity to RNAP.The motivating example used in this chapter is a genetic circuit for atoggle switch [14] shown in Figure 1. This genetic circuit is constructed fromtwo transcriptional units . The one on the left begins with the promoter P tet (shown as a bent arrow), followed by its RBS (shown as a half circle), a CDSthat codes for the protein LacI, and finally a terminator (shown as a (cid:62) ). Theone on the right begins with the P lac promoter, which initiates transcriptionof the CDSs for the TetR protein and the green fluorescent protein (GFP).GFP is a reporter, since the cells glow green when it is present. The switchlike behavior is created by mutual repression. Namely, the TetR protein bindsto P tet to repress LacI production, while the LacI protein binds to P lac torepress TetR production. The state of the switch is changed by adding smallmolecule chemical inducers . Namely, when the switch is OFF (i.e., LacI ispresent but no TetR or GFP is present), IPTG can be added, which binds toLacI forming the complex C1, which is unable to repress P lac . This situationleads to TetR and GFP being produced, which represses LacI production andthus changes the switch to the ON state. To change back to the OFF state,aTc can be added, which binds to TetR to form the complex C2, which isunable to repress P tet . This situation leads to LacI being produced, whichrepresses further production of TetR and GFP and thus the changes thegenetic toggle switch to the OFF state.One possible reaction-based model of the genetic toggle switch is shown inFigure 2. This model is derived from a more detailed model, using quasi-steady-state approximations and reaction-based abstractions as describedin [24, 35]. This model is composed of a species for each protein (i.e., LacI,TetR, and GFP) and each small molecule (i.e., IPTG and aTc). This modelalso includes a production reaction for each promoter, P tet and P lac , and adegradation reaction for each protein. The reactions are shown as boxes inthe diagram, with their propensity functions shown inside the boxes. The pa-rameters for these propensity functions are given in Table 1. Note that these Neupane et al. P tet tetRlacI GFP gfp
TetR
LacI P lac C1IPTG
C2aTc
Fig. 1: The genetic toggle switch. This switch is created using two repressors,LacI and TetR, which repress each others production, denoted by the ⊥ and (cid:62) arrows on promoters P tet and P lac . The small molecule IPTG can bind toLacI, effectively reducing LacI’s ability to repress TetR and GFP production.Similarly, the small molecule aTc can bind to TetR to reduce TetR’s ability torepress LacI’s production. To indicate the ON and OFF states of this switch,this circuit includes the reporter protein GFP to cause the cell to glow greenwhen it is present.are simply reasonable default parameters and not measured experimentally,and they can be easily updated if better information becomes available. Theedges are labeled to indicate reactants (r), species consumed by the reactions, products (p), species produced by the reactions, and modifiers (m), speciesneither produced or consumed. The stoichiometry , the number of moleculesproduced or consumed, for each reaction is assumed to be 1, unless indicatedotherwise (e.g., production reactions produce np molecules). The degradationreactions have a propensity that is just the degradation rate, k d , times thecurrent number of molecules of the species that is degrading. The productionreactions have a propensity that is the number of molecules produced, np ,times the rate of production, k p , times the proportion of promoters bound toRNAP in steady-state. This proportion is a function of the amount of repres-sor molecules present in free form (i.e., not bound to the corresponding smallmolecule inducer). Further details are outside the scope of this chapter, butthey can be found in [24, 35].Unlike an electronic circuit, the behavior of a genetic toggle switch circuitis extremely noisy due to the small molecule counts involved. It is, therefore,necessary to evaluate a genetic circuit’s behaviors using stochastic analyses. pproximation Techniques for Stochastic Analysis of Biological Systems 9 k d | LacI | (cid:79) (cid:79) r IPTG m LacI mn p k p | P lac | K o | RNAP | K o | RNAP | + (cid:16) K r | LacI | Kc | IPTG | (cid:17) nc np,p (cid:15) (cid:15) np,p (cid:40) (cid:40) n p k p | P tet | K o | RNAP | K o | RNAP | + (cid:16) K r | TetR | Kc | aTc | (cid:17) nc np,p (cid:104) (cid:104) GFP TetR m aTc mk d | GFP | (cid:15) (cid:15) r k d | TetR | (cid:15) (cid:15) r Fig. 2: Reaction graph adapted from [29] for the genetic toggle switch afterapplying reaction-based abstractions to the chemical reaction network.Table 1: List of parameters for the genetic toggle switch model.
Parameter Symbol Value UnitsDegradation rate k d − Complex formation equilibrium K c − Stoichiometry of binding n c K r − RNAP binding equilibrium K o − Open complex production rate k p − Stoichiometry of production np
10 dimensionlessNumber of RNAP molecules | RNAP |
30 moleculesNumber of P tet promoters | P tet | P lac promoters | P lac | Figure 3 shows the average output response of 100 stochastic simulation runsusing Gillespies stochastic simulation algorithm (SSA) [15]. These simula-tions start with the same initial state with 60 LacI molecules, and 0 for otherspecies. At time 5,000, 100 molecules of IPTG are applied, which activatesthe production of TetR and GFP to bring them to the high state, and re-presses LacI to slow down its production to allow its degradation to reduceits molecule count. When the input IPTG is removed at time 10,000 makingboth inputs absent, the outputs retain their current states. At time 15,000,applying inducer aTc causes the circuit to switch output states again. Re-moving aTc at time 20,000, once again leaves the outputs to hold their states.
It should be noted that this figure shows the average output responses of 100simulation runs, as an individual run may fail to exhibit meaningful logicalbehavior due to the noise in the circuit. This chapter aims to efficiently de-termine the probability of erroneous behavior induced by the inherent noisynature of genetic circuits.
Algorithms 1 to 3 describe the state space approximation procedures for agiven SCK model M = (cid:104) S , R , x (cid:105) with reaction-based abstractions. Notethat models with reaction-based abstractions utilize quasi steady-state ap-proximations [37], where extremely fast reactions are approximated as param-eters on propensity functions to prevent starvation of other slower reactionsduring stochastic analysis. The presented state space approximation methodassumes that probability mass is distributed on a finite and relatively smallnumber of states, and the probability mass does not distribute uniformly astime progresses.With a given SCK model M = (cid:104) S , R , x (cid:105) , state space generation startsby assigning the sole initial state x a 1 to the termination indicator ˆ κ , asshown in Algorithm 1. The termination indicator is a function ˆ κ : X → R + ,which indicates whether state search should terminate from a state onwards.The initial state graph G includes the initial state x as its set of states. Thesubsequent state graphs are then constructed and refined by Algorithm 2. Ingeneral, both state graphs G k − and G k are constructed based on the sameSCK model M and refined from the same initial state x . The differenceis that G k refines ˆ κ values for some explored states in G k − , which mayexpand G k − to include new states in G k . This process of expansion andrefinement is repeated until the size of the approximate state graph stabilizes,at which point an absorbing state is added to this state graph by Algorithm 3.Algorithm 1 terminates by returning the approximated state graph G k . Algorithm 1:
Construction of approximate state space.
Input:
An SCK model M = (cid:104) S , R , x (cid:105) . Output:
Approximated state graph G k = (cid:104) X k , δ k , x (cid:105) . G = (cid:104) X , δ , x (cid:105) , where X = { x } , δ = ∅ ; ˆ κ ( x ) := 1 , ˆ γ ( x ) := 0; k := 0; repeat k := k + 1; Construct finite state graph G k = (cid:104) X k , δ k , x (cid:105) for M using Algorithm 2. until |G k | = |G k − | ; Update G k by adding an an extra absorbing state x abs using Algorithm 3.pproximation Techniques for Stochastic Analysis of Biological Systems 11 F i g . : T h e a v e r ag e o f s t o c h a s t i c s i m u l a t i o n r un s o f t h e g e n e t i c t ogg l e s w i t c h c i r c u i t . T h e L a c I m o l ec u l ec o un t d r o p s t o l o w s oo n a f t e r t h e i n t r o du c t i o n o f I P T G a tt i m e , s , w h i c h i s t h e n r e m o v e d a tt i m e , s . I t s m o l ec u l ec o un t s h a r p l y r i s e s t o h i g h w h e n a T c i s a dd e d a tt i m e , s , w h i c h i s t h e n r e m o v e d a tt i m e , s ,l e a v i n g L a c I t o s t a y a t a h i g h m o l ec u l ec o un t . T e t R h a s t h e o pp o s i t e b e h a v i o r , w h i c h i s c l o s e l y f o ll o w e db y G F P . Algorithm 2 constructs the approximate finite state space based on a user-defined termination indicator κ . Starting with the initial state x , all possiblereactions are scheduled to be explored (line 6). For each such reaction R i , itsupdated state x (cid:48) is obtained by adding the state-change vector v i specifiedby R i to the current state x (line 7). It should be noted that since thestate search in each iteration k begins at the same initial state x , x (cid:48) maynot necessarily be a new state after this step. The termination indicatorvalue at the current state x is then compared against κ to determine if stateexploration should continue (line 11). If the former is lower, then x becomesa (partially) terminal state, whose state expansion only includes its outgoingtransitions leading to existing states in the current state set X k , but omitstransitions leading to states not in X k . Therefore, if x (cid:48) already exists in X k (line 12), the algorithm includes the new state-transition relation ( x , R i , x (cid:48) )and updates its termination indicator (lines 13 to 15). For every state x (cid:48) to beupdated, its predecessor set is constructed (line 14). Each element of this setis a pair of the predecessor state x and the reaction index i , in which a uniqueexisting state transition ( x , R i , x (cid:48) ) defines reachability of x (cid:48) from x throughreaction R i . Then the updated termination indicator ˆ γ ( x (cid:48) ) is determined byline 15. It should be noted that the updated termination indicator ˆ γ is notused to update termination indicator values for other states explored in thecurrent iteration k , and only becomes available at the end of the currentiteration, at which point it is assigned to the current termination indicatorˆ κ (line 26). For each predecessor state x of x (cid:48) , its contribution to ˆ γ ( x (cid:48) ) isthe product of its current state termination value ˆ κ ( x ) and the probability oftransitioning from x to x (cid:48) , defined as the ratio of propensity α i , evaluated at x , to the sum of all propensities at this state. Intuitively, ˆ γ ( x (cid:48) ) accumulatespath probabilities from all of its predecessor states that have been explored initeration k . On line 16, the function explored ( x (cid:48) , k ) checks whether state x (cid:48) has been either expanded or updated at the current iteration k . This state canbe a state discovered at the current iteration or any of the previous iterations.This state is only scheduled to be explored, if it has not been explored yet inthe current iteration. For the case where ˆ κ ( x ) (cid:62) κ , in addition to updatingthe state-transition relation and termination indicator for x (cid:48) , the algorithmincludes it in the state set X k (line 22). This is because x cannot be theterminal state due to its large termination indicator value, and therefore itssuccessor x (cid:48) becomes the potential candidate for a terminal state. This stateis then scheduled for exploration, if the current iteration has not explored it.The termination indicator update is performed every time a new incomingpath is added to a state. It is crucial to have frequent updates since a newincoming path can add its probability contribution to the state, potentiallybringing the termination indicator value above κ , which in turn changesa terminal state to be non-terminal. This update, therefore, guarantees toexplore a state with many incoming paths whose accumulative probabilitiesare significant, although each individual one might be low compared to κ . pproximation Techniques for Stochastic Analysis of Biological Systems 13 Algorithm 2:
State space construction and approximation usingbreadth-first search.
Input:
An approximated global state graph G k − = (cid:104) X k − , δ k − , x (cid:105) . Output:
Updated state graph G k = (cid:104) X k , δ k , x (cid:105) . X k := X k − ; δ k := δ k − ; Enqueue ( queue, x ); while queue (cid:54) = ∅ do x := Dequeue ( queue ); forall i ∈ { j | α j ( x ) > } do Determine the state after reaction R i : x (cid:48) := x + v i ; if x (cid:48) / ∈ X k then ˆ κ ( x (cid:48) ) := 0; ˆ γ ( x (cid:48) ) := 0; if ˆ κ ( x ) < κ then if x (cid:48) ∈ X k then δ k := δ k ∪ { ( x , R i , x (cid:48) ) } ; Pre ( x (cid:48) ) := { ( x , i ) | ( x , R i , x (cid:48) ) ∈ δ k , ∀ i ∈ (1 , · · · , m ) } ; ˆ γ ( x (cid:48) ) := (cid:80) ( x ,i ) ∈ Pre ( x (cid:48) ) (cid:18) ˆ κ ( x ) · α i ( x ) (cid:80) mj =1 α j ( x ) (cid:19) ; if ¬ explored ( x (cid:48) , k ) then Enqueue ( queue, x (cid:48) ); else δ k := δ k ∪ { ( x , R i , x (cid:48) ) } ; Pre ( x (cid:48) ) := { ( x , i ) | ( x , R i , x (cid:48) ) ∈ δ k , ∀ i ∈ (1 , · · · , m ) } ; ˆ γ ( x (cid:48) ) := (cid:80) ( x ,i ) ∈ Pre ( x (cid:48) ) (cid:18) ˆ κ ( x ) · α i ( x ) (cid:80) mj =1 α j ( x ) (cid:19) ; X k := X k ∪ { x (cid:48) } ; if ¬ explored ( x (cid:48) , k ) then Enqueue ( queue, x (cid:48) ); forall x ∈ X k do ˆ κ ( x ) := ˆ γ ( x ); The theoretical state space for the genetic toggle switch described in Sec-tion 4 is infinite. To analyze the model, the state space is truncated based onthe value of κ . This truncation, however, leads to probability leakage (i.e.,cumulative probabilities of reaching states not included in the explored statespace) during the CTMC analysis. To account for probability loss, an ab-sorbing state x abs is created as the sole successor state for all terminal stateson each truncated path, and is added by Algorithm 3 to the state space gen-erated by Algorithm 1. For all states in the global state set, each possiblereaction for state x are checked for exploration. For each reaction R i , if ithas not been explored, its updated state x (cid:48) is set to x abs (line 5). It is obvi-ous that all unexplored transitions from such a terminal state x lead to theabsorbing state. Algorithm 3:
Absorbing state update from approximated global stategraph.
Input:
An approximated global state graph G . Output:
Updated state graph G with an absorbing state x abs . X := X ∪ { x abs } ; forall x ∈ X do forall i ∈ { j | α j ( x ) > } do Determine the state after reaction R i : x (cid:48) := x + v i ; if ( x , R i , x (cid:48) ) / ∈ δ then δ := δ ∪ { ( x , R i , x abs ) } ; The state graph returned by Algorithm 1 is essentially a (sparse) repre-sentation of the transition rate matrix. A standard CTMC analysis can beapplied directly to it to compute the approximate probability distribution.It should be noted that the termination indicator value for each state is onlyused to determine terminal states, and is omitted for the CTMC analysis.With the addition of the absorbing state, the CTMC analysis provides aprobability bound [ l, u ], where 0 (cid:54) l < u (cid:54)
1, and ( u − l ) is the probabilityaccumulated in x abs . Assuming the actual probability to satisfy a CSL prop-erty φ is p , then it holds that l (cid:54) p (cid:54) u . Because the lower bound l doesnot account for probabilities from paths that, if were not truncated, wouldfeed probabilities back to the explored states, as is the case for calculating p .For the upper bound u , it is always greater or equal to p . Because u includesprobabilities accumulated by the absorbing state, of which probabilities fromtruncated paths that would lead to falsification of φ are counted, in additionto probabilities of those leading to the satisfaction of φ . Complexity : The size of generated state space models depends on the dis-tribution of probability over states and the termination threshold. Therefore,detailed characterization of state space complexity is challenging. Intuitively,the state space complexity increases as the termination threshold decreases.This is because a lower termination threshold would allow exploration ofstates with lower accumulated path probabilities that would otherwise be ig-nored with a higher termination threshold. Exploration of these states wouldlikely lead to other new states. Moreover, if the majority of probability isdistributed over a small number of states, a smaller number of states may beexplored compared with a more even probability distribution.The complexity is highly dependent on the user provided termination in-dicator κ . Determining a reasonable value of κ can be an iterative process.Initially, κ can be set to any value that satisfies 0 < κ <<
1, and a stategraph and probability bound [ l, u ] can be generated. The user can then de-crease the value of κ , if necessary, to tighten the probability bound window.The user can repeat the process until the probability bound returned is guar-anteed to prove or disprove the given CSL property. pproximation Techniques for Stochastic Analysis of Biological Systems 15 The presented algorithms in Section 5 are guaranteed to terminate undercertain conditions. This section provides a description of the terminationconditions for each algorithm, and presents a proof for termination.To facilitate the following proof, we first define finite paths of a stategraph and depth for breadth-first search. A finite path ρ of a state graph isa sequence x R −−→ x R −−→ . . . x n − R n − −−−→ x n such that for every 0 (cid:54) i < n ,( x i , R i , x i +1 ) ∈ δ holds for some R i . State x n is reachable in G if x n is reach-able from the initial state through a finite path included in G . Denote the setof all states with depth ı as ı X . At depth 0, X = { x } . At depth ı > ı X isobtained by collecting all newly created states resulted from the one-step BFSsearch on all states in ı − X . Therefore, the depth for a state is determinedwhen it is explored for the first time. Note that X ∩ X · · · ı − X ∩ ı X = ∅ .Termination condition for Algorithm 1 requires that, as both the depth ı and iteration k increase, the sum of termination indicator values for allstates of ı X k decreases, with possibly finitely many iterations where thissum remains constant. This is formulated by Theorem 1 below. Theorem 1 (Termination of Algorithm 1).
Algorithm 1 terminates aftera finite number of iterations with a given κ , where < κ << , if thestate graph G +1 satisfies the following condition: for each depth > , theremust exist depth (cid:54) ı (cid:54) such that x d R h −−→ x d +1 R i −→ . . . x d +( m − R l −→ x d + m is a finite path in G +1 , and x d ∈ ı X +1 , x d +( m − ∈ X +1 , x d + m ∈ X +1 ∪ X +1 ∪ · · · − X +1 ∪ X +1 , and m ∈ Z (cid:62) .Proof. Initially, X = { x } and ˆ κ ( x ) = 1. At iteration k = 1, during theconstruction of G (line 6 of Algorithm 1), each state at depth 1, x ∈ X , isdiscovered for the first time when x is explored (line 6 to 24 of Algorithm 2).Therefore, the current termination indicator ˆ κ ( x ) is assigned a 0 (line 9of Algorithm 2), but its next termination indicator ˆ γ ( x ) gets updated byˆ κ ( x ), so that 0 < ˆ γ ( x ) (cid:54)
1. Each new state x ∈ X generated from X , is ignored, since ˆ κ ( x ) = 0, which is less than κ , and x / ∈ X (line 11to 16 of Algorithm 2). Then at iteration k = 2, the sum of terminationindicators is ζ = (cid:80) x ∈ X ˆ κ ( x ), where each ˆ κ ( x ) is a fraction of ζ , and ζ = ˆ κ ( x ) = 1. Therefore, ζ is solely contributed from ζ . If a self-loop transition { x , R , x } exists, then ζ > ζ ; otherwise ζ = ζ .Therefore, ζ (cid:62) ζ . Similar to the previous iteration, the updated ˆ γ ( x )will be used in the next iteration.In general, state set ı X at depth ı is first obtained in iteration ı by col-lecting all the new states, i.e., states whose depth has not been determined,which are expanded from states in ı − X . The sum of all termination indica-tor values for states in ı X is calculated at iteration ı + 1 by either line 15or 21 of Algorithm 2. To differentiate the termination indicator function ˆ κ in different iterations, we denote ˆ κ ı ( x ) as the termination indicator value for state x at iteration ı . The sum of all termination indicators at iteration ı + 1is computed as follows: ı ζ ı +1 = (cid:88) x (cid:48) ∈ ı X ı +1 ˆ κ ı +1 ( x (cid:48) )= (cid:88) x (cid:48) ∈ ı X ı +1 (cid:88) ( x ,i ) ∈ Pre ( x (cid:48) ) (cid:32) ˆ κ ı ( x ) · α i ( x ) (cid:80) mj =1 α j ( x ) (cid:33) If (cid:83) x (cid:48) ∈ ı X ı +1 Pre ( x (cid:48) ) is equal to all transition firings of every state in ı − X ı ,then termination indicator values for all the states at depth ı − ı , and hence ı ζ ı +1 = ı − ζ ı . On the other hand, if there exists one ormore transition firings from ı − X ı to depth other than ı , then ı ζ ı +1 < ı − ζ ı .Therefore, ı − ζ ı (cid:62) ı ζ ı +1 .We can, therefore, establish the following conclusion:1 = ζ (cid:62) ζ (cid:62) · · · ı − ζ ı (cid:62) ı ζ ı +1 · · · − ζ (cid:62) ζ +1 From the termination condition stated in Theorem 1, the slowest termina-tion scenario, i.e., the maximal number of iterations required to terminateAlgorithm 1, is the following:1 = ζ = ζ = · · · = ı ζ ı +1 = · · · = − ζ > ζ +1 . The inequality − ζ > ζ +1 holds only if at least one state in − X exe-cutes a transition leading to a state in X +1 ∪ X +1 ∪ · · · − X +1 , but notin X +1 . State x d + m in Theorem 1 is such a state. Additionally, the termi-nation condition requires that at least ı ζ ı +1 = · · · = − ζ > ζ +1 holdsfor every depth . This requirement guarantees that the sum of terminationindicator values keeps decreasing, with possibly many (or zero) iterationswhere this sum remains unchanged. Therefore, after a finite number ξ of it-erations, ξ − ζ ξ < κ . Since ξ − ζ ξ is the sum of all individual terminationindicator values, in the next iteration ( ξ + 1), termination indicator ˆ κ ( ξ x ) isless than κ for all states in ξ X ξ +1 , and they become terminal states. Hence, |G ξ | = |G ξ +1 | .Finally, Algorithm 3 terminates, provided its input state graph generatedby Algorithm 2 is finite. This has been proven true, and hence Algorithm 3always terminates. Therefore, Theorem 1 is true. (cid:117)(cid:116) pproximation Techniques for Stochastic Analysis of Biological Systems 17 Algorithms 1, 2 and 3 were implemented in Java as a prototype tool withinthe iBioSim genetic design automation (GDA) tool [36, 30, 39]. This toolconstructs an approximate CTMC transition rate matrix, which is then an-alyzed with the help of the PRISM probabilistic model checking tool [27].Experiments were performed on a 3.2 GHz AMD Debian Linux PC with sixcores and 64 GB of RAM. The presented CTMC approximation method wasevaluated on several CSL properties for the genetic toggle switch describedin Section 4. This method is then applied to several benchmark examples,and the results are compared with those generated by the STAR tool.
An important metric for a toggle switch circuit is the response time. In thefirst set of experiments, the goal is to determine the genetic toggle switch’sresponse time (i.e., the time it takes to switch from the OFF state to theON state). The initial OFF state for the toggle switch has 60 LacI, 0 TetR,and 100 IPTG molecules, representing the circuit has just received the setinput to switch to the ON state. It should be noted that the input valueof 100 molecules is chosen to ensure that the circuit should switch to theON state, but any moderately large value of input could be used as IPTG isrepresented as a boundary condition species which means that its moleculecount is treated as non-depleting and that it is not consumed by any reactionsthat it occurs in. The CSL property, F ( t (cid:54) , LacI < ∧ TetR > ,
100 seconds (an approximation of the cell cycle in
E. coli [40]). The ON state is characterized by LacI dropping below 20 andTetR rising above 40 molecules.The termination indicator values are set to 10 − , 10 − , 10 − and 10 − .Approximate state space generation and CTMC analysis are performed foreach such value. In addition, intermediate verification results are generated ona time course from 0 to 2 ,
100 seconds with the increment of 100 seconds. Tomeasure the accuracy of the presented state space approximations with differ-ent termination indicator values, a reference finite-state SCK model is createdallowing both LacI and TetR to reach the upper bound of 300 molecules each,which is significantly higher than the upper bounds of TetR and LacI for allexperiments performed. The reference model, therefore, incurs significantlylarger state space with 90 ,
601 states, but provides accurate verification re-sults for comparison.Both the accuracy and performance results for the response rate verifi-cation are presented in Table 2. The column “ |G | ” lists results for the ap-proximate state graph size used for the CTMC analysis, respectively. The column “ (cid:15) ” reports the difference between minimum ( P min ) and maximum( P max ) final response rate probability, which can be taken as the uncertaintywindow. The columns labeled T P min and T P max , provide the CTMC analy-sis time taken by the PRISM tool to calculate the minimum and maximumprobability value, respectively.Table 2: Genetic toggle switch response rate results. κ |G| P min P max (cid:15) T P min T P max ref 90601 0 . − − . − − . . . × − .
492 0 . − . . . × − .
714 0 . − . . . × − .
811 0 . − . . . × − .
161 1 . As the table shows, reducing the termination indicator value improves theaccuracy of the final probability, at the price of increased performance cost.Furthermore, the final probability for t (cid:54) κ . As we decreasethe value for κ , from 10 − to 10 − , the error window becomes narrower from1 . × − to 5 . × − . The reference model has the final probabilityof 0 . κ = 10 − , its final probability already produces veryaccurate final probability with significantly smaller performance cost. The ap-proximate model only explores 6333 states, compared to 90 ,
601 states fromthe reference model, but it produces the minimum result 0 . . × − . The runtime for the CTMC analysison the reference model is 23 .
49 seconds, much longer than the runtime foranalyzing the approximate CTMC. As an additional comparison, the sametoggle switch model built with pre-determined thresholds of molecule countsfor LacI and TetR in [31] produces a state graph of 70 states, and CTMCanalysis with the same initial condition and CSL property reports a finalprobability of 98 . ,
100 seconds. Thisbehavior occurs if production of LacI erroneously and significantly inhibitsTetR’s production to let TetR degrade away and consequently switch state.The toggle switch is initialized to a starting state with 60 LacI molecules,and 0 molecules for all other species. The same CSL properties are verified pproximation Techniques for Stochastic Analysis of Biological Systems 19 and the results are summarized in Table 3. Similar to the above experiment,the final probability for t (cid:54) κ . Decreasing thevalue for κ from 10 − to 10 − decreases the error window from 4 . × − to 1 . × − . Figure 4b shows the time-series plot for the genetic toggleswitch failure rates with κ = 10 − and κ = 10 − .Table 3: Genetic toggle switch failure rate results. κ |G| P min P max (cid:15) T P min T P max ref 90601 0 . − − . − − . . . × − .
239 0 . − . . . × − .
287 0 . − . . . × − .
361 0 . − . . . × − .
560 0 . To illustrate the accuracy and efficiency of the presented method, we com-pared the probability distribution results with the STAR tool for the birth-death model and the presented toggle switch model. Table 4 summarizes thecomparison for a simple birth-death model, whose birth rate is 1 and deathrate is 0.1. Column “ t ” shows the time point at which the state probabilityis computed, and column labeled (cid:15) max shows the maximum absolute prob-ability difference for the same individual state obtained from the two tools,among all explored states. Columns labeled T iBioSim and T PRISM list runtimesin seconds to generate the state space in iBioSim and to analyze the modelin PRISM for each given time point, respectively. Column T STAR lists theruntime reported by STAR. The maximum probability difference reaches itspeak value of 2 . × − at time point t = 50. All other time points showsignificantly smaller errors. The run time to analyze the model in iBioSim and PRISM is less than a second as the generated state space is only 28states. The STAR tool also reports a similar run time.Table 5 shows a comparison of results for the aforementioned toggle switch.Our proposed method produces accurate results compared to those from theSTAR tool, as is indicated by the maximal probability difference ( (cid:15) max ).Columns T iBioSim and T PRISM list runtimes in seconds to generate the statespace in iBioSim and to analyze the model in PRISM for each given timepoint, respectively. Column labeled T STAR lists the runtime reported by STAR.The combined runtime to generate the state space and analyze the model for P e r c e n t Prob(Except being in Absorbing-State)Full Model(State-Space="300x300")Max Prob(Prob + p(Being in Absorbing-State)) (a) Genetic toggle switch failure rate with κ = 10 − P e r c e n t Prob(Except being in Absorbing-State)Full Model(State-Space="300x300")Max Prob(Prob + p(Being in Absorbing-State)) (b) Genetic toggle switch failure rate with κ = 10 − Fig. 4: Error window comparison for different values of κ .our method is less than 24 seconds for different time points and remainsalmost constant as the time point t increases. The STAR tool reports shorterruntime for smaller t but linear increase in runtime as the time point valuegets larger. pproximation Techniques for Stochastic Analysis of Biological Systems 21 Table 4: State probability comparison for birth-death model with κ = 10 − . t (cid:15) max T iBioSim T PRISM T STAR
10 1 . × − .
06 0 .
304 0 . . × − .
06 0 .
311 0 . . × − .
06 0 .
303 0 . . × − .
06 0 .
307 0 . . × − .
06 0 .
309 0 . Table 5: State probability comparison for switching rate experiment of toggle-switch model with κ = 10 − . t (cid:15) max T iBioSim T PRISM T STAR
400 1 . × − .
86 4 .
44 7 . . × − .
86 4 .
51 17 . . × − .
86 4 .
59 29 . . × − .
86 4 .
69 40 . . × − .
86 4 .
79 50 . This chapter presents a method that builds an approximate state space of ge-netic circuit models to analyze infinite-state continuous-time Markov chains.This approximation method iteratively expands from the initial state usinga breadth first search, computes and updates the termination indicator valuefor each state on-the-fly, based on the cumulative path probabilities for allincoming transitions to a state. The probability of each path segment is theratio of the propensity of a reaction to the sum of all propensities for anygiven state. Our state space approximation is determined by comparing thestate termination indicator to a user-provided termination threshold and onlyexploring states with a significant termination indicator value. This methodis capable of computing the approximate state space with no prior knowledgeand is completely automated.For future work, we plan to improve and optimize probability approxima-tion for re-convergent paths that close cycles during the state exploration inorder to achieve potentially faster termination of the state search. We willconsider different approaches to determining the termination indicator valueautomatically from the CSL property being analyzed. Additionally, we planto explore augmenting our technique with bi-simulation minimization andabstraction to further minimize the generated state space and better allowfor scalability. To improve performance of tool implementation, we plan toinvestigate tighter integration with the PRISM tool.
The authors thank Verena Wolf for providing benchmarks and the STARtool. The authors would also like to thank Dave Parker and Joachim Kleinfor providing assistance in interfacing with the PRISM tool. We also thankthe reviewers for their feedback on an earlier version of this paper. pproximation Techniques for Stochastic Analysis of Biological Systems 23
References [1] Abate A, Brim L, ˇCeˇska M, Kwiatkowska M (2015) Adaptive aggrega-tion of Markov chains: Quantitative analysis of chemical reaction net-works. In: Kroening D, P˘as˘areanu CS (eds) Computer Aided Verifica-tion, Springer International Publishing, Cham, pp 195–213[2] Aziz A, Sanwal K, Singhal V, Brayton R (2000) Model-checkingcontinuous-time Markov chains. ACM Trans Comput Logic 1:162–170[3] Baier C, Gr¨oßer M, Ciesinski F (2004) Partial order reduction forprobabilistic systems. In: 1st International Conference on Quan-titative Evaluation of Systems (QEST 2004), 27-30 September2004, Enschede, The Netherlands, IEEE Computer Society, pp 230–239, DOI 10.1109/QEST.2004.1348037, URL https://doi.org/10.1109/QEST.2004.1348037 [4] Baier C, D’Argenio P, Groesser M (2006) Partial order reductionfor probabilistic branching time. Electron Notes Theor Comput Sci153(2):97–116, DOI 10.1016/j.entcs.2005.10.034, URL http://dx.doi.org/10.1016/j.entcs.2005.10.034 [5] Burrage K, Hegland M, Macnamara S, Sidje R (2006) A Krylov-basedfinite state projection algorithm for solving the chemical master equa-tion arising in the discrete modelling of biological systems. In: LangvilleAN, Stewart WJ (eds) MAM 2006 : Markov Anniversary Meeting: aninternational conference to celebrate the 150th anniversary of the birthof A.A. Markov, Boston Books, Charleston, South Carolina, pp 21–38,URL http://eprints.qut.edu.au/46148/ [6] Clarke EM, Emerson EA, Sistla AP (1986) Automatic verification offinite-state concurrent systems using temporal logic specifications. ACMTrans Program Lang Syst 8(2):244–263, DOI 10.1145/5397.5399, URL http://doi.acm.org/10.1145/5397.5399 [7] Courcoubetis C, Yannakakis M (1988) Verifying temporal properties offinite state probabilistic programs. In: Proc. 29th Annual Symposium onFoundations of Computer Science (FOCS’88), IEEE Computer SocietyPress, pp 338–345[8] Courcoubetis C, Yannakakis M (1995) The complexity of probabilisticverification. Journal of the ACM 42(4):857–907[9] D´ıaz ´AF, Baier C, Earle CB, Fredlund L (2012) Static partial orderreduction for probabilistic concurrent systems. In: Ninth InternationalConference on Quantitative Evaluation of Systems, QEST 2012, London,United Kingdom, September 17-20, 2012, IEEE Computer Society, pp104–113, DOI 10.1109/QEST.2012.22, URL https://doi.org/10.1109/QEST.2012.22 [10] Fecher H, Leucker M, Wolf V (2006) Don’t know in probabilistic sys-tems. In: Proceedings of the 13th international conference on ModelChecking Software, Springer-Verlag, Berlin, Heidelberg, SPIN’06, pp 71–
88, DOI 10.1007/11691617 5, URL http://dx.doi.org/10.1007/11691617_5 [11] Fisler K, Vardi MY (1998) Bisimulation minimization in an automata-theoretic verification framework. In: FMCAD, Springer, Lecture Notesin Computer Science, vol 1522, pp 115–132[12] Fisler K, Vardi MY (1999) Bisimulation and model checking. In: PierreL, Kropf T (eds) Correct Hardware Design and Verification Methods,Springer Berlin Heidelberg, Berlin, Heidelberg, pp 338–342[13] Fisler K, Vardi MY (2002) Bisimulation minimization and symbolicmodel checking. Formal Methods in System Design 21(1):39–78, DOI10.1023/A:1016091902809, URL http://dx.doi.org/10.1023/A:1016091902809 [14] Gardner TS, Cantor CR, Collins JJ (2000) Construction of a genetictoggle switch in
Escherichia coli . Nature 403:339–342[15] Gillespie DT (1977) Exact stochastic simulation of coupled chemicalreactions. J Chem Phys 81(25):2340–2361[16] Grassmann W (1977) Transient solutions in markovian queue-ing systems. Computers & Operations Research 4(1):47 –53, DOI http://dx.doi.org/10.1016/0305-0548(77)90007-7, URL [17] Gross D, Miller DR (1984) The randomization technique as a model-ing tool and solution procedure for transient markov processes. OperRes 32(2):343–361, DOI 10.1287/opre.32.2.343, URL http://dx.doi.org/10.1287/opre.32.2.343 [18] Hansson H, Jonsson B (1994) A logic for reasoning about time and reli-ability. Formal Aspects of Computing 6(5):512–535[19] Heath J, Kwiatkowska M, Norman G, Parker D, Tymchyshyn O (2006)Probabilistic model checking of complex biological pathways. In: Pri-ami C (ed) Computational Methods in Systems Biology, Springer BerlinHeidelberg, Berlin, Heidelberg, pp 32–47[20] Henzinger TA, Mateescu M, Wolf V (2009) Sliding window abstractionfor infinite Markov chains. In: Bouajjani A, Maler O (eds) ComputerAided Verification, Springer Berlin Heidelberg, Berlin, Heidelberg, pp337–352[21] Hermanns H, Wachter B, Zhang L (2008) Probabilistic CEGAR. In:Proceedings of the 20th international conference on Computer AidedVerification, Springer-Verlag, Berlin, Heidelberg, CAV ’08, pp 162–175,DOI 10.1007/978-3-540-70545-1 16, URL http://dx.doi.org/10.1007/978-3-540-70545-1_16 [22] Katoen JP, Kemna T, Zapreev I, Jansen DN (2007) Bisimulation min-imisation mostly speeds up probabilistic model checking. In: Proceedingsof the 13th international conference on Tools and algorithms for the con-struction and analysis of systems, Springer-Verlag, Berlin, Heidelberg, pproximation Techniques for Stochastic Analysis of Biological Systems 25
TACAS’07, pp 87–101, URL http://dl.acm.org/citation.cfm?id=1763507.1763519 [23] Katoen JP, Klink D, Leucker M, Wolf V (2007) Three-valued abstrac-tion for continuous-time Markov chains. In: Proceedings of the 19th in-ternational conference on Computer aided verification, Springer-Verlag,Berlin, Heidelberg, CAV’07, pp 311–324, URL http://dl.acm.org/citation.cfm?id=1770351.1770401 [24] Kuwahara H, Myers C, Barker N, Samoilov M, Arkin A (2006) Auto-mated abstraction methodology for genetic regulatory networks. TransComp Syst Biol VI:150–175[25] Kwiatkowska M, Norman G, Parker D (2007) Stochastic model checking.In: Bernardo M, Hillston J (eds) Formal Methods for the Design of Com-puter, Communication and Software Systems: Performance Evaluation(SFM’07), Springer, LNCS (Tutorial Volume), vol 4486, pp 220–270[26] Kwiatkowska M, Norman G, Parker D, Qu H (2010) Assume-guaranteeverification for probabilistic systems. In: Esparza J, Majumdar R (eds)Tools and Algorithms for the Construction and Analysis of Systems, pp23–37[27] Kwiatkowska M, Norman G, Parker D (2011) PRISM 4.0: Verificationof probabilistic real-time systems. In: Gopalakrishnan G, Qadeer S (eds)Proc. 23rd International Conference on Computer Aided Verification(CAV’11), Springer, LNCS, vol 6806, pp 585–591[28] Lapin M, Mikeev L, Wolf V (2011) Shave: stochastic hybrid analy-sis of markov population models. In: Proceedings of the 14th ACMInternational Conference on Hybrid Systems: Computation and Con-trol, HSCC 2011, Chicago, IL, USA, April 12-14, 2011, ACM, pp311–312, URL https://publications.cispa.saarland/891/ ,pub id: 705 Bibtex: LaMiWo 11:Shave URL date: None[29] Madsen C (2013) Stochastic analysis of synthetic genetic circuits. PhDthesis, University of Utah[30] Madsen C, Myers CJ, Patterson T, Roehner N, Stevens JT, WinsteadC (2012) Design and test of genetic circuits using iBioSim. IEEE DesignTest of Computers 29(3):32–39[31] Madsen C, Zhang Z, Roehner N, Winstead C, Myers C (2014) Stochas-tic model checking of genetic circuits. J Emerg Technol ComputSyst 11(3):23:1–23:21, DOI 10.1145/2644817, URL http://doi.acm.org/10.1145/2644817 [32] Mikeev L, Sandmann W, Wolf V (2013) Numerical approximation of rareevent probabilities in biochemically reacting systems. In: Gupta A, Hen-zinger TA (eds) Computational Methods in Systems Biology, SpringerBerlin Heidelberg, Berlin, Heidelberg, pp 5–18[33] Munsky B, Khammash M (2006) The finite state projection algo-rithm for the solution of the chemical master equation. J Chem Phys124(4):044104, DOI 10.1063/1.2145882, URL https://doi.org/10.1063/1.2145882 , https://doi.org/10.1063/1.2145882 [34] Munsky B, Khammash M (2008) The finite state projection approach forthe analysis of stochastic noise in gene networks. IEEE Transactions onAutomatic Control 53(Special Issue):201–214, DOI 10.1109/TAC.2007.911361[35] Myers CJ (2009) Engineering Genetic Circuits, 1st edn. Chapman& Hall/CRC Mathematical and Computational Biology, Chapman &Hall/CRC[36] Myers CJ, Barker N, Jones K, Kuwahara H, Madsen C, Nguyen NPD(2009) iBioSim: a tool for the analysis and design of genetic circuits.Bioinformatics 25(21):2848–2849, DOI 10.1093/bioinformatics/btp457[37] Rao CV, Arkin AP (2003) Stochastic chemical kinetics and the quasi-steady-state assumption: Application to the Gillespie algorithm. JChem Phys 118(11):4999–5010, DOI 10.1063/1.1545446, URL https://doi.org/10.1063/1.1545446 , https://doi.org/10.1063/1.1545446 [38] Wachter B, Zhang L, Hermanns H (2007) Probabilistic model checkingmodulo theories. In: Fourth International Conference on the Quantita-tive Evaluation of Systems (QEST 2007), pp 129–140[39] Watanabe L, Nguyen T, Zhang M, Zundel Z, Zhang Z, Mad-sen C, Roehner N, Myers C (0) ibiosim 3: A tool for model-based genetic circuit design. ACS Synthetic Biology 0(0):null,DOI 10.1021/acssynbio.8b00078, URL https://doi.org/10.1021/acssynbio.8b00078 , pMID: 29944839, https://doi.org/10.1021/acssynbio.8b00078 [40] Zheng H, Ho PY, Jiang M, Tang B, Liu W, Li D, Yu X, KlecknerNE, Amir A, Liu C (2016) Interrogating the escherichia coli cell cycleby cell dimension perturbations. Proceedings of the National Academyof Sciences 113(52):15000–15005, DOI 10.1073/pnas.1617932114, URL ,,