- Email: [email protected]
Stochastic analysis of Chemical Reaction Networks using Linear Noise Approximation
Accepted Manuscript Title: Stochastic Analysis of Chemical Reaction Networks Using Linear Noise Approximation Author: Luca Cardelli Marta Kwiatkowska Luca Laurenti PII: DOI: Reference:
S0303-2647(16)30203-9 http://dx.doi.org/doi:10.1016/j.biosystems.2016.09.004 BIO 3702
To appear in:
BioSystems
Received date: Revised date: Accepted date:
27-11-2015 8-7-2016 1-9-2016
Please cite this article as: Luca Cardelli, Marta Kwiatkowska, Luca Laurenti, Stochastic Analysis of Chemical Reaction Networks Using Linear Noise Approximation, (2016), http://dx.doi.org/10.1016/j.biosystems.2016.09.004 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
*Manuscript Click here to view linked References
ip t
Stochastic Analysis of Chemical Reaction Networks Using Linear Noise ApproximationI Luca Cardellia,b,∗, Marta Kwiatkowskaa,∗, Luca Laurentia,∗ Department of Computer Science, University of Oxford b Microsoft Research Cambridge
us
cr
a
Abstract
Ac ce p
te
d
M
an
Stochastic evolution of Chemical Reactions Networks (CRNs) over time is usually analysed through solving the Chemical Master Equation (CME) or performing extensive simulations. Analysing stochasticity is often needed, particularly when some molecules occur in low numbers. Unfortunately, both approaches become infeasible if the system is complex and/or it cannot be ensured that initial populations are small. We develop a probabilistic logic for CRNs that enables stochastic analysis of the evolution of populations of molecular species. We present an approximate model checking algorithm based on the Linear Noise Approximation (LNA) of the CME, whose computational complexity is independent of the population size of each species and polynomial in the number of different species. The algorithm requires the solution of first order polynomial differential equations. We prove that our approach is valid for any CRN close enough to the thermodynamical limit. However, we show on four case studies that it can still provide good approximation even for low molecule counts. Our approach enables rigorous analysis of CRNs that are not analyzable by solving the CME, but are far from the deterministic limit. Moreover, it can be used for a fast approximate stochastic characterization of a CRN. Keywords: Chemical Reaction Networks, Linear Noise Approximation, Probabilistic
I
Authors order is alphabetical. Corresponding author Email addresses: [email protected] (Luca Cardelli), [email protected] (Marta Kwiatkowska), [email protected] (Luca Laurenti) ∗
Preprint submitted to BioSystems
July 8, 2016
Page 1 of 24
ip t
Logic, Model Checking 1. Introduction
Ac ce p
te
d
M
an
us
cr
Chemical reaction networks (CRNs) and mass action kinetics are well studied formalisms for modelling biochemical systems [1]. In recent years, CRNs have also been successfully used as a formal programming language for biochemical systems [2, 3, 4]. There are two well established approaches for analyzing chemical networks: deterministic and stochastic [5]. The deterministic approach models the kinetics of a CRN as a system of ordinary differential equations (ODEs) and represents average behaviour, valid in the thermodynamic limit, when the molecular population is sufficiently high [6]. The stochastic approach, on the other hand, is based on the Chemical Master Equation (CME) and models the CRN as a continuous-time Markov chain (CTMC) [7]. The stochastic behavior can be analyzed by stochastic simulation [5] or by exhaustive probabilistic model checking of the CTMC, which can be performed, for example, by using PRISM [8]. Exhaustive analysis of the CTMC is able to find the best- and worst-case scenarios and is correct for any population size, but suffers from the state-space explosion problem [9] and can only be used for relatively small systems. In contrast, deterministic methods are much more robust with respect to state-space explosion, but unable to represent stochastic fluctuations, which play a fundamental role when the system is not in thermodynamic equilibrium. 1.0.1. Contributions. In this paper we develop a novel approach for analysing the stochastic evolution of a CRN based on the Linear Noise Approximation (LNA) of the CME. We formulate SEL (Stochastic Evolution Logic), a probabilistic logic for CRNs that enables reasoning about probability, expectation and variance of linear combinations of populations of the species. Examples of properties that can be specified in our logic include “the maximum expected population size of species λ1 during the first 20 seconds is 75 molecules” and “the probability that the combined population of species λ1 and λ2 has degraded between 10 and 30 seconds is less than 0.1”. We propose an approximate model checking algorithm for the logic based on the LNA and implement it in Matlab and Java. We demonstrate that the complexity of model checking is polynomial in the initial number of species and independent of the initial molecule counts, thus ameliorating state-space explosion. Further, we 2
Page 2 of 24
an
us
cr
ip t
show that model checking is exact when approaching the thermodynamic limit. Though the algorithm may not be accurate for systems far from the deterministic limit, this generally happens when the populations are small, in which case the analysis can be performed by transient analysis of the induced CTMC [10]. Our approach is essential for CRNs that cannot be analyzed by (partial) state space exploration, because of large or infinite state spaces. Moreover, it is useful for a fast (approximate) stochastic characterization of CRNs, since solving the LNA is much faster than solving the CME [11]. We prove asymptotic correctness of LNA-based model checking and show on four examples that it is still possible to obtain very good approximations even for small population systems, compared to standard uniformisation [10] and statistical model checking implemented in PRISM [8].
Ac ce p
te
d
M
1.0.2. Related work The closest work to ours is by Bortolussi et al. [12], which uses the Central Limit Approximation (CLA) (essentially the same as the LNA) for checking restricted timed automata specifications, assuming a fixed population size. Wolf et al. [13] develop a sliding window method to approximately verify infinite-state CTMCs, which applies to cases where most of the probability mass is concentrated in a confined region of the state space. Recently, Finite State Projection algorithms (FSP algorithms) for the solution or approximation of the CME have been introduced [14]. Sliding window and FSP algorithms apply to the induced CTMC, but require at least partial exploration of the state space, and are thus not immune to state-space explosion. Moment closure techniques [15, 16] improve scalability by estimating the first k ∈ N moments of the distribution of the species over time. The LNA itself can be seen as a moment closure technique, as a Gaussian distribution is completely characterized by the first two moments. However, the LNA tells us more because it guarantees that, if certain conditions are satisfied, the distribution of the process is Gaussian. Continuous Stochastic Logic (CSL), originally introduced in [17] and extended by Baier et al. [18], is a logic widely used to perform model checking of continuous-time Markov chains. CSL combines temporal operators of the logic CTL with the probabilistic and steady-state operators, and is further extended with reward operators in [10]. CSL model checking is based on solving the CME and proceeds through uniformisation of the CTMC, essentially a time discretisation, and thus involves traversal of the full state space. This can be partially ameliorated by fast adaptive uniformisation [19] that does 3
Page 3 of 24
te
d
M
an
us
cr
ip t
not consider states with negligible probability. An alternative is statistical model checking (SMC) which involves A key operator of CSL is probabilistic reachability, that is, computing the probability that a particular region of the state space is reached over a given time interval. Although SEL is endowed with a probabilistic operator, this operator gives the average value of the probability over time and, if the time interval is not a singleton, this is not equivalent to probabilistic reachability. Nevertheless, as shown in [20], SEL and our approximate model checking algorithm can be extended to express reachability, but currently lacks reward operators. The CSL steady-state operator of CSL cannot be added to CSL because LNA is acccurate only for finite time. PRISM [8] implements CSL model checking using uniformisation, fast adaptive uniformisation and statistical model checking. Hybrid Automata Stochastic Logic HASL [21] is an expressive specification formalism for stochastic Petri nets based on linear hybrid automata that is employed by the tool Cosmos [22]. CRNs have a natural interpretation in terms of stochastic Petri nets, see e.g. [23]. The HASL formalism is more expressive than SEL and CSL, and can express CSL probabilistic reachability and expected reward properties. HASL model checking proceeds through statistical model checking of the product of a HASL specification automaton and the Petri net, and is implemented in Cosmos. In contrast, SEL model checking follows through approximating the solution of the CME with the Gaussian process induced by the LNA.
Ac ce p
1.0.3. Structure of the paper In Section 2 we summarise the deterministic and stochastic modelling approaches for CRNs, and in Section 3 we describe the Linear Noise Approximation method. Section 4 introduces the logic SEL and the corresponding model checking algorithm based on the LNA. In Section 5 we demonstrate our approach on four networks taken from the literature. Section 6 concludes the paper. 2. Chemical Reaction Networks
A chemical reaction network (CRN) C = (Λ, R) is a pair of finite sets, where Λ is the set of chemical species and R the set of reactions. |Λ| denotes the size of the set of species. A reaction τ ∈ R is a triple τ = (rτ , pτ , kτ ), where rτ , pτ ∈ N|Λ| and kτ ∈ R>0 . rτ and pτ represent the stoichiometry of reactants and products and kτ is the coefficient associated to the rate of the 4
Page 4 of 24
Ac ce p
te
d
M
an
us
cr
ip t
reaction; its dimension is s−1 . We often write reactions as λ1 + λ3 →k1 2λ2 instead of τ1 = ([1, 0, 1]T , [0, 2, 0]T , k1 ), where ·T indicates the transpose of a vector. We define the net change associated to a reaction τ by υτ = pτ − rτ . For example, for τ1 as above, we have υτ1 = [−1, 2, −1]T . We make the assumption that the system is well stirred, that is, the probability of the next reaction occurring between two molecules is independent of the location of those molecules. We consider fixed volume V and temperature; under these assumptions a configuration or state x ∈ N|Λ| of the system is given by the number of molecules of each species. We define [x] = Nx , the vector of the species concentration in x for a given N , where N = V · NA is the volumetric factor, V is the volume of the solution and NA is Avogadro’s number. The physical dimension of N is M ol−1 · L, where M ol indicates mole and L is litre. Given λi ∈ Λ then #λi x ∈ N represents the number of molecules of λi in x and [λi ] x ∈ R the concentration of λi in the same configuration. In some cases we elide x, and we simply write #λi and [λi ] i . The dimension instead of #λi x and [λi ] x. They are related by [λi ] = #λ N −1 of [λi ] is M ol · L . The propensity αn,τ of a reaction τ in terms of the number of molecules (here subscript n stands for the numebr of molecules) is a function of the current configuration of the system x such that αn,τ (x)dt is the probability that a reaction event occurs in the next infinitesimal interval dt. In this paper we assume as valid the stochastic form of the law of mass action, so the propensity rates are proportional to the number of molecules that participate in the reaction [7]. Stochastic models consider the system in terms of numbers of molecules, while deterministic ones, generally, in terms of concentrations, denoted αc,τ (x) where subscript c stands for concentrations, and the relationship is as follows. For a reaction τ = (rτ , pτ , kτ ), given the configuQ rτ,i ration x and rτ,i , the i-th component of rτ , then αc,τ (x) = kτ |Λ| i=1 ([λi ] x) is the propensity function expressed in terms of concentrations as given by the deterministic law of mass action. It is possible to show that, for any order of reaction, αn,τ (x) ≈ N αc,τ (x) if N is sufficiently large [24]. Note that αc,τ is independent of N . In this paper we are interested only in finite time horizon, because of the problematic character of studying solutions of ODEs for infinite time horizon [25]. Example 2.1. Consider the CRN C = ({λ1 , λ2 , λ3 }, R), where R = {(λ1 + λ2 →10 λ2 + λ2 ), (λ2 + λ3 →10 λ3 + λ3 )}, with initial conditions #λ1 = 98, #λ2 = 1, #λ3 = 1, for a system with N = 1000. Figure 1 plots the ex5
Page 5 of 24
ip t cr us
an
Figure 1: Expected number and standard deviation of species of the CRN of Example 2.1 for the given initial conditions, calculated by simulating the CME.
te
d
M
pectation and standard deviation of population sizes. We may wish to check if the maximum expected value of #λ2 remains smaller than 75 molecules during the first 2sec. However, the system is stochastic, so we also need to analyse whether the variance is limited enough when #λ2 reaches the maximum. Sometimes, analysis of first and second moments does not suffice, so it could be of interest to check the probability of some events, for instance, is the probability that #λ2 − (#λ1 + #λ3 ) > 0, between t1 = 0.5sec and t2 = 1.0sec, greater than 0.6?
Ac ce p
2.1. Deterministic semantics Let C = (Λ, R) be a CRN. The deterministic model approximates the concentration of the species of the system over time as a set of autonomous polynomial first order differential equations: dΦ(t) = F (Φ(t)) dt
(1)
P Q|Λ| where F (Φ(t)) = τ =(rτ ,pτ ,kτ )∈R υτ αc,τ (Φ(t)) and αc,τ (Φ(t)) = kτ i=1 Φi (t)rτ,i . Function Φ : R≥0 → R|Λ| describes the behaviour of the system as a set of deterministic equations assuming a continuous state-space semantics, and therefore Φ(t) ∈ R|Λ| is the vector of the species concentrations at time t. Assuming t0 = 0, the initial condition is Φ(0) = [x0 ], expressed as a concentration. Note that F (Φ(t)) is Lipschitz continuous, so Φ exists and is unique [26]. 6
Page 6 of 24
ip t
2.2. Stochastic semantics CRNs are well represented by CTMCs, whose transient analysis can be performed via the Chemical Master Equation (CME) [27].
us
cr
Definition 1. Given a CRN C = (Λ, R) and the volumetric factor N , we define a time-homogeneous CTMC [28, 29] (X N (t), t ∈ R≥0 ) with state space S ⊆ N|Λ| . Given x0 ∈ S, the initial configuration of the system, then P (X N (0) = x0 ) = P 1. The transition rate from state xi to state xj is defined as r(xi , xj ) = {τ ∈R|xj =xi +vτ } N αc,τ (xi ).
an
X N (t) describes the stochastic evolution of molecular populations of each species at time t. For x ∈ S, we define P (t) (x) = P (X N (t) = x|X N (0) = x0 ), where x0 is the initial configuration. The CME describes the time evolution of X N as:
M
X d P (t) (x) = {N αc,τ (x − υτ )P (t) (x − υτ ) − N αc,τ (x)P (t) (x)}. dt τ ∈R
(2)
Ac ce p
te
d
The CME can be equivalently defined in terms of the infinitesimal generator matrix [13], which admits computing an approximation of the CME using, for example, fast adaptive uniformisation [30, 19] or the sliding window method [13]. N We also define the CTMC ( X N(t) , t ∈ R≥0 ) with state space S ⊆ Q|Λ| . If N [x0 ]) = 1. The transition [x0 ] ∈ S is the initial configuration, then P ( X N(0) =P rate from state [xi ] to [xj ] is defined as r([xi ], [xj ]) = {τ ∈R|[xj ]=[xi ]+ vτ } N αc,τ (xi ). N
X N (t) N
is the random vector describing the system at time t in terms of conN 0 centrations. In [24, 26] it is proved that lim sup k X N(t ) − Φ(t0 )]k = 0 almost N →∞ t0 ≤t
surely for any time t. This explains the relationship between the two different semantics, where the deterministic solution can be viewed as a limit of the stochastic solution, valid when close enough to the thermodynamic limit. 3. Linear Noise Approximation The solution of the CME can be computationally expensive, or even infeasible, because the set of reachable states can be huge or infinite. The Linear Noise Approximation (LNA) has been introduced by Van Kampen as a second order approximation of the system size expansion of the CME 7
Page 7 of 24
an
us
cr
ip t
[27]. It permits a stochastic characterization of the evolution of a CRN, still maintaining scalability comparable to that of the deterministic models. In what follows, we introduce the LNA following the derivation in [31]. This approach clearly shows that, assuming mass action kinetics, the LNA is always accurate for any CRN, if N is large enough, at least for a limited time. In fact, the original derivation of [27] does not shield much light on the validity of such an approximation because truncating the expansion of the CME works fine only if terms of higher order are well behaved, and this is not always the case [31]. In order to derive the LNA, we first consider the following conditions, namely the leap conditions. Given a CRN C = (Λ, R), X N satisfies the leap conditions at time t, if for any τ ∈ R, it holds that there exists an infinitesimal time interval dt such that:
M
αn,τ (X N (t)) is constant in [t, t + dt] and αn,τ (X N (t)) · dt >>> 1.
te
d
In [32], Gillespie shows that if these conditions are satisfied then the CME can be approximated by the following Chemical Langevin Equation (CLE): X q X √ υτ αn,τ (X N (t))Nτ (0, 1) dt υτ αn,τ (X N (t))dt+ X N (t+dt) = X N (t)+ τ ∈R
τ ∈R
Ac ce p
(3) where Nτ (0, 1) are a set of independent normally distributed random variables with expected value 0 and variance 1. It is possible to show that, assuming mass action kinetics, for N large enough, the leap conditions can always be satisfied, and so Eqn (3) can be considered as a valid approximation of the real Markov process, at least for finite time: the Gaussian nature of the CLE makes it impossible to handle rare events. Since stochastic fluctuations depend on the volumetric factor, N , of the system, and, specifically, 1 for average concentrations are of the order of N 2 [33], we can assume that Eqn (3) has a solution of the form: 1
X N (t) = N Φ(t) + N 2 Z(t)
(4)
where Z(t) = (Z1 (t), Z2 (t), ..., Z|Λ| ) is a random vector, independent of N , representing the stochastic fluctuations at time t and Φ(t) is given by the solution of Eqn (1). This assumption can also be justified by considering the 8
Page 8 of 24
|Λ|
|Λ|
|Λ|
|Λ|
ip t
work of Either and Kurtz [26]. Assuming such a structure for the solution of Eqn (3), then the probability distribution of Z(t) can be approximated by the following linear Fokker-Plank equations [34]:
cr
X X ∂Fj (Φ(t)) ∂(Zj P (Z, t)) 1 X X ∂ 2 P (Z, t) ∂P (Z, t) Gi,j (Φ(t)) =− + ∂t ∂Φi ∂Zi 2 i=1 j=1 ∂Zi ∂Zj i=1 j=1
(5) where G(Φ(t)) = τ ∈R υτ υτ αc,τ (Φ(t)) and Fj (Φ(t)) is the j−th component of F (Φ(t)). The solution of Eqn (5) yields a Gaussian process. For every time t, Z(t) has a multivariate normal distribution whose expected value and covariance matrix are the solution of the following equations [33]: T
us
P
an
dE[Z(t)] = JF (Φ(t))E[Z(t)] dt
(6)
Ac ce p
te
d
M
dC[Z(t)] = JF (Φ(t))C[Z(t)] + C[Z(t)]J T F (Φ(t)) + G(Φ(t)) (7) dt where JF (Φ(t)) is the Jacobian of F (Φ(t)). We consider as initial conditions E[Z(0)] = 0 and C[Z(0)] = 0. This means that E[Z(t)] = 0 for every t. The LNA is therefore obtained as an approximation of the CLE, and so it yields all its conditions of validity. However, we need to check that hypothesis (4) is effectively satisfied, which implies the need to check that C[Z(t)] remains bounded. This ensures that, assuming mass action kinetics, the LNA can be an accurate approximation, at least for a limited time, for any CRN. Even if both LNA and CLE give rise to a Gaussian process, solving the LNA is much simpler than solving the CLE as explained in [31]. As a consequence, it can be used for a fast stochastic characterization of the stochastic semantics of any CRN. LNA is exact in the limit of high populations, but can also be used for small populations if the behaviour is not too far from the deterministic limit, taking into account the continuous nature of the approximation and Gaussian assumptions on the noise. To compute the LNA it is necessary to solve O(|Λ|2 ) first order differential equations, but the complexity is independent of the initial number of molecules of each species. Therefore, one can avoid the exploration of the state space that methods based on uniformization rely upon. 3.1. Probabilistic analysis of CRNs The LNA thus permits approximation of the probability distribution 1 of X N (t) with the probability distribution of Y N (t) = N Φ(t) + N 2 Z(t), 9
Page 9 of 24
an
us
cr
ip t
where Z, and hence Y N , are Gaussian processes. As a consequence, Y N (t) has a multivariate Gaussian distribution, so it is completely characterized by its expected value and covariance matrix, whose values are respectively 1 1 E[Y N (t)] = N Φ(t) and C[Y N (t)] = N 2 C[Z(t)]N 2 = N C[Z(t)]. Since Y N has a multivariate normal distribution, then every linear combination of its components is normally distributed. Therefore, given B = [b1 , b2 , · · · , b|Λ| ] where b1 , b2 , ..., b|Λ| ∈ Z, we can consider the random variable BY N (t), which defines a linear combination of the species at time t. For every t, BY N (t) is a normal random variable, whose expected value and variance are: E[BY N (t)] = BE[Y N (t)] (8) C[BY N (t)] = BC[Y N (t)]B T .
(9)
g(x|E[BY N (tk )], C[BY N (tk )])dx
(10)
d
ΩY N ,B,I (tk ) =
ui X Z
M
For a specific time tk , it is possible to calculate the probability that BY N (tk ) is within a set I of closed, disjoint real intervals [li , ui ], where li , ui ∈ R ∪ {+∞, −∞}. This probability ΩY N ,B,I (tk ) is given by:
[li ,ui ]∈I l i
Ac ce p
te
where g(x|EV, σ 2 ) is the Gaussian distribution with expected value EV and covariance σ 2 . We recall that it is possible to find numerical solution of Eqn (10) in constant time using the Z table [35].
Example 3.1. Consider the CRN of Example 2.1, then we can obtain the probability that #λ1 − 2#λ3 is at least 10 at time 20 by defining B 0 = [1, 0, −2], I 0 = {[10, +∞]} and calculating ΩY N ,B 0 ,I 0 (20). The following theorems are consequences of results in [31], which can be generalized for reactions with a finite number of reagents and products. They show asymptotic pointwise convergence of expected value, variance and probability.
Theorem 1. Let C = (Λ, R) be a CRN. Suppose the solution of Eqn (7) is bounded, then, approaching the thermodynamic limit, for any finite instant of time ti e X N ,B,I (ti )k = 0, lim kΩY N ,B,I (ti ) − Ω (11) N →∞
e X N ,B,I (ti ) is the probability that B(X N ) is within I at time ti . where Ω 10
Page 10 of 24
lim kC[BY N (tk )] − C[BX N (tk )]k = 0
N →∞
cr
lim kE[BY N (tk )] − E[BX N (tk )]k = 0.
ip t
Theorem 2. Suppose the solution of Eqn (7) is bounded, then, approaching the thermodynamic limit, for any finite instant of time tk :
N →∞
(12)
(13)
M
an
us
To solve the differential equations (6) and (7), it is necessary to use a numerical method such as the adaptive Runge-Kutta algorithm [36]. This yields the solution for a finite set of sampling times Σ = [t1 , ..., t|Σ| ] ∈ R|Σ| , where t1 ≤ ... ≤ tk ≤ ... ≤ t|Σ| and |Σ| is the sample size. Assuming Y N is separable, that is, it is possible to completely define the behavior of Y N by only considering a countable number of points, we can calculate ΩY N ,B,I for any point in Σ and, if points are dense enough, then this set exhaustively describes the probability that BX N is within I over time. This restriction is not a limitation since for any stochastic process there exists a separable modification of it [37].
d
4. Stochastic Evolution Logic (SEL)
Ac ce p
te
Let C = (Λ, R) be a CRN with initial state x0 , in a system of size N . We now define the logic SEL (Stochastic Evolution Logic) which enables evaluation of the probability, variance and expectation of linear combinations of populations of the species of C. The syntax of SEL is given by: η := P∼p [B, I][t1 ,t2 ]
| Q∼v [B][t1 ,t2 ]
| η1 ∧ η2
| η1 ∨ η2
where Q = {supV, inf V, supE, inf E}, ∼∈ {<, >}, p ∈ [0, 1], v ∈ R, B ∈ Z|Λ| , I = {[li , ui ] | li , ui ∈ R ∪ [+∞, −∞] ∧ [li , ui ] ∩ [lj , ui ] = ∅, i 6= j} and [t1 , t2 ] is a closed interval, with the constraint that t1 ≤ t2 and t1 , t2 ∈ R. If t1 = t2 the interval reduces to a singleton. Formulae η describe global properties of the stochastic evolution of the system. (B, I) specifies a linear combination of the species of C and a set of intervals, where B ∈ Z|Λ| is the vector defining the linear combination and I represents a set of disjoint closed real intervals. P∼p [B, I][t1 ,t2 ] is the probabilistic operator, which specifies the probability that the linear combination defined by B falls within the range I over the time interval [t1 , t2 ]. 11
Page 11 of 24
ip t
supE, inf E, inf V, supV respectively yield the supremum and infimum of expected value and variance of the random variables associated to B within the specified time interval.
an
us
cr
Example 4.1. Consider the CRN of Example 2.1. Checking if the variance of #λ1 remains smaller than K1 within [tj , tk ] can be expressed as supV0.95 [[1, −1, 0], ([K3 , K4 ], [K2 , ∞])][tj ,tk ] . Equivalently, instead of writing B, we write directly the linear combination it defines. For example, in the latter case we have P>0.95 [(#λ1 −#λ2 ), ([K3 , K4 ], [K2 , ∞])][tj ,tk ] .
d
M
We now comment about expressiveness of SEL in relation to CSL, the logic typically employed to specify properties of the CTMCs induced by CRNs. Though SEL includes the probabilistic operator P∼p [B, I][t1 ,t2 ] , this is different from the probabilistic reachability operator of CSL. As shown in [20], under some restrictions SEL can be endowed with the probabilistic reachability operator, but reward operators for SEL have not been studied. The steady-state operator of CSL cannot be handled by LNA because it is accurate only for finite time.
Ac ce p
te
4.1. Semantics Given a CRN C = (Λ, R) with initial configuration x0 in a system of fixed volumetric factor N , its stochastic behaviour is described by the CTMC X N of Definition 1. We define a path of CTMC X N as a sequence ω = x0 t1 x1 t1 x2 ... where xi is a state and ti ∈ R>0 is the time spent in the state xi . A path is finite if there is a state xk that is absorbing. ω ⊗ t is the state of the path at time t. P ath(X N , x0 ) is the set of all (finite and infinite) paths of the CTMC starting in x0 . We work with the standard probability measure P rob over paths P ath(X N , x0 ) defined using cylinder sets [10]. We first define when a path ω satisfies (B, I) at time t: ω, t |= (B, I)
↔
∃[li , ui ] ∈ I . li ≤ B(ω ⊗ t) ≤ ui .
Note that B(ω ⊗ t) is well defined because ω ⊗ t ∈ N|Λ| . XN We now define P rB,I (t) = P rob{ω ∈ P ath(X N , x0 ) | ω, t |= (B, I)}, then if the time interval is a singleton the satisfaction relation for the probabilistic operator is: X N , x0 |= P∼p [B, I][t1 ,t1 ]
↔
N
X P rB,I (t1 ) ∼ p
12
Page 12 of 24
Instead, for t1 < t2 we have: X , x0 |= P∼p [B, I][t1 ,t2 ]
↔
1 t2 − t1
Z
t2
N
X P rB,I (t) dt ∼ p
t1
N
ip t
N
X N , x0 |= inf V∼v [B][t1 ,t2 ]
↔
X N , x0 |= supE∼v [B][t1 ,t2 ]
↔
X N , x0 |= inf E∼v [B][t1 ,t2 ]
inf (C[B(X N )], [t1 , t2 ]) ∼ v sup(E[B(X N )], [t1 , t2 ]) ∼ v
↔
↔
inf (E[B(X N )], [t1 , t2 ]) ∼ v
X N , x0 |= η1 ∧ X N , x0 |= η2
d
X N , x0 |= η1 ∧ η2
sup(C[B(X N )], [t1 , t2 ]) ∼ v
an
↔
M
X N , x0 |= supV∼v [B][t1 ,t2 ]
us
cr
X P rB,I (t) is the probability of the set of paths of X N such that the linear combination of the species defined by B falls within I. It is well defined since we have previously defined the probability measure P rob on P ath(X N , x0 ). To define the satisfaction relation of the probabilistic operator we simply take XN the average value of P rB,I (t) during the interval [t1 , t2 ]. For the remaining operators the satisfaction relation is defined as:
↔
X N , x0 |= η1 ∨ X N , x0 |= η2
te
X N , x0 |= η1 ∨ η2
Ac ce p
where inf (·, [t1 , t2 ]) and sup(·, [t1 , t2 ]) respectively denote the infimum and supremum within [t1 , t2 ]. 4.2. LNA-based Approximate Model Checking for CRNs Stochastic model checking of CRNs is usually achieved by transient analysis of the CTMC X N [10], which involves solving the CME and thus suffers from the state-space explosion problem. We propose an approximate model checking algorithm based on LNA. The inputs are a SEL formula η, the stochastic process X N induced by the CRN and initial state x0 . The output is true in case the formula is verified, and otherwise f alse. The algorithm proceeds by induction on the structure of formula η, successively computing whether each subformula is satisfied or not. We assume that Eqn (6) and (7) are solved numerically where Σ is the finite set of sample points on which their solution is defined and that t0 , initial time, and tmax , final time, are always sampling points.
13
Page 13 of 24
ip t
4.2.1. Probabilistic operator. To evaluate P∼p [(B, I)][t1 ,t2 ] we construct the function P rob(B,I) (t) = ΩY N ,B,I (ti ) for t ∈ [ti , ti+1 ), ti , ti+1 ∈ Σ (alternatively, it can be constructed as the interpolation of the values of ΩY N ,B,I over Σ points).
cr
Lemma 1. P rob(B,I) is integrable on R≥0 .
us
Proof. P rob(B,I) is a bounded function with at most |MI | discontinuities, where |MI | ∈ N>0 . Therefore, the set of discontinuities is a countable set, and countable sets have measure 0. Hence, the conditions of the Lebesgue criterion for integrability holds. This concludes the proof.
M
an
Theorem 1 guarantees the pointwise correctness of P rob(B,I) and its integrability allows us to compute the following approximation, Rthen compare to t2 1 XN P rB,I (t) dt ≈ threshold p to decide the truth value. If t2 6= t1 then t2 −t t 1 1 R t2 1 XN P robB,I (t)dt else if t1 = t2 then P rB,I (t1 ) ≈ P robB,I (t1 ). t2 −t1 t1
Ac ce p
te
d
4.2.2. Expectation and variance operators To evaluate sup(C[B(X N )], [t1 , t2 ]), inf (C[B(X N )], [t1 , t2 ]), sup(E[B(X N )], [t1 , t2 ]) and inf (E[B(X N )], [t1 , t2 ]) we use the LNA, namely, compute the expected value and variance of Eqn (9) and (8). Theorem 2 guarantees the quality of the approximation. We can now compute the following approximations, then compare to the threshold v: sup(C[B(X N )], [t1 , t2 ]) ≈
max{C[BY N (tk )] | (tk ∈ Σ ∧ t1 ≤ tk ≤ t2 ) ∨ (tk ∈ L[t1 ,t2 ] )}
inf (C[B(X N )], [t1 , t2 ]) ≈
min{C[BY N (tk )] | (tk ∈ Σ ∧ t1 ≤ tk ≤ t2 ) ∨ (tk ∈ L[t1 ,t2 ] )}
and similarly for the expected value. L[t1 ,t2 ] = {ti ∈ Σ | @tj ∈ Σ . |t1 − tj | < |t1 − ti |} ensures that for any time interval there is at least one sampling point, even if the interval is a singleton. Note that, for each sub-formula, the algorithm involves the calculation of some quantity, so one can define a quantitative semantics for SEL as in [38]. In our implementation, the syntax for obtaining this numerical value is by replacing the bounds in the Q and P operators with “=?” as shown in the case studies in Section 5. 14
Page 14 of 24
us
cr
ip t
LNA-based model checking can also be used for systems far from the thermodynamic limit, at a cost of some loss of precision. LNA assumes continuous state space, and it is not possible to justify this assumption for very small populations. However, if the distributions of interest are not multi-modal and the noise term is finite and approximated by a Gaussian distribution, then LNA gives very good approximation even for quite small systems. It is clear that model checking accuracy increases as N grows. We emphasise that the model checking algorithm we have presented is also able to handle CRNs whose stochastic semantics is an infinite CTMC, which occur frequently in biological models.
te
d
M
an
4.2.3. Complexity of LNA-based approximate model checking The time complexity for model checking formula η against a CRN C = (Λ, R) is linear in |η|. In the worst case, analysis of a single operator requires the solution of O(|Λ|2 ) polynomial differential equations for a bounded time. However, an efficient implementation can solve the O(|Λ|2 ) ODEs only once for the interval [0, tmax ], and then reuse this result for every operator, where tmax is the greatest (finite) time of interest. Note that ODEs are solved in terms of concentrations (a value between 0 and 1 by convention), ensuring independence from the number of molecules of each species, although stiffness can slow down the solution of the LNA.
Ac ce p
5. Experimental Results
We implemented the methods in a framework based on Matlab and Java. The experiments were run on an Intel Dual Core i7 machine with 8 GB of RAM. To solve the differential equations, we use M atlab ode45, a variable step Runge-Kutta algorithm. We employ LNA-based model checking for the analysis of four biological reaction networks: a Phosphorelay Network [39], a Gene Expression Model [40, 41], the FGF pathway [42] and the GW network [43]. For every network, the CRN and parameters have been taken from the referenced papers. We coded the same CRNs in PRISM [8] in order to compare accuracy and time of execution with standard uniformisation of the CME [10] and statistical model checking (SMC) techniques (confidence interval method) as implemented in PRISM. For the FGF and GW case studies, global analysis and SMC cannot be used, because the state space is too large for direct analysis, and SMC requires many time-consuming simulations to obtain good accuracy. 15
Page 15 of 24
us
cr
ip t
5.1. Phosphorelay Network We consider a three-layer phosphorelay network whose structure is derived from [39]. Phosphorelay networks are extended two-component signalling systems found in diverse bacteria, lower eukaryotes and plants. Each layer of the network, (L1, L2, L3), can be found in phosphorylate form (L1p, L2p, L3p). We consider the initial condition #L1p = #L2p = #L3p = 0, #L1 = #L2p = #L3p =Init, where Init∈ N. Then we analyse the ligand B, whose initial condition is #B = 3 ∗ Init. We are interested in checking the following SEL property: P>0.7 [(#L1p − #L3p), [0, +∞]][0,100] ∧ P>0.98 [(#L3p − #L1p), [0, +∞]][300,600]
Ac ce p
te
d
M
an
which is verified if, in the first interval, the probability that #L1p is greater than #L3p is > 0.7 and if, between 300 and 600, with probability > 0.98, #L3p is greater than #L1p. We evaluate this formula in three different initial conditions, firstly Init = 32 and N = 5000, then Init= 64 and N = 10000, and finally Init= 100 and N = 15625, so the same concentration but different numbers of molecules. In all cases, the LNA-based model checking evaluates the formula as true. To understand the quality of the approximation, we check the following quantitative formula P=? [(#L3p − #L1p, [0, +∞])][T,T ] for T ∈ [0, 600] (recall that in our implementation =? gives the quantity calculated by model checking the operator). We compare the results with the evaluation of the corresponding CSL formula using standard uniformisation (U nif ) with error 10−7 [10]. The following table shows the results. MaxErr is the maximum error computed by LNA-based approach compared to standard uniformisation and AvgErr is the average error; Time(·) stands for execution time. Init 20 32 64 100
Time (LNA) 0.22 sec 0.23 sec 0.26 sec 0.3 sec
Time (Unif ) 2 min 5 min > 2 hr > 2 hr
MaxErr 0.0675 0.059 0.0448 0.03
AvgErr 0.0519 0.02 0.0027 0.0011
Note that as Init increases the error of our method decreases, while the execution time is practically independent of the molecular count. LNA-based algorithms are faster in all cases. Thus our approach can be used even for quite small population systems, giving fast approximate stochastic characterization. 16
Page 16 of 24
Time (Simul) 75 sec 198 sec 337 sec 483 sec
ExpVal (LNA) 100.17 142.15 159.73 167.1
M
Time (LNA) 0.52 sec 0.54 sec 0.54 sec 0.56 sec
ExpVal (Simul) 100.14 ± 0.1 142.11 ± 0.1 159.74 ± 0.1 167.1 ± 0.1
d
T 300 600 900 1200
an
us
cr
ip t
5.2. Gene Expression We consider a simple CRN that models the transcription of a gene into an mRNA molecule, and the translation of the latter into a protein. The CRN, rates and initial conditions are the same as in [41]. The stochastic semantics of the reaction network is an infinite CTMC, and we use this model to show that our method can handle infinite state-space processes. We consider the quantitative property supE=? [#mRN A][T,T ] , which gives the number of molecules of mRN A in the system at time T . We compare our method with SMC estimation of the same property by using 50000 simulations, for T = {300, 600, 900, 1200}, and in the following tables we compare the results in terms of execution time (T ime(·)) and estimated expected value of #mRN A (ExpV al(·)). LNA-based model checking is several orders of magnitude faster without loss of accuracy.
Ac ce p
te
5.3. Fibroblast Growth Factor (FGF) pathway Fibroblast Growth Factors (FGF) are a family of proteins which play a key role in the process of cell signalling in a variety of contexts, like wound healing and skeletal development. We consider the model of FGF signalling pathway developed in [42], which is composed of more than 50 reactions and species. We consider the system with initially 105 molecules for species with non-zero initial concentration. Analysis of the model reveals that the phosphorylated form of F RS2 can bind the protein Src, and then this new complex, Src:F RS2, can relocate out. We want to check if the expected value of #Src:F RS2 during the first 3000 seconds reaches a maximum value greater than 40. We do that by checking the property supE>40 [#Src:F RS2][0,3000] . The formula evaluates to true, and in Figure 2 we analyze the expected value and standard deviation of #Src:F RS2. We obtain these values directly from the logic considering the quantitative interpretation of supE=? [#Src:F RS2][T,T ] and supV=? [#Src:F RS2][T,T ] for T ∈ [0, 3000]. It is possible to see that, after an initial peak, relocation causes exponential decay.
17
Page 17 of 24
ip t cr us an
M
Figure 2: Expected number and standard deviation of species of #Src:F RS2 in the FGF pathway during the first 8000 seconds estimated by our method is compared with a stochastic simulation of the same species.
Ac ce p
te
d
In the same figure we show a single stochastic simulation of the system for the same initial conditions, confirming our evaluation. Moreover, the approximation can be justified theoretically. #Src:F RS2 converges to zero necessarily and this demonstrates the unimodality of the distribution of the species; we note that the variance is finite, so Eqn (4) holds. 5.4. DNA strand displacement of the GW network GW is a network related to the G2-M cell cycle switch [44]. Under particular initial conditions, it has been shown that GW can emulate the Approximate Majority algorithm [43]. Here, we consider the two-domain DNA strand-displacement implementation of GW [3]. The corresponding CRN is composed of 340 species and 240 reactions. For our analysis the species of interest are R and P , whose initial conditions are #R = 90 and #P = 10. These species model the switch for the activating phosphatase Cdc25; initial conditions of other species are taken from the referenced papers. We check the property P>0.6 [#R − #P, [65, +∞]][0,T ] for T = {1000, 2000, 3500, 5000}, in a system of size N = 45000. The results are reported in the following table.
18
Page 18 of 24
Quantitative Value 0.4297 0.5313 0.6535 0.7349
Qualitative Value false false true true
ip t
Time Execution 420 sec 780 sec 1380 sec 2120 sec
cr
T 1000 2000 3500 5000
6. Concluding Remarks
Ac ce p
te
d
M
an
us
We presented a novel probabilistic logic (SEL) for analysing stochastic behaviour of CRNs and proposed an approximate model checking algorithm of the CME based on the LNA. We have implemented the algorithm and demonstrated on four non-trivial examples that LNA-based model checking enables analysis of CRNs with hundreds of species, and even infinite CTMCs, at a cost of some loss of accuracy. It would be interesting to find bounds on the approximation error when the system is far from the thermodynamic limit. However, the error is not only dependent on the value of N , but also on the structure of the CRN, the rates, and the property. For some systems the leap conditions, which guarantee the quality of the approximation if using the LNA, could be satisfied only for a subset of reactions. As recently shown in [45], it is possible to formulate a stochastic hybrid approach, in which the LNA is used only for a subset of species, namely, those for which the leap conditions are satisfied. Other species are treated as a discrete-state Markov process. This improves precision of the stochastic analysis of CRNs when multimodality is present. It would also be interesting to extend SEL with reward operators as in, e.g., [46, 10], with which one can express properties such as the expected number of molecules and variance of a species when certain events happen. We also intend to release a software tool based on LBS [47]. Acknowledgements
The authors would like to thank Luca Bortolussi for helpful discussions.This research is supported by a Royal Society Research Professorship and ERC AdG VERIWARE.
19
Page 19 of 24
References
ip t
[1] V. Chellaboina, S. Bhat, M. Haddad, D. S. Bernstein, Modeling and analysis of mass-action kinetics, Control Systems, IEEE 29 (4) (2009) 60–78.
us
cr
[2] D. Soloveichik, G. Seelig, E. Winfree, DNA as a universal substrate for chemical kinetics, Proceedings of the National Academy of Sciences 107 (12) (2010) 5393–5398.
[3] L. Cardelli, Two-domain DNA strand displacement, Mathematical Structures in Computer Science 23 (02) (2013) 247–271.
M
an
[4] Y.-J. Chen, N. Dalchau, N. Srinivas, A. Phillips, L. Cardelli, D. Soloveichik, G. Seelig, Programmable chemical controllers made from DNA, Nature Nanotechnology 8 (10) (2013) 755–762.
d
[5] D. T. Gillespie, A. Hellander, L. R. Petzold, Perspective: Stochastic algorithms for chemical kinetics, The Journal of Chemical Physics 138 (17) (2013) 170901.
te
[6] D. T. Gillespie, Deterministic limit of stochastic chemical kinetics, The Journal of Physical Chemistry B 113 (6) (2009) 1640–1644.
Ac ce p
[7] L. Cardelli, On process rate semantics, Theoretical Computer Science 391 (3) (2008) 190–215. [8] M. Kwiatkowska, G. Norman, D. Parker, Prism 4.0: Verification of probabilistic real-time systems, in: International Conference on Computer Aided Verification, Springer, 2011, pp. 585–591. [9] M. Kwiatkowska, C. Thachuk, Probabilistic model checking for biology, Software Systems Safety 36 (2014) 165.
[10] M. Kwiatkowska, G. Norman, D. Parker, Stochastic model checking, in: International School on Formal Methods for the Design of Computer, Communication and Software Systems, Springer, 2007, pp. 220–270. [11] J. Elf, M. Ehrenberg, Fast evaluation of fluctuations in biochemical networks with the linear noise approximation, Genome Research 13 (11) (2003) 2475–2484. 20
Page 20 of 24
ip t
[12] L. Bortolussi, R. Lanciani, Model checking Markov population models by central limit approximation, in: International Conference on Quantitative Evaluation of Systems, Springer, 2013, pp. 123–138.
cr
[13] V. Wolf, R. Goel, M. Mateescu, T. A. Henzinger, Solving the chemical master equation using sliding windows, BMC Systems Biology 4 (1) (2010) 42.
us
[14] B. Munsky, M. Khammash, The finite state projection algorithm for the solution of the chemical master equation, The Journal of chemical physics 124 (4) (2006) 044104.
an
[15] A. Singh, J. P. Hespanha, Lognormal moment closures for biochemical reactions, in: Proceedings of the 45th IEEE Conference on Decision and Control, IEEE, 2006, pp. 2063–2068.
M
[16] J. Hespanha, Moment closure for biochemical networks, in: Communications, Control and Signal Processing, 2008. ISCCSP 2008. 3rd International Symposium on, IEEE, 2008, pp. 142–147.
te
d
[17] A. Aziz, K. Sanwal, V. Singhal, R. Brayton, Model-checking continuoustime Markov chains, ACM Transactions on Computational Logic 1 (1) (2000) 162–170.
Ac ce p
[18] C. Baier, B. Haverkort, H. Hermanns, J.-P. Katoen, Model-checking algorithms for continuous-time Markov chains, IEEE Transactions on Software Engineering 29 (6) (2003) 524–541. [19] F. Dannenberg, E. M. Hahn, M. Kwiatkowska, Computing cumulative rewards using fast adaptive uniformization, ACM Transactions on Modeling and Computer Simulation (TOMACS) 25 (2) (2015) 9. [20] L. Bortolussi, L. Cardelli, M. Kwiatkowska, L. Laurenti, Approximation of probabilistic reachability for chemical reaction networks using the linear noise approximation, in: Proc. 13th International Conference on Quantitative Evaluation of SysTems (QEST 2016), LNCS, Springer, 2016, to appear. [21] P. Ballarini, H. Djafri, M. Duflot, S. Haddad, N. Pekergin, COSMOS: A statistical model checker for the hybrid automata stochastic logic, in:
21
Page 21 of 24
ip t
Proc. 8th Int. Conf. Quantitative Evaluation of SysTems (QEST’11), IEEE Computer Society Press, 2011, pp. 143–144.
cr
[22] P. Ballarini, H. Djafri, M. Duflot, S. Haddad, N. Pekergin, Cosmos: a statistical model checker for the hybrid automata stochastic logic, in: Quantitative Evaluation of Systems (QEST), 2011 Eighth International Conference on, IEEE, 2011, pp. 143–144.
an
us
[23] B. Barbot, M. Kwiatkowska, On quantitative modelling and verification of DNA walker circuits using stochastic petri nets, in: R. Devillers, A. Valmari (Eds.), Application and Theory of Petri Nets and Concurrency, Vol. 9115 of Lecture Notes in Computer Science, Springer International Publishing, 2015, pp. 1–32.
M
[24] D. F. Anderson, T. G. Kurtz, Continuous time Markov chain models for chemical reaction networks, in: Design and Analysis of Biomolecular Circuits, Springer, 2011, pp. 3–42.
d
[25] L. Bortolussi, J. Hillston, D. Latella, M. Massink, Continuous approximation of collective system behaviour: A tutorial, Performance Evaluation 70 (5) (2013) 317–349.
te
[26] S. N. Ethier, T. G. Kurtz, Markov processes: characterization and convergence, Vol. 282, John Wiley & Sons, 2009.
Ac ce p
[27] N. G. Van Kampen, Stochastic processes in physics and chemistry, Vol. 1, Elsevier, 1992. [28] E. Cinlar, Introduction to stochastic processes, Courier Corporation, 2013. [29] M. Pinsky, S. Karlin, An introduction to stochastic modeling, Academic press, 2010. [30] F. Didier, T. A. Henzinger, M. Mateescu, V. Wolf, Fast adaptive uniformization of the chemical master equation, in: High Performance Computational Systems Biology, 2009. HIBI’09. International Workshop on, IEEE, 2009, pp. 118–127. [31] E. Wallace, D. Gillespie, K. Sanft, L. Petzold, Linear noise approximation is valid over limited times for any chemical system that is sufficiently large, IET Systems Biology 6 (4) (2012) 102–115. 22
Page 22 of 24
ip t
[32] D. T. Gillespie, The chemical Langevin equation, The Journal of Chemical Physics 113 (1) (2000) 297–306.
cr
[33] J. Elf, M. Ehrenberg, Fast evaluation of fluctuations in biochemical networks with the linear noise approximation, Genome research 13 (11) (2003) 2475–2484. [34] H. Risken, Fokker-Planck Equation, Springer, 1984.
us
[35] J. K. Patel, C. B. Read, Handbook of the normal distribution, Vol. 150, CRC Press, 1996.
an
[36] J. C. Butcher, The numerical analysis of ordinary differential equations: Runge-Kutta and general linear methods, Wiley-Interscience, 1987.
M
[37] K. It¯o, Essentials of stochastic processes, Vol. 231, American Mathematical Soc., 2006. [38] A. Donz´e, O. Maler, Robust satisfaction of temporal logic over realvalued signals, Springer, 2010.
te
d
[39] A. Csik´asz-Nagy, L. Cardelli, O. S. Soyer, Response dynamics of phosphorelays suggest their potential utility in cell signalling, Journal of The Royal Society Interface 8 (57) (2011) 480–488.
Ac ce p
[40] M. Thattai, A. Van Oudenaarden, Intrinsic noise in gene regulatory networks, Proceedings of the National Academy of Sciences 98 (15) (2001) 8614–8619. [41] M.-E.-C. Mateescu, Propagation models for biochemical reaction networks, Ph.D. thesis, EPFL (2011). [42] J. Heath, M. Kwiatkowska, G. Norman, D. Parker, O. Tymchyshyn, Probabilistic model checking of complex biological pathways, Theoretical Computer Science 391 (3) (2008) 239–257. [43] L. Cardelli, Morphisms of reaction networks that couple structure to function, BMC Systems Biology 8 (1) (2014) 84. [44] B. Novak, J. J. Tyson, Numerical analysis of a comprehensive model of M-phase control in Xenopus oocyte extracts and intact embryos, Journal of Cell Science 106 (4) (1993) 1153–1168. 23
Page 23 of 24
ip t
[45] L. Cardelli, M. Kwiatkowska, L. Laurenti, A stochastic hybrid approximation for chemical kinetics based on the linear noise approximation, in: Proc. 14th International Conference on Computational Methods in Systems Biology (CMSB 2016), LNCS, Springer, 2016, to appear.
us
cr
[46] C. Baier, B. Haverkort, H. Hermanns, J.-P. Katoen, On the logical characterisation of performability properties, in: International Colloquium on Automata, Languages, and Programming, Springer, 2000, pp. 780– 792.
Ac ce p
te
d
M
an
[47] M. Pedersen, G. D. Plotkin, A language for biochemical systems: design and formal specification, in: Transactions on computational systems biology XII, Springer, 2010, pp. 77–145.
24
Page 24 of 24
S0303-2647(16)30203-9 http://dx.doi.org/doi:10.1016/j.biosystems.2016.09.004 BIO 3702
To appear in:
BioSystems
Received date: Revised date: Accepted date:
27-11-2015 8-7-2016 1-9-2016
Please cite this article as: Luca Cardelli, Marta Kwiatkowska, Luca Laurenti, Stochastic Analysis of Chemical Reaction Networks Using Linear Noise Approximation, (2016), http://dx.doi.org/10.1016/j.biosystems.2016.09.004 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
*Manuscript Click here to view linked References
ip t
Stochastic Analysis of Chemical Reaction Networks Using Linear Noise ApproximationI Luca Cardellia,b,∗, Marta Kwiatkowskaa,∗, Luca Laurentia,∗ Department of Computer Science, University of Oxford b Microsoft Research Cambridge
us
cr
a
Abstract
Ac ce p
te
d
M
an
Stochastic evolution of Chemical Reactions Networks (CRNs) over time is usually analysed through solving the Chemical Master Equation (CME) or performing extensive simulations. Analysing stochasticity is often needed, particularly when some molecules occur in low numbers. Unfortunately, both approaches become infeasible if the system is complex and/or it cannot be ensured that initial populations are small. We develop a probabilistic logic for CRNs that enables stochastic analysis of the evolution of populations of molecular species. We present an approximate model checking algorithm based on the Linear Noise Approximation (LNA) of the CME, whose computational complexity is independent of the population size of each species and polynomial in the number of different species. The algorithm requires the solution of first order polynomial differential equations. We prove that our approach is valid for any CRN close enough to the thermodynamical limit. However, we show on four case studies that it can still provide good approximation even for low molecule counts. Our approach enables rigorous analysis of CRNs that are not analyzable by solving the CME, but are far from the deterministic limit. Moreover, it can be used for a fast approximate stochastic characterization of a CRN. Keywords: Chemical Reaction Networks, Linear Noise Approximation, Probabilistic
I
Authors order is alphabetical. Corresponding author Email addresses: [email protected] (Luca Cardelli), [email protected] (Marta Kwiatkowska), [email protected] (Luca Laurenti) ∗
Preprint submitted to BioSystems
July 8, 2016
Page 1 of 24
ip t
Logic, Model Checking 1. Introduction
Ac ce p
te
d
M
an
us
cr
Chemical reaction networks (CRNs) and mass action kinetics are well studied formalisms for modelling biochemical systems [1]. In recent years, CRNs have also been successfully used as a formal programming language for biochemical systems [2, 3, 4]. There are two well established approaches for analyzing chemical networks: deterministic and stochastic [5]. The deterministic approach models the kinetics of a CRN as a system of ordinary differential equations (ODEs) and represents average behaviour, valid in the thermodynamic limit, when the molecular population is sufficiently high [6]. The stochastic approach, on the other hand, is based on the Chemical Master Equation (CME) and models the CRN as a continuous-time Markov chain (CTMC) [7]. The stochastic behavior can be analyzed by stochastic simulation [5] or by exhaustive probabilistic model checking of the CTMC, which can be performed, for example, by using PRISM [8]. Exhaustive analysis of the CTMC is able to find the best- and worst-case scenarios and is correct for any population size, but suffers from the state-space explosion problem [9] and can only be used for relatively small systems. In contrast, deterministic methods are much more robust with respect to state-space explosion, but unable to represent stochastic fluctuations, which play a fundamental role when the system is not in thermodynamic equilibrium. 1.0.1. Contributions. In this paper we develop a novel approach for analysing the stochastic evolution of a CRN based on the Linear Noise Approximation (LNA) of the CME. We formulate SEL (Stochastic Evolution Logic), a probabilistic logic for CRNs that enables reasoning about probability, expectation and variance of linear combinations of populations of the species. Examples of properties that can be specified in our logic include “the maximum expected population size of species λ1 during the first 20 seconds is 75 molecules” and “the probability that the combined population of species λ1 and λ2 has degraded between 10 and 30 seconds is less than 0.1”. We propose an approximate model checking algorithm for the logic based on the LNA and implement it in Matlab and Java. We demonstrate that the complexity of model checking is polynomial in the initial number of species and independent of the initial molecule counts, thus ameliorating state-space explosion. Further, we 2
Page 2 of 24
an
us
cr
ip t
show that model checking is exact when approaching the thermodynamic limit. Though the algorithm may not be accurate for systems far from the deterministic limit, this generally happens when the populations are small, in which case the analysis can be performed by transient analysis of the induced CTMC [10]. Our approach is essential for CRNs that cannot be analyzed by (partial) state space exploration, because of large or infinite state spaces. Moreover, it is useful for a fast (approximate) stochastic characterization of CRNs, since solving the LNA is much faster than solving the CME [11]. We prove asymptotic correctness of LNA-based model checking and show on four examples that it is still possible to obtain very good approximations even for small population systems, compared to standard uniformisation [10] and statistical model checking implemented in PRISM [8].
Ac ce p
te
d
M
1.0.2. Related work The closest work to ours is by Bortolussi et al. [12], which uses the Central Limit Approximation (CLA) (essentially the same as the LNA) for checking restricted timed automata specifications, assuming a fixed population size. Wolf et al. [13] develop a sliding window method to approximately verify infinite-state CTMCs, which applies to cases where most of the probability mass is concentrated in a confined region of the state space. Recently, Finite State Projection algorithms (FSP algorithms) for the solution or approximation of the CME have been introduced [14]. Sliding window and FSP algorithms apply to the induced CTMC, but require at least partial exploration of the state space, and are thus not immune to state-space explosion. Moment closure techniques [15, 16] improve scalability by estimating the first k ∈ N moments of the distribution of the species over time. The LNA itself can be seen as a moment closure technique, as a Gaussian distribution is completely characterized by the first two moments. However, the LNA tells us more because it guarantees that, if certain conditions are satisfied, the distribution of the process is Gaussian. Continuous Stochastic Logic (CSL), originally introduced in [17] and extended by Baier et al. [18], is a logic widely used to perform model checking of continuous-time Markov chains. CSL combines temporal operators of the logic CTL with the probabilistic and steady-state operators, and is further extended with reward operators in [10]. CSL model checking is based on solving the CME and proceeds through uniformisation of the CTMC, essentially a time discretisation, and thus involves traversal of the full state space. This can be partially ameliorated by fast adaptive uniformisation [19] that does 3
Page 3 of 24
te
d
M
an
us
cr
ip t
not consider states with negligible probability. An alternative is statistical model checking (SMC) which involves A key operator of CSL is probabilistic reachability, that is, computing the probability that a particular region of the state space is reached over a given time interval. Although SEL is endowed with a probabilistic operator, this operator gives the average value of the probability over time and, if the time interval is not a singleton, this is not equivalent to probabilistic reachability. Nevertheless, as shown in [20], SEL and our approximate model checking algorithm can be extended to express reachability, but currently lacks reward operators. The CSL steady-state operator of CSL cannot be added to CSL because LNA is acccurate only for finite time. PRISM [8] implements CSL model checking using uniformisation, fast adaptive uniformisation and statistical model checking. Hybrid Automata Stochastic Logic HASL [21] is an expressive specification formalism for stochastic Petri nets based on linear hybrid automata that is employed by the tool Cosmos [22]. CRNs have a natural interpretation in terms of stochastic Petri nets, see e.g. [23]. The HASL formalism is more expressive than SEL and CSL, and can express CSL probabilistic reachability and expected reward properties. HASL model checking proceeds through statistical model checking of the product of a HASL specification automaton and the Petri net, and is implemented in Cosmos. In contrast, SEL model checking follows through approximating the solution of the CME with the Gaussian process induced by the LNA.
Ac ce p
1.0.3. Structure of the paper In Section 2 we summarise the deterministic and stochastic modelling approaches for CRNs, and in Section 3 we describe the Linear Noise Approximation method. Section 4 introduces the logic SEL and the corresponding model checking algorithm based on the LNA. In Section 5 we demonstrate our approach on four networks taken from the literature. Section 6 concludes the paper. 2. Chemical Reaction Networks
A chemical reaction network (CRN) C = (Λ, R) is a pair of finite sets, where Λ is the set of chemical species and R the set of reactions. |Λ| denotes the size of the set of species. A reaction τ ∈ R is a triple τ = (rτ , pτ , kτ ), where rτ , pτ ∈ N|Λ| and kτ ∈ R>0 . rτ and pτ represent the stoichiometry of reactants and products and kτ is the coefficient associated to the rate of the 4
Page 4 of 24
Ac ce p
te
d
M
an
us
cr
ip t
reaction; its dimension is s−1 . We often write reactions as λ1 + λ3 →k1 2λ2 instead of τ1 = ([1, 0, 1]T , [0, 2, 0]T , k1 ), where ·T indicates the transpose of a vector. We define the net change associated to a reaction τ by υτ = pτ − rτ . For example, for τ1 as above, we have υτ1 = [−1, 2, −1]T . We make the assumption that the system is well stirred, that is, the probability of the next reaction occurring between two molecules is independent of the location of those molecules. We consider fixed volume V and temperature; under these assumptions a configuration or state x ∈ N|Λ| of the system is given by the number of molecules of each species. We define [x] = Nx , the vector of the species concentration in x for a given N , where N = V · NA is the volumetric factor, V is the volume of the solution and NA is Avogadro’s number. The physical dimension of N is M ol−1 · L, where M ol indicates mole and L is litre. Given λi ∈ Λ then #λi x ∈ N represents the number of molecules of λi in x and [λi ] x ∈ R the concentration of λi in the same configuration. In some cases we elide x, and we simply write #λi and [λi ] i . The dimension instead of #λi x and [λi ] x. They are related by [λi ] = #λ N −1 of [λi ] is M ol · L . The propensity αn,τ of a reaction τ in terms of the number of molecules (here subscript n stands for the numebr of molecules) is a function of the current configuration of the system x such that αn,τ (x)dt is the probability that a reaction event occurs in the next infinitesimal interval dt. In this paper we assume as valid the stochastic form of the law of mass action, so the propensity rates are proportional to the number of molecules that participate in the reaction [7]. Stochastic models consider the system in terms of numbers of molecules, while deterministic ones, generally, in terms of concentrations, denoted αc,τ (x) where subscript c stands for concentrations, and the relationship is as follows. For a reaction τ = (rτ , pτ , kτ ), given the configuQ rτ,i ration x and rτ,i , the i-th component of rτ , then αc,τ (x) = kτ |Λ| i=1 ([λi ] x) is the propensity function expressed in terms of concentrations as given by the deterministic law of mass action. It is possible to show that, for any order of reaction, αn,τ (x) ≈ N αc,τ (x) if N is sufficiently large [24]. Note that αc,τ is independent of N . In this paper we are interested only in finite time horizon, because of the problematic character of studying solutions of ODEs for infinite time horizon [25]. Example 2.1. Consider the CRN C = ({λ1 , λ2 , λ3 }, R), where R = {(λ1 + λ2 →10 λ2 + λ2 ), (λ2 + λ3 →10 λ3 + λ3 )}, with initial conditions #λ1 = 98, #λ2 = 1, #λ3 = 1, for a system with N = 1000. Figure 1 plots the ex5
Page 5 of 24
ip t cr us
an
Figure 1: Expected number and standard deviation of species of the CRN of Example 2.1 for the given initial conditions, calculated by simulating the CME.
te
d
M
pectation and standard deviation of population sizes. We may wish to check if the maximum expected value of #λ2 remains smaller than 75 molecules during the first 2sec. However, the system is stochastic, so we also need to analyse whether the variance is limited enough when #λ2 reaches the maximum. Sometimes, analysis of first and second moments does not suffice, so it could be of interest to check the probability of some events, for instance, is the probability that #λ2 − (#λ1 + #λ3 ) > 0, between t1 = 0.5sec and t2 = 1.0sec, greater than 0.6?
Ac ce p
2.1. Deterministic semantics Let C = (Λ, R) be a CRN. The deterministic model approximates the concentration of the species of the system over time as a set of autonomous polynomial first order differential equations: dΦ(t) = F (Φ(t)) dt
(1)
P Q|Λ| where F (Φ(t)) = τ =(rτ ,pτ ,kτ )∈R υτ αc,τ (Φ(t)) and αc,τ (Φ(t)) = kτ i=1 Φi (t)rτ,i . Function Φ : R≥0 → R|Λ| describes the behaviour of the system as a set of deterministic equations assuming a continuous state-space semantics, and therefore Φ(t) ∈ R|Λ| is the vector of the species concentrations at time t. Assuming t0 = 0, the initial condition is Φ(0) = [x0 ], expressed as a concentration. Note that F (Φ(t)) is Lipschitz continuous, so Φ exists and is unique [26]. 6
Page 6 of 24
ip t
2.2. Stochastic semantics CRNs are well represented by CTMCs, whose transient analysis can be performed via the Chemical Master Equation (CME) [27].
us
cr
Definition 1. Given a CRN C = (Λ, R) and the volumetric factor N , we define a time-homogeneous CTMC [28, 29] (X N (t), t ∈ R≥0 ) with state space S ⊆ N|Λ| . Given x0 ∈ S, the initial configuration of the system, then P (X N (0) = x0 ) = P 1. The transition rate from state xi to state xj is defined as r(xi , xj ) = {τ ∈R|xj =xi +vτ } N αc,τ (xi ).
an
X N (t) describes the stochastic evolution of molecular populations of each species at time t. For x ∈ S, we define P (t) (x) = P (X N (t) = x|X N (0) = x0 ), where x0 is the initial configuration. The CME describes the time evolution of X N as:
M
X d P (t) (x) = {N αc,τ (x − υτ )P (t) (x − υτ ) − N αc,τ (x)P (t) (x)}. dt τ ∈R
(2)
Ac ce p
te
d
The CME can be equivalently defined in terms of the infinitesimal generator matrix [13], which admits computing an approximation of the CME using, for example, fast adaptive uniformisation [30, 19] or the sliding window method [13]. N We also define the CTMC ( X N(t) , t ∈ R≥0 ) with state space S ⊆ Q|Λ| . If N [x0 ]) = 1. The transition [x0 ] ∈ S is the initial configuration, then P ( X N(0) =P rate from state [xi ] to [xj ] is defined as r([xi ], [xj ]) = {τ ∈R|[xj ]=[xi ]+ vτ } N αc,τ (xi ). N
X N (t) N
is the random vector describing the system at time t in terms of conN 0 centrations. In [24, 26] it is proved that lim sup k X N(t ) − Φ(t0 )]k = 0 almost N →∞ t0 ≤t
surely for any time t. This explains the relationship between the two different semantics, where the deterministic solution can be viewed as a limit of the stochastic solution, valid when close enough to the thermodynamic limit. 3. Linear Noise Approximation The solution of the CME can be computationally expensive, or even infeasible, because the set of reachable states can be huge or infinite. The Linear Noise Approximation (LNA) has been introduced by Van Kampen as a second order approximation of the system size expansion of the CME 7
Page 7 of 24
an
us
cr
ip t
[27]. It permits a stochastic characterization of the evolution of a CRN, still maintaining scalability comparable to that of the deterministic models. In what follows, we introduce the LNA following the derivation in [31]. This approach clearly shows that, assuming mass action kinetics, the LNA is always accurate for any CRN, if N is large enough, at least for a limited time. In fact, the original derivation of [27] does not shield much light on the validity of such an approximation because truncating the expansion of the CME works fine only if terms of higher order are well behaved, and this is not always the case [31]. In order to derive the LNA, we first consider the following conditions, namely the leap conditions. Given a CRN C = (Λ, R), X N satisfies the leap conditions at time t, if for any τ ∈ R, it holds that there exists an infinitesimal time interval dt such that:
M
αn,τ (X N (t)) is constant in [t, t + dt] and αn,τ (X N (t)) · dt >>> 1.
te
d
In [32], Gillespie shows that if these conditions are satisfied then the CME can be approximated by the following Chemical Langevin Equation (CLE): X q X √ υτ αn,τ (X N (t))Nτ (0, 1) dt υτ αn,τ (X N (t))dt+ X N (t+dt) = X N (t)+ τ ∈R
τ ∈R
Ac ce p
(3) where Nτ (0, 1) are a set of independent normally distributed random variables with expected value 0 and variance 1. It is possible to show that, assuming mass action kinetics, for N large enough, the leap conditions can always be satisfied, and so Eqn (3) can be considered as a valid approximation of the real Markov process, at least for finite time: the Gaussian nature of the CLE makes it impossible to handle rare events. Since stochastic fluctuations depend on the volumetric factor, N , of the system, and, specifically, 1 for average concentrations are of the order of N 2 [33], we can assume that Eqn (3) has a solution of the form: 1
X N (t) = N Φ(t) + N 2 Z(t)
(4)
where Z(t) = (Z1 (t), Z2 (t), ..., Z|Λ| ) is a random vector, independent of N , representing the stochastic fluctuations at time t and Φ(t) is given by the solution of Eqn (1). This assumption can also be justified by considering the 8
Page 8 of 24
|Λ|
|Λ|
|Λ|
|Λ|
ip t
work of Either and Kurtz [26]. Assuming such a structure for the solution of Eqn (3), then the probability distribution of Z(t) can be approximated by the following linear Fokker-Plank equations [34]:
cr
X X ∂Fj (Φ(t)) ∂(Zj P (Z, t)) 1 X X ∂ 2 P (Z, t) ∂P (Z, t) Gi,j (Φ(t)) =− + ∂t ∂Φi ∂Zi 2 i=1 j=1 ∂Zi ∂Zj i=1 j=1
(5) where G(Φ(t)) = τ ∈R υτ υτ αc,τ (Φ(t)) and Fj (Φ(t)) is the j−th component of F (Φ(t)). The solution of Eqn (5) yields a Gaussian process. For every time t, Z(t) has a multivariate normal distribution whose expected value and covariance matrix are the solution of the following equations [33]: T
us
P
an
dE[Z(t)] = JF (Φ(t))E[Z(t)] dt
(6)
Ac ce p
te
d
M
dC[Z(t)] = JF (Φ(t))C[Z(t)] + C[Z(t)]J T F (Φ(t)) + G(Φ(t)) (7) dt where JF (Φ(t)) is the Jacobian of F (Φ(t)). We consider as initial conditions E[Z(0)] = 0 and C[Z(0)] = 0. This means that E[Z(t)] = 0 for every t. The LNA is therefore obtained as an approximation of the CLE, and so it yields all its conditions of validity. However, we need to check that hypothesis (4) is effectively satisfied, which implies the need to check that C[Z(t)] remains bounded. This ensures that, assuming mass action kinetics, the LNA can be an accurate approximation, at least for a limited time, for any CRN. Even if both LNA and CLE give rise to a Gaussian process, solving the LNA is much simpler than solving the CLE as explained in [31]. As a consequence, it can be used for a fast stochastic characterization of the stochastic semantics of any CRN. LNA is exact in the limit of high populations, but can also be used for small populations if the behaviour is not too far from the deterministic limit, taking into account the continuous nature of the approximation and Gaussian assumptions on the noise. To compute the LNA it is necessary to solve O(|Λ|2 ) first order differential equations, but the complexity is independent of the initial number of molecules of each species. Therefore, one can avoid the exploration of the state space that methods based on uniformization rely upon. 3.1. Probabilistic analysis of CRNs The LNA thus permits approximation of the probability distribution 1 of X N (t) with the probability distribution of Y N (t) = N Φ(t) + N 2 Z(t), 9
Page 9 of 24
an
us
cr
ip t
where Z, and hence Y N , are Gaussian processes. As a consequence, Y N (t) has a multivariate Gaussian distribution, so it is completely characterized by its expected value and covariance matrix, whose values are respectively 1 1 E[Y N (t)] = N Φ(t) and C[Y N (t)] = N 2 C[Z(t)]N 2 = N C[Z(t)]. Since Y N has a multivariate normal distribution, then every linear combination of its components is normally distributed. Therefore, given B = [b1 , b2 , · · · , b|Λ| ] where b1 , b2 , ..., b|Λ| ∈ Z, we can consider the random variable BY N (t), which defines a linear combination of the species at time t. For every t, BY N (t) is a normal random variable, whose expected value and variance are: E[BY N (t)] = BE[Y N (t)] (8) C[BY N (t)] = BC[Y N (t)]B T .
(9)
g(x|E[BY N (tk )], C[BY N (tk )])dx
(10)
d
ΩY N ,B,I (tk ) =
ui X Z
M
For a specific time tk , it is possible to calculate the probability that BY N (tk ) is within a set I of closed, disjoint real intervals [li , ui ], where li , ui ∈ R ∪ {+∞, −∞}. This probability ΩY N ,B,I (tk ) is given by:
[li ,ui ]∈I l i
Ac ce p
te
where g(x|EV, σ 2 ) is the Gaussian distribution with expected value EV and covariance σ 2 . We recall that it is possible to find numerical solution of Eqn (10) in constant time using the Z table [35].
Example 3.1. Consider the CRN of Example 2.1, then we can obtain the probability that #λ1 − 2#λ3 is at least 10 at time 20 by defining B 0 = [1, 0, −2], I 0 = {[10, +∞]} and calculating ΩY N ,B 0 ,I 0 (20). The following theorems are consequences of results in [31], which can be generalized for reactions with a finite number of reagents and products. They show asymptotic pointwise convergence of expected value, variance and probability.
Theorem 1. Let C = (Λ, R) be a CRN. Suppose the solution of Eqn (7) is bounded, then, approaching the thermodynamic limit, for any finite instant of time ti e X N ,B,I (ti )k = 0, lim kΩY N ,B,I (ti ) − Ω (11) N →∞
e X N ,B,I (ti ) is the probability that B(X N ) is within I at time ti . where Ω 10
Page 10 of 24
lim kC[BY N (tk )] − C[BX N (tk )]k = 0
N →∞
cr
lim kE[BY N (tk )] − E[BX N (tk )]k = 0.
ip t
Theorem 2. Suppose the solution of Eqn (7) is bounded, then, approaching the thermodynamic limit, for any finite instant of time tk :
N →∞
(12)
(13)
M
an
us
To solve the differential equations (6) and (7), it is necessary to use a numerical method such as the adaptive Runge-Kutta algorithm [36]. This yields the solution for a finite set of sampling times Σ = [t1 , ..., t|Σ| ] ∈ R|Σ| , where t1 ≤ ... ≤ tk ≤ ... ≤ t|Σ| and |Σ| is the sample size. Assuming Y N is separable, that is, it is possible to completely define the behavior of Y N by only considering a countable number of points, we can calculate ΩY N ,B,I for any point in Σ and, if points are dense enough, then this set exhaustively describes the probability that BX N is within I over time. This restriction is not a limitation since for any stochastic process there exists a separable modification of it [37].
d
4. Stochastic Evolution Logic (SEL)
Ac ce p
te
Let C = (Λ, R) be a CRN with initial state x0 , in a system of size N . We now define the logic SEL (Stochastic Evolution Logic) which enables evaluation of the probability, variance and expectation of linear combinations of populations of the species of C. The syntax of SEL is given by: η := P∼p [B, I][t1 ,t2 ]
| Q∼v [B][t1 ,t2 ]
| η1 ∧ η2
| η1 ∨ η2
where Q = {supV, inf V, supE, inf E}, ∼∈ {<, >}, p ∈ [0, 1], v ∈ R, B ∈ Z|Λ| , I = {[li , ui ] | li , ui ∈ R ∪ [+∞, −∞] ∧ [li , ui ] ∩ [lj , ui ] = ∅, i 6= j} and [t1 , t2 ] is a closed interval, with the constraint that t1 ≤ t2 and t1 , t2 ∈ R. If t1 = t2 the interval reduces to a singleton. Formulae η describe global properties of the stochastic evolution of the system. (B, I) specifies a linear combination of the species of C and a set of intervals, where B ∈ Z|Λ| is the vector defining the linear combination and I represents a set of disjoint closed real intervals. P∼p [B, I][t1 ,t2 ] is the probabilistic operator, which specifies the probability that the linear combination defined by B falls within the range I over the time interval [t1 , t2 ]. 11
Page 11 of 24
ip t
supE, inf E, inf V, supV respectively yield the supremum and infimum of expected value and variance of the random variables associated to B within the specified time interval.
an
us
cr
Example 4.1. Consider the CRN of Example 2.1. Checking if the variance of #λ1 remains smaller than K1 within [tj , tk ] can be expressed as supV
d
M
We now comment about expressiveness of SEL in relation to CSL, the logic typically employed to specify properties of the CTMCs induced by CRNs. Though SEL includes the probabilistic operator P∼p [B, I][t1 ,t2 ] , this is different from the probabilistic reachability operator of CSL. As shown in [20], under some restrictions SEL can be endowed with the probabilistic reachability operator, but reward operators for SEL have not been studied. The steady-state operator of CSL cannot be handled by LNA because it is accurate only for finite time.
Ac ce p
te
4.1. Semantics Given a CRN C = (Λ, R) with initial configuration x0 in a system of fixed volumetric factor N , its stochastic behaviour is described by the CTMC X N of Definition 1. We define a path of CTMC X N as a sequence ω = x0 t1 x1 t1 x2 ... where xi is a state and ti ∈ R>0 is the time spent in the state xi . A path is finite if there is a state xk that is absorbing. ω ⊗ t is the state of the path at time t. P ath(X N , x0 ) is the set of all (finite and infinite) paths of the CTMC starting in x0 . We work with the standard probability measure P rob over paths P ath(X N , x0 ) defined using cylinder sets [10]. We first define when a path ω satisfies (B, I) at time t: ω, t |= (B, I)
↔
∃[li , ui ] ∈ I . li ≤ B(ω ⊗ t) ≤ ui .
Note that B(ω ⊗ t) is well defined because ω ⊗ t ∈ N|Λ| . XN We now define P rB,I (t) = P rob{ω ∈ P ath(X N , x0 ) | ω, t |= (B, I)}, then if the time interval is a singleton the satisfaction relation for the probabilistic operator is: X N , x0 |= P∼p [B, I][t1 ,t1 ]
↔
N
X P rB,I (t1 ) ∼ p
12
Page 12 of 24
Instead, for t1 < t2 we have: X , x0 |= P∼p [B, I][t1 ,t2 ]
↔
1 t2 − t1
Z
t2
N
X P rB,I (t) dt ∼ p
t1
N
ip t
N
X N , x0 |= inf V∼v [B][t1 ,t2 ]
↔
X N , x0 |= supE∼v [B][t1 ,t2 ]
↔
X N , x0 |= inf E∼v [B][t1 ,t2 ]
inf (C[B(X N )], [t1 , t2 ]) ∼ v sup(E[B(X N )], [t1 , t2 ]) ∼ v
↔
↔
inf (E[B(X N )], [t1 , t2 ]) ∼ v
X N , x0 |= η1 ∧ X N , x0 |= η2
d
X N , x0 |= η1 ∧ η2
sup(C[B(X N )], [t1 , t2 ]) ∼ v
an
↔
M
X N , x0 |= supV∼v [B][t1 ,t2 ]
us
cr
X P rB,I (t) is the probability of the set of paths of X N such that the linear combination of the species defined by B falls within I. It is well defined since we have previously defined the probability measure P rob on P ath(X N , x0 ). To define the satisfaction relation of the probabilistic operator we simply take XN the average value of P rB,I (t) during the interval [t1 , t2 ]. For the remaining operators the satisfaction relation is defined as:
↔
X N , x0 |= η1 ∨ X N , x0 |= η2
te
X N , x0 |= η1 ∨ η2
Ac ce p
where inf (·, [t1 , t2 ]) and sup(·, [t1 , t2 ]) respectively denote the infimum and supremum within [t1 , t2 ]. 4.2. LNA-based Approximate Model Checking for CRNs Stochastic model checking of CRNs is usually achieved by transient analysis of the CTMC X N [10], which involves solving the CME and thus suffers from the state-space explosion problem. We propose an approximate model checking algorithm based on LNA. The inputs are a SEL formula η, the stochastic process X N induced by the CRN and initial state x0 . The output is true in case the formula is verified, and otherwise f alse. The algorithm proceeds by induction on the structure of formula η, successively computing whether each subformula is satisfied or not. We assume that Eqn (6) and (7) are solved numerically where Σ is the finite set of sample points on which their solution is defined and that t0 , initial time, and tmax , final time, are always sampling points.
13
Page 13 of 24
ip t
4.2.1. Probabilistic operator. To evaluate P∼p [(B, I)][t1 ,t2 ] we construct the function P rob(B,I) (t) = ΩY N ,B,I (ti ) for t ∈ [ti , ti+1 ), ti , ti+1 ∈ Σ (alternatively, it can be constructed as the interpolation of the values of ΩY N ,B,I over Σ points).
cr
Lemma 1. P rob(B,I) is integrable on R≥0 .
us
Proof. P rob(B,I) is a bounded function with at most |MI | discontinuities, where |MI | ∈ N>0 . Therefore, the set of discontinuities is a countable set, and countable sets have measure 0. Hence, the conditions of the Lebesgue criterion for integrability holds. This concludes the proof.
M
an
Theorem 1 guarantees the pointwise correctness of P rob(B,I) and its integrability allows us to compute the following approximation, Rthen compare to t2 1 XN P rB,I (t) dt ≈ threshold p to decide the truth value. If t2 6= t1 then t2 −t t 1 1 R t2 1 XN P robB,I (t)dt else if t1 = t2 then P rB,I (t1 ) ≈ P robB,I (t1 ). t2 −t1 t1
Ac ce p
te
d
4.2.2. Expectation and variance operators To evaluate sup(C[B(X N )], [t1 , t2 ]), inf (C[B(X N )], [t1 , t2 ]), sup(E[B(X N )], [t1 , t2 ]) and inf (E[B(X N )], [t1 , t2 ]) we use the LNA, namely, compute the expected value and variance of Eqn (9) and (8). Theorem 2 guarantees the quality of the approximation. We can now compute the following approximations, then compare to the threshold v: sup(C[B(X N )], [t1 , t2 ]) ≈
max{C[BY N (tk )] | (tk ∈ Σ ∧ t1 ≤ tk ≤ t2 ) ∨ (tk ∈ L[t1 ,t2 ] )}
inf (C[B(X N )], [t1 , t2 ]) ≈
min{C[BY N (tk )] | (tk ∈ Σ ∧ t1 ≤ tk ≤ t2 ) ∨ (tk ∈ L[t1 ,t2 ] )}
and similarly for the expected value. L[t1 ,t2 ] = {ti ∈ Σ | @tj ∈ Σ . |t1 − tj | < |t1 − ti |} ensures that for any time interval there is at least one sampling point, even if the interval is a singleton. Note that, for each sub-formula, the algorithm involves the calculation of some quantity, so one can define a quantitative semantics for SEL as in [38]. In our implementation, the syntax for obtaining this numerical value is by replacing the bounds in the Q and P operators with “=?” as shown in the case studies in Section 5. 14
Page 14 of 24
us
cr
ip t
LNA-based model checking can also be used for systems far from the thermodynamic limit, at a cost of some loss of precision. LNA assumes continuous state space, and it is not possible to justify this assumption for very small populations. However, if the distributions of interest are not multi-modal and the noise term is finite and approximated by a Gaussian distribution, then LNA gives very good approximation even for quite small systems. It is clear that model checking accuracy increases as N grows. We emphasise that the model checking algorithm we have presented is also able to handle CRNs whose stochastic semantics is an infinite CTMC, which occur frequently in biological models.
te
d
M
an
4.2.3. Complexity of LNA-based approximate model checking The time complexity for model checking formula η against a CRN C = (Λ, R) is linear in |η|. In the worst case, analysis of a single operator requires the solution of O(|Λ|2 ) polynomial differential equations for a bounded time. However, an efficient implementation can solve the O(|Λ|2 ) ODEs only once for the interval [0, tmax ], and then reuse this result for every operator, where tmax is the greatest (finite) time of interest. Note that ODEs are solved in terms of concentrations (a value between 0 and 1 by convention), ensuring independence from the number of molecules of each species, although stiffness can slow down the solution of the LNA.
Ac ce p
5. Experimental Results
We implemented the methods in a framework based on Matlab and Java. The experiments were run on an Intel Dual Core i7 machine with 8 GB of RAM. To solve the differential equations, we use M atlab ode45, a variable step Runge-Kutta algorithm. We employ LNA-based model checking for the analysis of four biological reaction networks: a Phosphorelay Network [39], a Gene Expression Model [40, 41], the FGF pathway [42] and the GW network [43]. For every network, the CRN and parameters have been taken from the referenced papers. We coded the same CRNs in PRISM [8] in order to compare accuracy and time of execution with standard uniformisation of the CME [10] and statistical model checking (SMC) techniques (confidence interval method) as implemented in PRISM. For the FGF and GW case studies, global analysis and SMC cannot be used, because the state space is too large for direct analysis, and SMC requires many time-consuming simulations to obtain good accuracy. 15
Page 15 of 24
us
cr
ip t
5.1. Phosphorelay Network We consider a three-layer phosphorelay network whose structure is derived from [39]. Phosphorelay networks are extended two-component signalling systems found in diverse bacteria, lower eukaryotes and plants. Each layer of the network, (L1, L2, L3), can be found in phosphorylate form (L1p, L2p, L3p). We consider the initial condition #L1p = #L2p = #L3p = 0, #L1 = #L2p = #L3p =Init, where Init∈ N. Then we analyse the ligand B, whose initial condition is #B = 3 ∗ Init. We are interested in checking the following SEL property: P>0.7 [(#L1p − #L3p), [0, +∞]][0,100] ∧ P>0.98 [(#L3p − #L1p), [0, +∞]][300,600]
Ac ce p
te
d
M
an
which is verified if, in the first interval, the probability that #L1p is greater than #L3p is > 0.7 and if, between 300 and 600, with probability > 0.98, #L3p is greater than #L1p. We evaluate this formula in three different initial conditions, firstly Init = 32 and N = 5000, then Init= 64 and N = 10000, and finally Init= 100 and N = 15625, so the same concentration but different numbers of molecules. In all cases, the LNA-based model checking evaluates the formula as true. To understand the quality of the approximation, we check the following quantitative formula P=? [(#L3p − #L1p, [0, +∞])][T,T ] for T ∈ [0, 600] (recall that in our implementation =? gives the quantity calculated by model checking the operator). We compare the results with the evaluation of the corresponding CSL formula using standard uniformisation (U nif ) with error 10−7 [10]. The following table shows the results. MaxErr is the maximum error computed by LNA-based approach compared to standard uniformisation and AvgErr is the average error; Time(·) stands for execution time. Init 20 32 64 100
Time (LNA) 0.22 sec 0.23 sec 0.26 sec 0.3 sec
Time (Unif ) 2 min 5 min > 2 hr > 2 hr
MaxErr 0.0675 0.059 0.0448 0.03
AvgErr 0.0519 0.02 0.0027 0.0011
Note that as Init increases the error of our method decreases, while the execution time is practically independent of the molecular count. LNA-based algorithms are faster in all cases. Thus our approach can be used even for quite small population systems, giving fast approximate stochastic characterization. 16
Page 16 of 24
Time (Simul) 75 sec 198 sec 337 sec 483 sec
ExpVal (LNA) 100.17 142.15 159.73 167.1
M
Time (LNA) 0.52 sec 0.54 sec 0.54 sec 0.56 sec
ExpVal (Simul) 100.14 ± 0.1 142.11 ± 0.1 159.74 ± 0.1 167.1 ± 0.1
d
T 300 600 900 1200
an
us
cr
ip t
5.2. Gene Expression We consider a simple CRN that models the transcription of a gene into an mRNA molecule, and the translation of the latter into a protein. The CRN, rates and initial conditions are the same as in [41]. The stochastic semantics of the reaction network is an infinite CTMC, and we use this model to show that our method can handle infinite state-space processes. We consider the quantitative property supE=? [#mRN A][T,T ] , which gives the number of molecules of mRN A in the system at time T . We compare our method with SMC estimation of the same property by using 50000 simulations, for T = {300, 600, 900, 1200}, and in the following tables we compare the results in terms of execution time (T ime(·)) and estimated expected value of #mRN A (ExpV al(·)). LNA-based model checking is several orders of magnitude faster without loss of accuracy.
Ac ce p
te
5.3. Fibroblast Growth Factor (FGF) pathway Fibroblast Growth Factors (FGF) are a family of proteins which play a key role in the process of cell signalling in a variety of contexts, like wound healing and skeletal development. We consider the model of FGF signalling pathway developed in [42], which is composed of more than 50 reactions and species. We consider the system with initially 105 molecules for species with non-zero initial concentration. Analysis of the model reveals that the phosphorylated form of F RS2 can bind the protein Src, and then this new complex, Src:F RS2, can relocate out. We want to check if the expected value of #Src:F RS2 during the first 3000 seconds reaches a maximum value greater than 40. We do that by checking the property supE>40 [#Src:F RS2][0,3000] . The formula evaluates to true, and in Figure 2 we analyze the expected value and standard deviation of #Src:F RS2. We obtain these values directly from the logic considering the quantitative interpretation of supE=? [#Src:F RS2][T,T ] and supV=? [#Src:F RS2][T,T ] for T ∈ [0, 3000]. It is possible to see that, after an initial peak, relocation causes exponential decay.
17
Page 17 of 24
ip t cr us an
M
Figure 2: Expected number and standard deviation of species of #Src:F RS2 in the FGF pathway during the first 8000 seconds estimated by our method is compared with a stochastic simulation of the same species.
Ac ce p
te
d
In the same figure we show a single stochastic simulation of the system for the same initial conditions, confirming our evaluation. Moreover, the approximation can be justified theoretically. #Src:F RS2 converges to zero necessarily and this demonstrates the unimodality of the distribution of the species; we note that the variance is finite, so Eqn (4) holds. 5.4. DNA strand displacement of the GW network GW is a network related to the G2-M cell cycle switch [44]. Under particular initial conditions, it has been shown that GW can emulate the Approximate Majority algorithm [43]. Here, we consider the two-domain DNA strand-displacement implementation of GW [3]. The corresponding CRN is composed of 340 species and 240 reactions. For our analysis the species of interest are R and P , whose initial conditions are #R = 90 and #P = 10. These species model the switch for the activating phosphatase Cdc25; initial conditions of other species are taken from the referenced papers. We check the property P>0.6 [#R − #P, [65, +∞]][0,T ] for T = {1000, 2000, 3500, 5000}, in a system of size N = 45000. The results are reported in the following table.
18
Page 18 of 24
Quantitative Value 0.4297 0.5313 0.6535 0.7349
Qualitative Value false false true true
ip t
Time Execution 420 sec 780 sec 1380 sec 2120 sec
cr
T 1000 2000 3500 5000
6. Concluding Remarks
Ac ce p
te
d
M
an
us
We presented a novel probabilistic logic (SEL) for analysing stochastic behaviour of CRNs and proposed an approximate model checking algorithm of the CME based on the LNA. We have implemented the algorithm and demonstrated on four non-trivial examples that LNA-based model checking enables analysis of CRNs with hundreds of species, and even infinite CTMCs, at a cost of some loss of accuracy. It would be interesting to find bounds on the approximation error when the system is far from the thermodynamic limit. However, the error is not only dependent on the value of N , but also on the structure of the CRN, the rates, and the property. For some systems the leap conditions, which guarantee the quality of the approximation if using the LNA, could be satisfied only for a subset of reactions. As recently shown in [45], it is possible to formulate a stochastic hybrid approach, in which the LNA is used only for a subset of species, namely, those for which the leap conditions are satisfied. Other species are treated as a discrete-state Markov process. This improves precision of the stochastic analysis of CRNs when multimodality is present. It would also be interesting to extend SEL with reward operators as in, e.g., [46, 10], with which one can express properties such as the expected number of molecules and variance of a species when certain events happen. We also intend to release a software tool based on LBS [47]. Acknowledgements
The authors would like to thank Luca Bortolussi for helpful discussions.This research is supported by a Royal Society Research Professorship and ERC AdG VERIWARE.
19
Page 19 of 24
References
ip t
[1] V. Chellaboina, S. Bhat, M. Haddad, D. S. Bernstein, Modeling and analysis of mass-action kinetics, Control Systems, IEEE 29 (4) (2009) 60–78.
us
cr
[2] D. Soloveichik, G. Seelig, E. Winfree, DNA as a universal substrate for chemical kinetics, Proceedings of the National Academy of Sciences 107 (12) (2010) 5393–5398.
[3] L. Cardelli, Two-domain DNA strand displacement, Mathematical Structures in Computer Science 23 (02) (2013) 247–271.
M
an
[4] Y.-J. Chen, N. Dalchau, N. Srinivas, A. Phillips, L. Cardelli, D. Soloveichik, G. Seelig, Programmable chemical controllers made from DNA, Nature Nanotechnology 8 (10) (2013) 755–762.
d
[5] D. T. Gillespie, A. Hellander, L. R. Petzold, Perspective: Stochastic algorithms for chemical kinetics, The Journal of Chemical Physics 138 (17) (2013) 170901.
te
[6] D. T. Gillespie, Deterministic limit of stochastic chemical kinetics, The Journal of Physical Chemistry B 113 (6) (2009) 1640–1644.
Ac ce p
[7] L. Cardelli, On process rate semantics, Theoretical Computer Science 391 (3) (2008) 190–215. [8] M. Kwiatkowska, G. Norman, D. Parker, Prism 4.0: Verification of probabilistic real-time systems, in: International Conference on Computer Aided Verification, Springer, 2011, pp. 585–591. [9] M. Kwiatkowska, C. Thachuk, Probabilistic model checking for biology, Software Systems Safety 36 (2014) 165.
[10] M. Kwiatkowska, G. Norman, D. Parker, Stochastic model checking, in: International School on Formal Methods for the Design of Computer, Communication and Software Systems, Springer, 2007, pp. 220–270. [11] J. Elf, M. Ehrenberg, Fast evaluation of fluctuations in biochemical networks with the linear noise approximation, Genome Research 13 (11) (2003) 2475–2484. 20
Page 20 of 24
ip t
[12] L. Bortolussi, R. Lanciani, Model checking Markov population models by central limit approximation, in: International Conference on Quantitative Evaluation of Systems, Springer, 2013, pp. 123–138.
cr
[13] V. Wolf, R. Goel, M. Mateescu, T. A. Henzinger, Solving the chemical master equation using sliding windows, BMC Systems Biology 4 (1) (2010) 42.
us
[14] B. Munsky, M. Khammash, The finite state projection algorithm for the solution of the chemical master equation, The Journal of chemical physics 124 (4) (2006) 044104.
an
[15] A. Singh, J. P. Hespanha, Lognormal moment closures for biochemical reactions, in: Proceedings of the 45th IEEE Conference on Decision and Control, IEEE, 2006, pp. 2063–2068.
M
[16] J. Hespanha, Moment closure for biochemical networks, in: Communications, Control and Signal Processing, 2008. ISCCSP 2008. 3rd International Symposium on, IEEE, 2008, pp. 142–147.
te
d
[17] A. Aziz, K. Sanwal, V. Singhal, R. Brayton, Model-checking continuoustime Markov chains, ACM Transactions on Computational Logic 1 (1) (2000) 162–170.
Ac ce p
[18] C. Baier, B. Haverkort, H. Hermanns, J.-P. Katoen, Model-checking algorithms for continuous-time Markov chains, IEEE Transactions on Software Engineering 29 (6) (2003) 524–541. [19] F. Dannenberg, E. M. Hahn, M. Kwiatkowska, Computing cumulative rewards using fast adaptive uniformization, ACM Transactions on Modeling and Computer Simulation (TOMACS) 25 (2) (2015) 9. [20] L. Bortolussi, L. Cardelli, M. Kwiatkowska, L. Laurenti, Approximation of probabilistic reachability for chemical reaction networks using the linear noise approximation, in: Proc. 13th International Conference on Quantitative Evaluation of SysTems (QEST 2016), LNCS, Springer, 2016, to appear. [21] P. Ballarini, H. Djafri, M. Duflot, S. Haddad, N. Pekergin, COSMOS: A statistical model checker for the hybrid automata stochastic logic, in:
21
Page 21 of 24
ip t
Proc. 8th Int. Conf. Quantitative Evaluation of SysTems (QEST’11), IEEE Computer Society Press, 2011, pp. 143–144.
cr
[22] P. Ballarini, H. Djafri, M. Duflot, S. Haddad, N. Pekergin, Cosmos: a statistical model checker for the hybrid automata stochastic logic, in: Quantitative Evaluation of Systems (QEST), 2011 Eighth International Conference on, IEEE, 2011, pp. 143–144.
an
us
[23] B. Barbot, M. Kwiatkowska, On quantitative modelling and verification of DNA walker circuits using stochastic petri nets, in: R. Devillers, A. Valmari (Eds.), Application and Theory of Petri Nets and Concurrency, Vol. 9115 of Lecture Notes in Computer Science, Springer International Publishing, 2015, pp. 1–32.
M
[24] D. F. Anderson, T. G. Kurtz, Continuous time Markov chain models for chemical reaction networks, in: Design and Analysis of Biomolecular Circuits, Springer, 2011, pp. 3–42.
d
[25] L. Bortolussi, J. Hillston, D. Latella, M. Massink, Continuous approximation of collective system behaviour: A tutorial, Performance Evaluation 70 (5) (2013) 317–349.
te
[26] S. N. Ethier, T. G. Kurtz, Markov processes: characterization and convergence, Vol. 282, John Wiley & Sons, 2009.
Ac ce p
[27] N. G. Van Kampen, Stochastic processes in physics and chemistry, Vol. 1, Elsevier, 1992. [28] E. Cinlar, Introduction to stochastic processes, Courier Corporation, 2013. [29] M. Pinsky, S. Karlin, An introduction to stochastic modeling, Academic press, 2010. [30] F. Didier, T. A. Henzinger, M. Mateescu, V. Wolf, Fast adaptive uniformization of the chemical master equation, in: High Performance Computational Systems Biology, 2009. HIBI’09. International Workshop on, IEEE, 2009, pp. 118–127. [31] E. Wallace, D. Gillespie, K. Sanft, L. Petzold, Linear noise approximation is valid over limited times for any chemical system that is sufficiently large, IET Systems Biology 6 (4) (2012) 102–115. 22
Page 22 of 24
ip t
[32] D. T. Gillespie, The chemical Langevin equation, The Journal of Chemical Physics 113 (1) (2000) 297–306.
cr
[33] J. Elf, M. Ehrenberg, Fast evaluation of fluctuations in biochemical networks with the linear noise approximation, Genome research 13 (11) (2003) 2475–2484. [34] H. Risken, Fokker-Planck Equation, Springer, 1984.
us
[35] J. K. Patel, C. B. Read, Handbook of the normal distribution, Vol. 150, CRC Press, 1996.
an
[36] J. C. Butcher, The numerical analysis of ordinary differential equations: Runge-Kutta and general linear methods, Wiley-Interscience, 1987.
M
[37] K. It¯o, Essentials of stochastic processes, Vol. 231, American Mathematical Soc., 2006. [38] A. Donz´e, O. Maler, Robust satisfaction of temporal logic over realvalued signals, Springer, 2010.
te
d
[39] A. Csik´asz-Nagy, L. Cardelli, O. S. Soyer, Response dynamics of phosphorelays suggest their potential utility in cell signalling, Journal of The Royal Society Interface 8 (57) (2011) 480–488.
Ac ce p
[40] M. Thattai, A. Van Oudenaarden, Intrinsic noise in gene regulatory networks, Proceedings of the National Academy of Sciences 98 (15) (2001) 8614–8619. [41] M.-E.-C. Mateescu, Propagation models for biochemical reaction networks, Ph.D. thesis, EPFL (2011). [42] J. Heath, M. Kwiatkowska, G. Norman, D. Parker, O. Tymchyshyn, Probabilistic model checking of complex biological pathways, Theoretical Computer Science 391 (3) (2008) 239–257. [43] L. Cardelli, Morphisms of reaction networks that couple structure to function, BMC Systems Biology 8 (1) (2014) 84. [44] B. Novak, J. J. Tyson, Numerical analysis of a comprehensive model of M-phase control in Xenopus oocyte extracts and intact embryos, Journal of Cell Science 106 (4) (1993) 1153–1168. 23
Page 23 of 24
ip t
[45] L. Cardelli, M. Kwiatkowska, L. Laurenti, A stochastic hybrid approximation for chemical kinetics based on the linear noise approximation, in: Proc. 14th International Conference on Computational Methods in Systems Biology (CMSB 2016), LNCS, Springer, 2016, to appear.
us
cr
[46] C. Baier, B. Haverkort, H. Hermanns, J.-P. Katoen, On the logical characterisation of performability properties, in: International Colloquium on Automata, Languages, and Programming, Springer, 2000, pp. 780– 792.
Ac ce p
te
d
M
an
[47] M. Pedersen, G. D. Plotkin, A language for biochemical systems: design and formal specification, in: Transactions on computational systems biology XII, Springer, 2010, pp. 77–145.
24
Page 24 of 24