Analysis and Verification of Qualitative Models of Genetic Regulatory Networks: A Model-Checking Approach

1 Introduction

Qualitative simulation aims at making predictions of the behavior of dynamical systems in the absence of detailed, quantitative information on parameter values and functional relations [Kuipers, 1994]. The idea of qualitative simulation has attracted much interest in the context of a biological application, the analysis of genetic regulatory networks, that is, the networks of regulatory interactions between genes, proteins, and small molecules underlying the development and functioning of all living organisms. On the one hand, mathematical methods supported by computer tools are indispensable for the analysis of genetic regulatory networks: most networks of interest involve many genes connected through complex feedback loops, thus making an intuitive understanding of their dynamics difficult to obtain. On the other hand, numerical simulation methods are difficult to apply, because only a few networks are well-understood on the molecular level, and quantitative information on the interactions is seldom available.

We have previously developed a method for the qualitative simulation of genetic regulatory networks [de Jong et al., 2004b]. This method differs from traditional approaches towards qualitative simulation in that it has been tailored to a class of piecewise-linear (PL) differential equations with favorable mathematical properties [Glass and Kauffman, 1973]. This allows it to deal with large and complex networks of regulatory interactions. The qualitative simulation method has been implemented in a publicly-available computer tool, called Genetic Network Analyzer (GNA) [de Jong et al., 2003].

GNA has been used to analyze several genetic regulatory networks of biological interest (e.g., [de Jong et al., 2004a]).

Given a model of a genetic regulatory network, the above qualitative simulation method produces a state transition graph, consisting of qualitative states and transitions between qualitative states. The graph describes the possible qualitative behaviors of the network. When simulating large and complex genetic regulatory networks, the state transition graph usually consists of hundreds or even thousands of states, as a consequence of which the analysis of the graph by visual inspection alone becomes error-prone or even practically infeasible. In order to analyze large state transition graphs, previous work in qualitative simulation has proposed the use of model-checking techniques [Brajnik and Clancy, 1998; Shults and Kuipers, 1997]. These techniques allow for the efficient verification of properties of discrete transition systems [Clarke et al., 1999].

The aim of this paper, extending preliminary ideas in [Batt et al., 2003], is to demonstrate the effectiveness and practical applicability of the combined use of qualitative simulation and model checking for the analysis of large state transition graphs produced by the qualitative simulation of genetic regulatory networks. In particular, we show how state transition graphs can be reformulated as so-called Kripke structures on which properties of the behavior of the network, expressed in temporal logic, can be verified. In addition, we present a new implementation of GNA that connects the qualitative simulator to state-of-the-art model checkers.

We have applied the above approach to the analysis of the complex regulatory network controlling the nutritional stress response in the bacterium Escherichia coli. Among other results, the use of model-checking tools has allowed the characterization of cycles in a large state transition graph obtained through qualitative simulation of the network. The cycles correspond to biologically-important phenomena that are currently being experimentally tested in our laboratory.

In the next section of this paper, we briefly review the qualitative modeling and simulation of genetic regulatory networks. This will set the stage for a discussion of the combined use of qualitative simulation and model-checking techniques in the third section. The application of this approach to the analysis of the nutritional stress response network is the subject of the next section. We finish with a discussion of the approach in the context of related work.

2 Qualitative modeling and simulation of genetic regulatory networks

The dynamics of genetic regulatory networks can be modeled by a class of piecewise-linear (PL) differential equations of the following general form [Glass and Kauffman, 1973]:

where x = (x1, . . . , xn)′ ∈ Ω is a vector of cellular protein concentrations, f =(f1, . . . , fn)′ , g =diag(g1, . . . , gn), and Ω ⊂ Rn≥0 is a bounded n-dimensional phase space box. The rate of change of each protein concentration xi, i ∈ {1..n}, is thus defined as the difference of the rate of synthesis fi(x) and the rate of degradation gi(x) xi of the protein.

The function fi : Ω → R≥0 expresses how the rate of synthesis of the protein encoded by gene i depends on the concentrations x of the proteins in the cell. More specifically, the function fi is defined as a sum of terms having the general form κli bli (x), where κli > 0 is a rate parameter, and bli : Ω → {0, 1} a piecewise-continuous regulation function defined in terms of step functions s+ and s-:

where xi is an element of the state vector x and θi a constant denoting a threshold concentration for xi.

The function gi expressing the regulation of protein degradation is defined analogously, except that we demand that gi is strictly positive. Notice that with the above definitions, h is a piecewise-linear (PL) vector-valued function. An example of a simple genetic regulatory network and its PL model is shown in Figure 1(a)-(b).

The use of step functions s±(xi, θi) in (1) gives rise to complications, because the step functions are discontinuous at xi = θi, and therefore h is discontinuous on Θ = ⋃i∈[1..n],li∈[1..pi]{x ∈ Ω | xi = θli i }, the union of the threshold hyperplanes. In order to deal with this problem, following an approach widely used in control theory, we extend the differential equation (1), defined on Ω \ Θ, to the differential inclusion

defined on Ω, where H(x) is the smallest closed hyperrectangular set containing the set of all limit values of h(y), for y ̸∈ Θ and y → x [Gouz´e and Sari, 2002; de Jong et al., 2004b; Batt et al., 2005]. A solution of the differential inclusion (2) on a time interval I is an absolutely-continuous vector-valued function ξ(t) such that ˙ξ(t) ∈ H(ξ(t)) almost everywhere on I .

Formally, we define the PL system Σ as the triple (Ω,Θ,H), that is, the set-valued function H given by (2), defined on the n-dimensional phase space box Ω, with Θ the union of the threshold hyperplanes. Also, we denote by ΞΣ the set of all solutions ξ(t), on a finite or infinite time interval, of the differential inclusion (2), that reach and leave a threshold hyperplane finitely-many times.

The mathematical analysis of PL systems has shown that mere knowledge of the relative order of threshold parameters and ratios of synthesis and degradation parameters, describing the local strength of gene expression, is sufficient for the definition of a finite partition D of the phase space Ω into hyperrectangular regions, called domains. This partition is such that the derivatives of the solutions in each domain have a unique sign pattern (Figure 1(b)-(d)). The definition of the sign pattern S(x) for a point x in a domain D is complicated by the fact that the solution of a differential inclusion is in general not unique. In particular, if x lies on a threshold hyperplane, H(x) is multi-valued, and several solutions ξ(t) may pass through x with different derivative signs. We therefore define S(x) = {sign( ˙ξ(tx)) | ξ ∈ ΞΣ in D, ξ(tx) = x, and ˙ξ(tx) ∈ H(ξ(tx))}. We have proven that S(x) is the same for every x ∈ D [Batt et al., 2005]. 1

Using the domain partition D of the phase space Ω, together with the above qualitative characterization of the dynamics in each of the domains, we can discretize the continuous dynamics. In the resulting abstract description, the state of the system is represented by a domain, and there exists a transition from a domain D to a domain D′ , if and only if there exists a solution reaching D′ from D, without leaving D ∪ D′ . This leads to the introduction of the so-called qualitative transition system, Σ-QTS = (D,→,Π, |=Π), where D is the set of all domains, →⊆ D × D is a transition relation describing all transitions between the domains, Π = {Dsign = S | S ∈ 2{-1,0,1}n } is a set of propositions describing the derivative sign patterns, and |=Π ⊆ D × Π is a satisfaction relation that associates to each domain a qualitative description of the dynamics of the system in the domain, defined by D |=Π Dsign = S iff for all x ∈ D: S(x) = S.

The graph representation of the qualitative transition system, G = (D,→), is called a state transition graph. The domains are also called qualitative states. A qualitative state D is called instantaneous, if no solution remains in D longer than a time instant, and persistent, otherwise. Moreover, it is a qualitative equilibrium state, if D |=Π Dsign = S, with 0 ∈ S, since then for every x ∈ D, ξ(t) = x is a solution. The state transition graph for the example network is shown in Figure 1(e). A path in the graph describes a possible qualitative behavior of the system (Figure 1(f)).

In [Batt et al., 2005], the qualitative transition system is defined as the discrete quotient of a continuous transition system having the same reachability properties as the original PL system Σ, given the equivalence relation induced by the partition D. Using standard results from hybrid systems theory [Alur et al., 2000], it directly follows that the qualitative transition system Σ-QTS is a conservative approximation of Σ, in the sense that if a solution ξ of the PL systemΣ, defined on some time interval I, passes through the time-ordered sequence of domains (D0 , . . . ,Dm), then (D0 , . . . ,Dm) is a path in Σ-QTS. On the other hand, it may happen that no solution of Σ passes through a sequence of domains corresponding to a path in Σ-QTS. 2

Simple rules exploiting the favorable mathematical properties of the class of PL models have been formulated such that the qualitative transition system can be symbolically computed from a PL model of the network and a set of parameter inequalities defining the relative order of the parameters [Batt et al., 2005]. The implementation of these rules has resulted in a new version of the computer tool Genetic Network Analyzer (GNA), available at http://www-helix.inrialpes.fr/gna.

1 In previous work [de Jong et al., 2004b], a coarser-grained partition of Ω is used, for which this property does not hold.

2 In the terminology of qualitative simulation, the method is sound but incomplete.

3 Analysis of genetic regulatory networks by model checking

Models of genetic regulatory networks of biological interest may produce large state transition graphs, thus making manual analysis of dynamical properties error-prone or even practically infeasible. To address this problem, we propose to combine the qualitative modeling and simulation method outlined above with techniques for model checking [Clarke et al., 1999]. These techniques allow for the efficient verification of properties of the behavior of discrete transition systems, expressed as formulae in some temporal logic.

Various model-checking frameworks exist, differing by their expressiveness, user-friendliness, and computational efficiency. For the sake of clarity of exposition, we focus on one particular framework in this paper, in which the discrete transition system takes the form of a Kripke structure, and the behavioral properties are expressed in Computation Tree Logic (CTL) [Clarke et al., 1999]. However, we emphasize that our approach is not restricted to CTL model-checking, and allows other, more expressive temporal logics to be used as well (Section 3.3).

3.1 Translate qualitative transition system into Kripke structure

As a preliminary step, we introduce a set of atomic propositions, AP, to describe the state of the system. More precisely, the set of atomic propositions we use consists of simple expressions describing a protein concentration (e.g., value xi = θi, value xi < κi/γi), the sign of its derivative (e.g., dsign xi = {1}), or the type of a state (type = pers, type = inst, type = eq, for persistent, instantaneous, and equilibrium state, respectively). The set of atomic propositions for the two-gene example is shown in Figure 2(a).

We define the Kripke structure over AP corresponding to Σ-QTS as Σ-KS = (D,→, |=AP ), where |=AP⊆ D × 2AP is the satisfaction relation associating to every domain D the atomic propositions that hold in D. The existence of at least one solution ξ on some time interval [0, τ ], τ > 0, with initial condition ξ(0) = x0, for every x0 ∈ Ω, guarantees that → is total (i.e., every state has a successor), which is required for Kripke structures. The Kripke structure for the two-gene example is shown in Figure 2(b).

3.2 Express dynamical properties in temporal logic

CTL formulae are interpreted with respect to a Kripke structure. A CTL formula is built upon atomic propositions. The usual operators from propositional logic, such as negation (¬), disjunction (∨), conjunction (∧), and implication (→), can also be used. In addition, CTL provides two types of op-

• (a) AP = { value xa= 0, value xa>0, . . . , dsign xa={-1}, dsign xa={0}, . . . , type=inst, type=eq }.
• (b) Σ-KS = (D,→, |=AP ), with
• D = {D1 ,D2 ,D3 , . . . ,D27 },
• →= { (D1 ,D1 ), (D1 ,D2 ), (D1 ,D11 ), (D1 ,D12 ), . . . , (D27 ,D26 ) },
• |=AP= { (D1 , value xa≥ 0), (D1 , value xa< θ1a), . . . , (D1 , dsign xa={1}), (D1 , dsign xb={1}), (D1 , type=pers), (D2 , value xa= θ1a), . . . , (D2 , type =inst), . . . , (D27 , value xa> θ2a), . . . , (D27 , type =pers) }.

erators: path quantifiers, E and A, and temporal operators, such as F and G. Path quantifiers are used to specify that a property p holds for some (Ep) or every (Ap) path starting from a given state in the Kripke structure. Temporal operators are used to specify that, given a state and a path starting from that state, a property p holds for some (Fp) or for every (Gp) state of the path. Each path quantifier must be paired with a temporal operator [Clarke et al., 1999].

A variety of biologically-interesting properties of the dynamics of a genetic regulatory network can be expressed in CTL [Chabrier-Rivier et al., 2004]. Typical questions concern the possibility or the inevitability for the network to reach some state(s), the absence or the universality of some property, and response or precedence properties (i.e., a property is always followed or preceded by some other property). Some examples are given for the two-gene network in Figure 2(c). More generally, these properties can be used for testing the validity of a model given experimental data, or understanding how the dynamics of a system emerges from the interactions between its components.

3.3 Check if model satisfies dynamical properties

The model-checking problem consists in verifying whether a temporal-logic formula holds in the discrete transition system at the initial state(s). Highly-efficient algorithms have been developed and implemented in model checkers for automatically solving this problem. Two widely-used model checkers are NuSMV [Cimatti et al., 2002] and Evaluator, a component of the CADP toolbox [Mateescu and Sighireanu, 2003]. NuSMV is an efficient, state-of-the-art model checker for CTL, whereas Evaluator is an on-the-fly model checker for the regular alternation-free µ-calculus, a temporal logic based on regular expressions. The algorithmic complexity of these tools is linear in the size of the transition system and the formula. In addition to a yes/no answer, the tools return a diagnostic, either a witness or a counterexample, depending on whether the property holds or not.

In order to combine our qualitative simulator with modelchecking tools, we have integrated export functionalities in GNA. The text files generated by GNA can be imported in the model checkers, after which the verification of the properties of interest continues in the environment of the latter tools. Figure 2(c) illustrates the verification of some properties of the two-gene network. Fairness conditions [Clarke et al., 1999] have been used to ensure that attractors (equilibrium states, strongly connected components) are eventually reached (not shown).

Recall that the qualitative abstraction used to obtain Σ-QTS from Σ is a conservative approximation. Can the conclusions drawn from the verification of a property on Σ-QTS (or more precisely on Σ-KS) be transposed to Σ? If every path in Σ-QTS satisfies a property, then from the conservativeness, a corresponding property holds for every solution of Σ. Conversely, if a property holds for some path in Σ-QTS, then, since this path may be spurious, the corresponding property may hold for no solution of Σ. One can define the set of temporal-logic formulae such that, if a property holds in Σ-QTS, it does so in Σ. In the CTL model-checking framework, this set is the set of formulae that can be expressed in the temporal logic called ∀CTL [Grumberg and Long, 1994].

4 Application: Analysis of the nutritional stress response in E. coli

In case of nutritional stress, an E. coli population abandons exponential growth and enters a non-growth state called stationary phase. This growth-phase transition is accompanied by numerous physiological changes in the bacteria [Huisman et al., 1996]. On the molecular level, this transition is controlled by a complex genetic regulatory network integrating various environmental signals. A network including six genes believed to play a key role in the nutritional stress response is shown in Figure 3(a) [Ropers et al., 2004].

A PL model of seven variables has been constructed, one protein concentration variable for each gene and one input variable representing the presence or absence of a nutritional stress signal. Seven differential equations and forty parameter inequalities describe the dynamics of the system (Figure 3(b)). Using the new version of GNA, we have simulated two phenomena, namely nutrient starvation, causing the

transition from exponential to stationary phase, and nutrient upshift, leading to the reentry into exponential phase. The simulation results have been analyzed by means of the modelchecker NuSMV, using the export functionality of GNA. Below, we give some examples.

The simulation of a nutrient starvation has given rise to a state transition graph of 66 states (27 of which are persistent), computed in less than 1 s. on a PC (800 MHz, 256 Mb). The graph contains a single equilibrium state corresponding to stationary-phase conditions. Are the predictions obtained from the model consistent with experimental data? In [Azam et al., 1999], the measured concentration of the global regulator Fis is shown to decrease and become steady in stationary phase, which is characterized by a low concentration xrrn of stable RNAs:

The verification of (3) takes a fraction of a second to complete and shows that the observed temporal evolution of the Fis concentration is reproduced by the model.

We have also studied the reentry into exponential phase after a nutrient upshift. Using the same model as above, but starting the simulation from the qualitative state characterizing stationary-phase conditions and with the nutritional stress signal switched off, qualitative simulation results in a state transition graph of 1143 states (202 of which are persistent), generated in 1.7 s. The graph contains several strongly connected components (SCCs). One of these SCC, containing 222 states, can be reached from all states in the graph. The qualitative transition system satisfies the property

AG(statesInSCC → AGstatesInSCC ), (4) where the predicate statesInSCC is satisfied by all and only states in the SCC. That is, if the system has reached this SCC, it always remains in it. Checking this property by visual inspection is tedious and error-prone, while it takes NuSMV only 3.3 s. Further mathematical analysis has revealed that the cyclic paths correspond to solutions spiraling inwards to an equilibrium point [Ropers et al., 2004]. In other words, during the reentry into stationary phase, the concentrations of some of the proteins oscillate towards a new equilibrium level. This is a surprising result, which is currently subject to experimental verification in our laboratory. We also checked a dozen of other properties, two of which are shown in Figure 3(c).

5 Discussion

We propose a combination of qualitative simulation and model-checking techniques for the analysis and verification of qualitative models of genetic regulatory networks. Model checkers help in dealing with the problem that the state transition graphs generated by qualitative simulation may become prohibitively large. They permit complex dynamical properties, whose verification defeats manual analysis capabilities, to be efficiently and reliably checked. The analysis of the genetic regulatory network composed of key genes involved in the nutrient stress response in E. coli illustrates the applicability of the approach.

Model-checking techniques have been used before for analyzing biological networks. Most approaches start from discrete models, such as concurrent transition systems [Chabrier-Rivier et al., 2004] and Boolean networks [Bernot et al., 2004]. In this paper we show that model-checking techniques can also be used for more conventional continuous models, in particular differential equation models, when using qualitative abstractions to discretize the dynamics of the system. In comparison with ideas along the same line

[Brajnik and Clancy, 1998; Shults and Kuipers, 1997], our approach is adapted to a particular class of PL differential equations with favorable mathematical properties, allowing the development of tailored algorithms that scale up well to models of large and complex genetic regulatory networks.

The model-validation approach of this paper has been illustrated in the context of CTL model checking. While CTL allows a variety of biologically-meaningful properties to be expressed, some properties fall outside its scope. For instance, in Section 4 we would have liked to be able to express the occurence of oscillations in some of the protein concentrations after a nutrient upshift. The formula EG(p→F¬p ∧ ¬p→Fp) expresses this property, where p means that the concentration of some protein is above a certain threshold. Unfortunately it is not a CTL formula and it does not admit any CTL equivalent [Clarke and Draghicescu, 1988]. However, the property can be expressed in µ-calculus and evaluated using XTL, a component of the CADP toolbox [Mateescu and Garavel, 1998]. The capability of GNA to generate export files for different model checkers, allows one to take advantage from the specific strengths of each of these.