^ IN KI A INSTITUT NATIONAL DE RECHERCHE EN INFORMATIQUE ET EN AUTOMATIQUE Qualitative Simulation of Genetic Regulatory Networks Using Piecewise-Linear Models Hidde de Jong — Jean-Luc Gouze — Celine Hernandez — Michel Page — Tewfik Sari — Johannes Geiselmann N° 4407 Janvier 2002 THEMES 4 et 3 RHONE-ALPES Qualitative Simulation of Genetic Regulatory Networks Using Piecewise-Linear Models Hidde de Jong , Jean-Luc Gouze , Celine Hernandez , Michel Page , Tewfik Sari* , Johannes Geiselmann''' Themes 4 et 3 — Simulation et optimisation de systemes complexes — Interaction homme-machine, images, donnees, connaissances Projets COMORE, HELIX Rapport de recherche n° 4407 — Janvier 2002 — 31 pages Abstract: In order to cope with the large amounts of data that have become available in genomics, mathematical tools for the analysis of networks of interactions between genes, proteins, and other molecules are indispensable. We present a method for the qualitative simulation of genetic regulatory networks, based on a class of piecewise-linear (PL) differential equations that has been well-studied in mathematical biology. The simulation method is well-adapted to state-of-the-art measurement techniques in genomics which often provide qualitative and coarse-grained descriptions of genetic regulatory networks. The method is able to deal with nontrivial mathematical problems induced by the discontinuous right-hand sides of the differential equations. Furthermore, it guarantees that the simulation covers all possible solutions of quantitative PL models corresponding to the qualitative PL model used by the method. The qualitative simulation method has been implemented in Java. Key-words: genetic regulatory networks, mathematical modeling, piecewise-linear differential equations, discontinuous differential equations, qualitative simulation, mathematical biology, bioinformatics * Laboratoire de Mathématiques, Universitě de Mulhouse. t Laboratoire Plasticitě et Expression des Génomes Microbiens, Universitě Joseph Fourier, Grenoble. Unité de recherche INRIA Rhone-Alpes 655, avenue de l'Europe, 38330 Montbonnot-St-Martin (France) Telephone : +33 4 76 61 52 00 — Télécopie +33 4 76 61 52 52 Simulation Qualitative de Reseaux de Regulation Genique en Utilisant des Modeies Lineaires par Morceaux Resume : Afin de traiter les grandes quantites de donnees genomiques disponibles aujourd'hui, des outils mathematiques pour l'analyse de reseaux d'interactions entre genes, proteines, et d'autres molecules sont indispensables. Nous presentons une methode pour la simulation qualitative de reseaux de regulation genique, basee sur une classe d'equations differentielles lineaires par morceaux (LM) qui a ete bien etudiee en biologie mathematique. La methode de simulation est bien adaptee ä l'etat de l'art des techniques experimentales, qui donnent souvent des descriptions qualitatives et agregees des reseaux de regulation genique. La methode est capable de traiter des problemes mathematiques non triviaux induits par des discontinuities dans le second membre des equations differentielles. En outre, eile garantit que la simulation couvre toutes les solutions possibles des modeles LM quantitatifs correspondants au modele LM qualitatif utilise par la methode. La methode de simulation a ete implementee en Java. Mots-cles : reseaux de regulation genique, modelisation mathematique, equations differentielles lineaires par morceaux, equations differentielles discontinues, simulation qualitative, biologie mathematique, bioinformatique Qualitative Simulation of Genetic Regulatory Networks 3 1 Introduction Recent progress in genomics has provided us with experimental tools that hold great promises for unraveling the networks of regulatory interactions between genes, proteins, and small molecules which underlie the functioning of living organisms. On the one hand, these techniques allow protein-DNA and protein-protein interactions to be identified, thus providing insight into the structure of genetic regulatory systems (e.g., [47, 53]). On the other hand, they allow the evolution of the state of the system to be characterized, by large-scale measurement of the level of gene expression and protein activity across time (e.g., [40, 68]). In order to cope with the large amounts of data that have thus become available, formal methods for the representation and analysis of genetic regulatory networks are indispensable. Mathematical models allow networks of interactions to be described in a precise and unambiguous manner, while a large variety of analysis and simulation techniques exists to systematically derive behavioral predictions from the models. The application of formal methods, especially when supported by computer tools, may lead to a comprehension of the structure and functioning of large and complex networks of interactions that cannot be obtained through more intuitive approaches alone [14, 41]. The use of formal methods to study regulatory networks is currently subject to two major constraints [7]. First of all, the biochemical reaction mechanisms underlying the interactions are usually not or incompletely known. This prevents the formulation of detailed kinetic models, such as those developed for the genetic switch controlling phage A growth [42] or the feedback mechanisms regulating tryptophan synthesis in E. coli [55]. A second constraint arises from the general absence of quantitative information on kinetic parameters and molecular concentrations. As a consequence, traditional methods for numerical analysis are difficult to apply. Few of the modeling and simulation methods that have been developed thus far are capable of handling the above constraints. A notable exception is formed by approaches based on a class of piecewise-linear (PL) differential equation models originally proposed by Glass and Kauffman [23]. The state variables in the PL models correspond to the concentrations of proteins encoded by genes in the network, while the differential equations represent the interactions arising from the regulatory influence of some proteins on the synthesis and degradation of others. The regulatory interactions are modeled by means of step functions, which gives rise to the piecewise-linear structure of the differential equations. The use of step functions is motivated by the nonlinear, switch-like character of many of the interactions in gene expression and proteolysis [52, 67]. The PL models provide a coarse-grained description of genetic regulatory networks, well-adapted to state-of-the-art measurement techniques in genomics. Furthermore, the models have mathematical properties that favour qualitative analysis of the steady-state and transient behavior of regulatory systems [11, 12, 19, 20, 24, 25, 38, 44, 45, 48, 56, 57]. On a formal level, the PL models are related to a class of asynchronous logical models proposed by Thomas and colleagues [62, 63]. PL models and their logical relatives have been used for the study of a number of prokaryotic and eukaryotic regulatory networks [1, 10, 18, 43, 46, 51, 54, 61]. In addition, they have been used for modeling food webs [49], neural networks [39], and biological computers [5]. The discontinuous nature of the step functions in PL models brings about some nontrivial mathematical problems. Existing approaches either avoid these problems by restricting the analysis to a subclass of regulatory networks, or adopt solutions that have undesirable consequences for the pre-dictiveness of the method. Recently, it has been shown that an approach capable of dealing with differential equations with discontinuous right-hand sides, widely used in control theory, allows the above-mentioned problems to be resolved in a mathematically proper and practically useful manner [27]. This approach, originally proposed by Filippov [15], is based on the generalization of the differential equations to differential inclusions. In this paper we present a method for the qualitative simulation of genetic regulatory networks described by the generalized PL models. The method is obtained by formulating the analysis of PL RR n° 4407 4 made de Jong et al. models in terms of concepts developed for the qualitative simulation of dynamical systems [10, 35, 36]. In particular, we supplement the PL differential equations with qualitative constraints on the parameter values in the form of algebraic inequalities that can be inferred from biological data. The method exploits the constraints to determine all qualitative states that are reachable from some initial qualitative state through one or more transitions. Each qualitative state corresponds to a region of the phase space where the system behaves in a qualitatively distinct way. The graph of qualitative states and transitions predicted by the method is guaranteed to cover all solutions of the differential equations consistent with the parameter constraints. The simulation method has been implemented in a publicly-available computer tool, called Genetic Network Analyzer (GNA) [9]. In a companion paper we use the method and the tool to analyze a genetic regulatory network of biological interest, consisting of the genes and interactions that regulate the initiation of sporulation in B. subtilis [8]. The application shows that the simulation method can help to gain insight into the qualitative dynamics of complex regulatory networks involving dozens of genes. In the next two sections of the paper, the mathematical framework underlying the simulation method will be reviewed. Sections 4 and 5 introduce the notions of qualitative PL model and qualitative state and behavior, respectively. Using a qualitative PL model of the regulatory network, supplemented by qualitative initial conditions, the simulation algorithm described in section 6 generates a graph of qualitative states and transitions between qualitative states. The soundness and completeness of the algorithm are investigated in section 7. In the final section of the paper, the method is discussed in the context of related work. 2 Piecewise-linear models of genetic regulatory networks Figure 1 shows an example of a simple genetic regulatory network. The genes a and b, transcribed from separate promoters, encode the proteins A and B, each of which controls the expression of both genes.1 Proteins A and B repress gene a and b at different concentrations. Repression of the genes is achieved by binding of the proteins to regulatory sites overlapping with the promoters. The pattern of interactions gives rise to one positive and two negative feedback loops. O A B O Figure 1: Example of a genetic regulatory network of two genes (a and b) coding for a regulatory protein (A and B). The notation follows, in a somewhat simplified form, the graphical conventions proposed by Kohn [34]. The dynamics of genetic regulatory networks can be modeled by a class of differential equations proposed by Mestl et al. [44], extending previous work by Glass and Kauffman [23] (see also the work of Snoussi and Thomas [56, 62] and Ratner and Tchuarev [59]). The equations have the general form &i = fi{x) - 9i{x) Xi, > 0, 1 < i < n, (1) where x = (xi,... ,xn)' is a vector of cellular protein concentrations. The state equations (1) define the rate of change of each concentration X{ as the difference of the rate of synthesis fi(x) and the rate 1As a notational convention, names of genes are printed in italic and names of proteins start with a capital. INRIA Qualitative Simulation of Genetic Regulatory Networks 5 of degradation gi{x)xi of the protein. In vector notation, the system of differential equations (1) is written as x = f(x) - g(x)x, (2) with / = (/i,..., /„)' and g = dmg(gu .. .,gn). The function fi : M"0 —» M>o is defined as /«(«) = J3 ««/&«/(«)> (3) lei where ku > 0 is a (bounded) rate parameter, bu : M"0 —> {0,1} a regulation function defined in terms of step functions, and L a possibly empty set of indices of regulation functions. A regulation function bu is the arithmetic equivalent of a Boolean function expressing the logic of gene regulation [50, 56]. In the simplest case, fi(x) = K{ s+(xj,6j), where the step function s+ : M2 —» {0,1} is defined as follows (see also figure 7): sirj.Oj) J'' Xj>6j> (4) J' " \o, .r, < Oj. As a consequence, gene i is not expressed below a threshold concentration 6j > 0, whereas above this threshold it is expressed at a rate kj. If protein J is a negative regulator of gene i, we have fi(x) = K{ s~(xj,6j), with s~(xj,6j) = 1 — s+(xj,6j). More complex regulation functions can express the combined effect of several regulatory proteins. The use of step functions in gene regulation models has been motivated by the observation that the activation of a gene, as a function of the concentration of a regulatory protein, often follows a steep sigmoidal curve [52, 67]. The function gi allows the regulation of protein degradation to be modeled. It is defined analogously to (3), except that we demand that gi(x) is strictly positive. In addition, in order to formally distinguish degradation rates from synthesis rates, we will denote the former by 7 instead of k. Notice that with the above definitions of fi and gi, the state equations (1) and (2) are piecewise-linear (PL). The PL models can be extended to take into account input variables u = (ui,..., um)', representing the concentration of proteins and small molecules whose synthesis and degradation are regulated outside the svstem. This leads to models of the form: x = f(x,u) -g(x,u)x. (5) In what follows, we will assume that the input variables are constant, i.e., u = 0. As a consequence, (5) can be reduced to (2) without loss of generality, by prior evaluation of the step function expressions in which input variables occur. In figure 2 the state equations for the example network are shown. Gene a is expressed at a rate Ka, if the concentration of protein A is below its threshold d\ and the concentration of protein B below its threshold 6\, that is, if s^(xa, 6%) s^(xt,, 6\) = 1. Analogously, gene b is expressed at a rate Kb, if the concentration of protein A is below the threshold 6\ and the concentration of protein B below the threshold 02. Degradation of the proteins A and B is assumed to be spontaneous, which gives rise to regulation functions having the value 1, independent of the concentrations of the proteins. The dynamical properties of PL models of the form (2) can be analyzed in the n-dimensional phase space box fi = fii x ... x fi„, where every fij, 1 < i < n, is defined as &i = {xi £ IR>o I 0 < Xi < maxi}. (6) maxi is a parameter denoting a maximum concentration for the protein. It can be shown, by generalizing the argument in [21], that if we choose maxi = maxa,>Q fi(x)/gi(x), all trajectories starting RR n° 4407 6 made de Jong et al. State equation for gene a: Xa = Ka S (xa, 0a) S (Xf), Ofo) — 7a %a State equation for gene b: Xh = KbS~(xa,0l) S~(xb,0%) ^JbXb Figure 2: State equations for the network of figure 1. inside fi will remain in it, while trajectories starting outside will enter the phase space box at some point. In general, a protein encoded by a gene will be involved in different interactions at different threshold concentrations, which after ordering are denoted by 6},... ,6f{. The n — 1-dimensional hyperplanes Xi = Of, 1 < ki < pi, divide fi into hyperrectangular regions that are called regulatory domains. Within each region, the concentration of the proteins is bounded by thresholds. More precisely, a regulatory domain D C fi is defined by D = D\ x ... x Dn, where every Di, 1 < i < n, is defined by one of the equations below: Dt ={xi | 0 < Xi < 6}}, Di ={xt | 0} 0, if x(0) = Xq and for all t G [0, r[ it holds that x(t) G D and x(t) = fj,D — vD x(t). For initial values Xq G D, there exists a function x(t) and a r > 0, such that x(t) is the unique solution for (8) on [0, r[. Let 4>i be a function from Ar to fij, defined as i(D) = fif/vf. Analysis of (8) shows that all solution trajectories in D monotonically tend towards a target equilibrium, a stable equilibrium given INRIA Qualitative Simulation of Genetic Regulatory Networks 7 xb Xb maxi-, Kb/Jb — 61 D1 max a $1 i maxa Kahla Xa (a) (b) Figure 3: (a) Phase space box fi for the genetic regulatory network in figure 2. (b) Phase space box with the target equilibrium for the regulatory domain D1 = {(xa, Xb) € ffi2 | 0 < xa < 6^, 0 < xy, < 6\}. The variables xa and Xf, tend towards equilibrium values Ka/ja and «5/75. The target equilibrium is assumed to lie in the upperright domain. by x = 4>{D), with ip = (i,..., n)' [12, 20, 44, 48]. Intuitively speaking, the target equilibrium level i(D) of X{ gives an indication of the strength of gene expression in the regulatory domain. If tf>(D) € D, then for t —» 00 all trajectories in D approach a target equilibrium in D, which represents a so-called regular steady state [57] of the system. If tf>(D) 0 D, all trajectories will at some point leave D. In what follows, we will make the generic assumption that target equilibria are not located in threshold planes. In the example, as can be easily checked from the state equations, we have Ma = {0, Ka}, Na = {ja} for protein A, and Mj, = {0, Nj, = {75} for protein B. In the regulatory domain D1 in figure 3(b), the state equations simplify to Xa = Ka ~ la, Xa, i Xb = Kb ^ Jb Xb- As a consequence, the target equilibrium 4>{Dl) of D1 equals (Ka/ja, Kb/jb)'> which lies outside Dl. The trajectories in D1 will therefore leave the domain at some point. Different regulatory domains generally have different target equilibria. For instance, in the regulatory domain defined by 0 < xa < 6\ and 9\ < Xb < 6%, the target equilibrium is given by (0, Kb/jb)'. 3 Extension of piecewise-linear models to deal with discontinuities The global solution of (2) could be obtained by piecing together the local solutions in regulatory domains, in such a way as to guarantee continuity of the global solution across the threshold hyperplanes [12, 56]. This works fine as long as trajectories arriving at a threshold hyperplane can be continued in another regulatory domain, e.g., trajectories arriving at the threshold hyperplane Xb = Q\ from the regulatory domain D1 (figure 4(a)). However, when the trajectories on both sides of a threshold hyperplane evolve towards this plane, as in the case of trajectories arriving at the threshold hyperplane Xb = 9"l from regulatory domain D3 or D5 (figure 4(b)), mathematical perplexities arise. In the RR n° 4407 8 made de Jong et, al. framework of the previous section, there is no indication on how the local solutions in D3 and D5 can be continued.2 Xb Xb maxb (D3) 01 01 'ßr) D3 \ r>'> r yf LJ- Lj D1 D maxb (D3) el o\ , IP [ n4 \D D3 (D;>) (a) (b) Figure 4: Examples of the behavior of the system of figure 1 at threshold hyperplanes. The figures show the regulatory domains D1,..., D5, and the target equilibria = 02}, which separates the regulatory domains D3 and D5. Notice that xy, is the (only) switching variable in D4. Let D be a switching domain of order k. Let C be the hyperplane of dimension n — k containing D. The boundary of D in C is the set B(D) of all points x G C, such that each ball .Bc^a^e) in C of center x and radius e > 0 intersects both D and C\D [33]. In the case that D is a regulatory domain, C equals fi. Now, for every D G A we define the sets A(£>) = {£>' G A | D' C £(£>)}, and i?(D) = {D' G A | D' regulatory domain, D C B(D')}- A(D) contains the domains in the boundary of D, whereas R(D) contains the regulatory domains that have D in their boundary. In the case of the regulatory domain D1 in figure 4, we find A(Dl) = {-D2, -D6, D7}, while A(D2) = {D7}. Furthermore, R(Dl) = {} and R{D2) = {D1, D3}. The basic idea of the Filippov approach is to extend the differential equations (2) into differential inclusions x G H(x), (10) where H : fi ^ 5(fi) is a set-valued function.4 For a; G D, and D a regulatory domain, we define H(x) simply as H{x) = {fxD ^vDx). (11) Notice that, since the set H(x) contains a single element, the extension of the PL system agrees with the original system in the regulatory domains. If D is a switching domain, H(x) is defined by H{x) = cd({fj,D' ^vD'x | D' G R(D)}). (12) The smallest closed convex set cd(E) of a set E is the intersection of all closed convex sets containing E [15]. In the case of switching domains, H(x) will not generally be single-valued. An absolutely continuous function x(t) = 0, a5o, 0, k,7) is a solution of (10) in the sense of Filippov on [0, t[, t > 0, if x(0) = Xq and for almost all t G [0, r[ it holds that x(t) G H(x(t)) [15]. The qualification 'for almost all t G [0, r[' means that the set of time-points for which the condition does not hold is of measure 0. In particular, the condition is not satisfied at time-points when the solution arrives at a switching domain D, or leaves a switching domain D. If no misunderstanding is possible, we will often simply speak of 'a solution of (10),' instead of 'a solution of (10) in the sense of Filippov.' For all initial values Xq G O there exists a solution of (10) on some [0, r[ [15]. However, this solution is not guaranteed to be unique, due to the generalization of the differential equations to differential inclusions. In order to get an intuitive feeling for the meaning of the above concepts, consider again the examples in figure 4. In the first case, D2 is the switching domain bounding the regulatory domains D1 and D3. If a; € D2, then H(x) is the smallest closed convex set including the end-points of the vectors fxDl — vDl x and fj,D3 — vDi x starting at x. This set is graphically represented in figure 5(a) by the linear segment connecting the end-points of the vectors. H(x) is multiple-valued but, as can be seen, xy, > 0 for any x G H(x), so that a solution trajectory arriving at x crosses the threshold plane instantaneously. In the second example of figure 4, H(x) is a linear segment connecting the end-points of the vectors (j,1*3 — vDi x and fj,D'' — uD'' x. However, this time H(x) intersects with 4For a set E, S(E) represents the set of subsets of E. RR n° 4407 10 made de Jong et, al. D5 Xb (a) (b) Figure 5: Behavior of the system of figure 1 at a point a; in a threshold hyperplane, when the differential equations are generalized into differential inclusions by the method of Filippov. Whereas in (a) solution trajectories cross D2 instantaneously, in (b) they slide along D4. the threshold plane (figure 5(b)). Because the vector fields in D3 and D5 are directed towards D4, a solution trajectory arriving at x continues by sliding along the threshold plane with a derivative vector given by the intersection of H(x) and D4, displaying a so-called sliding mode behavior. For every domain D, a so-called target equilibrium set #(-D) can be defined. If D is a regulatory domain, then <&(£>) = {<£(£>)}. (13) If D is a switching domain, the definition is a little bit more complicated. Let D be a switching domain of order k, contained in the n — fc-dimensional hyperplane C. Then $(D) = C r\cd({(f>(D') | D' € R(D)}). (14) That is, #(-D) is the smallest closed convex set of the target equilibria of regulatory domains D' having D in their boundary, intersected with the hyperplane containing D. As illustrated in figure 5, a solution may (a) instantaneously cross a switching domain or (b) remain in the domain for some time r > 0, sliding along the threshold hyperplane containing the domain. Obviously, the latter case is more interesting when analyzing the behaviour of regulatory systems. Gouze and Sari [27] have shown that the latter sliding mode solutions exist in a switching domain D, iff <&(£>) ^{}. (15) The sliding mode solutions monotonically tends towards the target equilibrium set #(-D) [27]. Because #(-D) does not generally include a single point, the behavior of the system is not uniquely determined by the differential inclusion (10). For example, if a solution trajectory rests in D on the interval [0, r[, it might tend towards one target equilibrium in (D) on the subinterval [0, r'[, and to another target equilibrium in (D) on the subinterval [t',t[ (0 < t' < r). In the special case that #(D)flD ^ {}, there exist solutions in D that asymptotically approach a target equilibrium in D. The equilibrium represents a so-called singular steady state [57] of the system. Whether the equilibrium is stable or unstable must be determined through further analysis. Consider the examples in figure 6. The target equilibrium set #(D2) of the switching domain D2 is defined, following (14), by the intersection of co({(D5)}), the linear segment connecting the points (0, Kb/lb)' and (0,0)', and the threshold boundary Xf, = d\. Consequently, #(D4) equals {(0,6%)'}, and there exists a (unique) sliding mode solution in D4, tending towards (0,6$). Because the target equilibrium lies inside D4, it is also a singular steady state of the system. Closer analysis reveals that the equilibrium (0,02)' is stable, since all trajectories in its neighbourhood end up in this point. xb maxb (D3) 01 01 ) D3 9 u D1 Oz max a Xb maxb (D3) 1>(D4) D4 D3 53(#(D 3)^(ß5)] ) 0 0(ßä) gl Ol maxa Figure 6: Determination of the target equilibrium sets (a) #(D ) = {} and (b) #(D ) = {(0,0^)'}. In order to appreciate the intuitive validity of the Filippov approach, suppose that the step functions s+(xj,6j) in (2) are replaced by sigmoid functions h+(xj,6j,a), defined by h+(xj, 0j,a) x". 0? :i6) 3 3 The step functions and sigmoid functions (called Hill functions) are compared in figure 7. Replacing step functions by sigmoid functions results in differential equation models that are nonlinear but continuous. For steep sigmoid functions, the behavior of the continuous, nonlinear system resembles the behavior of the discontinuous, piecewise-linear system, as shown in figure 8 for the example system (see also [22, 23, 48]). In the continuous case, in the neighbourhood of the plane x\, = 6^, the solution trajectories tend towards an equilibrium located near (0,02)'. A mathematically rigorous comparison of the behavior of discontinuous systems and their continuous homologues can be found in the book of Filippov [15]. (b) Figure 7: Different threshold functions: (a) step function s+(xj,6j) and (b) sigmoid function h+(xj,6j,a). The sigmoid function approaches the step function as a ^s* oo. RR n° 4407 12 made de Jong et, al. Xb D5 D3 xb (h) Figure 8: Behavior of the system of figure 1 near the threshold plane x\> = Q\. In (a) the system is described by the state equations in figure 2, whereas in (b) the system is described by the same state equations, but with the step functions replaced by the sigmoid functions. 4 Qualitative constraints on parameter values Most of the time, precise numerical values for the threshold and rate parameters in (2) are not available. As a consequence, it is not possible to numerically simulate the behavior of a genetic regulatory network. Rather than numerical values, we will specify qualitative constraints on the parameter values. These constraints, having the form of algebraic inequalities, can usually be inferred from biological data. They are exploited by the simulation method to predict the qualitative dynamics of the regulatory system. First of all, we can order the pi threshold concentrations of gene i, yielding the so-called threshold inequalities 0 < 0\ < ... < 9f < maxi. (17) In the case of protein A, there are two threshold concentrations: 6\ is the threshold for the repression of gene b, while d\ is the threshold for the repression of gene a. Assuming the first to be lower than the second, we obtain the threshold inequalities 0 < d\ < d\ < maxa. The ordering of the thresholds of protein B can be determined likewise, giving rise to 0 < 9\ < d\ < maxi,, with 9\ being the threshold for gene a repression and df the threshold for gene b autorepression (figure 9). Second, the possible target equilibrium levels of X{ in different regulatory domains D G ar can be ordered with respect to the threshold concentrations. The resulting equilibrium inequalities define the strength of gene expression in the domain in a qualitative way, on the scale of ordered threshold concentrations. More precisely, for every m € A/;. V{ € iVj, and m,Vi ^ 0, we specify one of the following pairs of inequalities: 0 < m I vi < d}, 6} < ml vi < el (18) 0f < < maxi. The equilibrium inequalities for the example model are shown in figure 9. In the absence of protein B (s~(xb,8l) = 1), while protein A has not yet reached its highest level (s^(xa^6^) = 1), gene a is expressed at a rate na. The corresponding target equilibrium value Ka/ja of xa must be above the second threshold 6%, otherwise the concentration of the protein would not be able to reach or maintain a level at which the observed negative autoregulation of gene a occurs (i.e., d\ < Ka/ja < maxa). INRIA Qualitative Simulation oj Genetic Regulatory Networks 13 State equation for gene a: Xa — Ka s (xai 0a) s a Threshold inequalities: 0 < 0\ < d\ < max a Equilibrium inequalities: d\ < — < max la a State equation for gene b: xh = Khs (xa,9a)s Threshold inequalities: 0 < 6\ < d\ < maxt, xb,&b) - lbXb Equilibrium inequalities: < — < maxb lb Figure 9: State equations, threshold inequalities, and equilibrium inequalities for the network of figure 1. In a similar way, the target equilibrium value Kb/lb is positioned above 6%, again to ensure that the negative autoregulation of gene b at high concentrations of protein B can occur. A quantitative PL model of a genetic regulatory network consists of the state equations (2) and numerical parameter values 0, «,7. In a qualitative PL model, on the other hand, the state equations are supplemented by parameter inequalities (17) and (18). Every quantitative PL model can be abstracted into a unique qualitative PL model. The values of 0 can be uniquely ordered, giving rise to the threshold inequalities (17), while ratios of values of k and 7 can be uniquely ordered with respect to the threshold values, yielding the equilibrium inequalities. Conversely, a qualitative PL model corresponds to a set of quantitative PL models. 5 Qualitative states and behaviors An intuitive qualitative description of the state of a regulatory system consists of the domain in which the system resides, supplemented by the position with respect to this domain of the target equilibrium set to which the state of the system tends. A qualitative behavior is then given by the sequence of qualitative states traversed by the system. This idea will be elaborated below. We first define a function v : AxfJ—»{ —1,0,1}" that maps a domain D and a point e to a sign vector v(D,e) describing the relative position of D and e. Let v(D,e)i denote the ith component of this vector. If X{ is a non-switching variable, then 0 1 1 if Si > sup Di, if inf Di < ei < sup Di if Si < inf D(. % 1 On the other hand, if Xi is a switching variable, then Di = {Of}, and RR n° 4407 14 made de Jong et, al. v(D, e) if ei > Q\, ifOj. if r, < 0>. Generalizing this definition, we obtain the set function V : A x S(Q) domain D and a set E to a set of sign vectors: 1,0,1}") that maps a F(D, £") = {v(D,e) | e G £"}. :i9) Let a?(£) = £(£, 0, a?o,0,K,7) be the solution of a quantitative PL model describing a regulatory network on the time-interval [0, r[. Now suppose that for some t, 0 < i < r, we have a?(£) G D, D G A. The point a?(£) corresponds to a qualitative state of the system defined by QS(aj,t) = <£>,^(£>, $(£>))>. (20) The solution a?(£) on [0, r[ remains in fi (section 2), and hence passes through a sequence of domains D°,..., Dm. The corresponding sequence of qualitative states is called the qualitative behavior of the system on the time-interval. More specifically, a qualitative behavior of the system is defined by QB(x, 0, r) = ((D°, V(D°, $(£>°))>,..., <£>m, V(Dm, $(Dm)))). (21) Consider the solution trajectory in figure 10, obtained for given parameter values, which moves from an initial state in D1 towards a stable equilibrium in D4. Following the above definitions, the solution can be abstracted into a qualitative behavior, as shown in the figure. In the regulatory domain D1, we have #(DX) = {(Ka/ja, Kb/lb)'}- F°r the parameter values in figure 10, we find Ka/ja > Gg, > ®a an^ Kb/lb > > ^s a consequence, F(D1,#(D1)) = {(1,1)}, and we have a qualitative state QSl = (D1, {(1,1)}). In the case of the switching domain D2, the smallest closed convex set of the target equilibria in D1 and D3 consists of the linear segment connecting (Ka/ja, Kb/lb)' and (0, Kb/lb)'■ For the parameter values in figure 10, this segment does not intersect with Xb = 0\, so that F(D2, #(D2)) = {}. The resulting qualitative state is QS2 = (D2, {}). The target equilibrium set of the regulatory domain D3 is given by $(D3) = {(0, Kb/lb)'}, which leads to V(D?\ $(D3)) = {(0,1)}, and QS3 = (D3, {(0,1)}). In the case of D4, the smallest closed convex set containing <^>(-D3) and 4>{Db) consists of the linear segment connecting (0, Kb/lb)' and (0,0)'. This set intersects with Xb = 0\ a* (0,£2)', so that F(D4,$(D4)) = {(0,0)}. We have a qualitative state QS4 = (D4, {(0,0)}). maxb 01 01 (D3) 4>(DV) ■* v U D:s\x r)2 \ S D1 x0 0 ^(ßä) 01 QB(x,0,20) =(QSl, QS2, QS3, QS4) QSl={Dl, {(1,1)}) QS2={D2,{}) QS3={D3, {(0,1)}) QS4={D4, {(0,0)}) 01 max a Figure 10: Abstraction of a solution of the PL model in figure 9 into a qualitative behavior. The solution trajectory has been obtained for the parameter values 6\ = 4, 92 = 8, d\ = 3, 62 = 7, Ka = 20, Kb = 16, 7a = 2, 75 = 2, and the initial conditions Xq = (1.25, 2) on the time-interval [0, 20[. INRIA Qualitative Simulation oj Genetic Regulatory Networks 15 The example shows how a numerical solution of the system can be abstracted into a qualitative behavior consisting of a sequence of qualitative states. More generally, every solution a; of a quantitative PL model on a time-interval [0, r[ can be abstracted into a unique qualitative behavior QB(x, 0, r). The solution traverses a unique sequence of domains D°,..., Dm, while to every domain DJ, 0 < j < m, there corresponds a unique set of sign vectors F(DJ, #(DJ)). Conversely, a single qualitative behavior usually corresponds to a set of solutions. Some properties of the set of solutions abstracting to a qualitative behavior can be directly inferred from the behavior description. Let QB be a qualitative behavior containing a qualitative state (D,V(D,(D))) associated with a domain D G A. Proposition 1 Let D be a switching domain. There exist sliding mode solutions in D, iff V(D, ®(D)) ^ {}• Proof. The condition V(D,(D)) ^ {} means that #(-D) ^ {}. This is exactly condition (15) for the existence of sliding mode solutions in switching domains. □ Proposition 2 There exist solutions in D that asymptotically approach an equilibrium point in D, iff 0 G F(D,$(D)). Proof. If D is a regulatory domain, then (D) = {(D) D D ^ {}. This is the case, iff for all X{, inf D{ < i(D) < supDj, and hence v(D,(D)) from the inequalities (17) and (18) constraining the parameter values. This is a difficult problem in general, because #(-D) may be a complex polyhedron in fi, and the relative position of #(-D) with respect to the threshold planes underdetermined by the parameter inequalities. RR n° 4407 16 made de Jong et, al. Instead of V(D,$>(D)), we will calculate V(D,^(D)), where ^f(D) C fi is the smallest closed hyperrectangle including #(-D). Like <&(-D), ^(D) is a convex set. If D is a regulatory domain, then ^(D) = {(f>(D)}. On the other hand, if D is a switching domain, we obtain V(D) = C Drect ({4>(D') \ D' G R(D) }, (22) where C is the n^fc-dimensional hyperplane containing D. rect(E) is the smallest closed hyperrectangle containing the set E. Although in general ^(D) will be an overapproximation of <&(-D), V(D^(D)) is straightforward to compute from the parameter inequalities and uniquely determined by the latter. In the next section, we will discuss the consequences of the overapproximation of #(-D) by ^(D). In the special case that D is a regulatory domain, ^(D) is a single point in fi, and V(D, ^(D)) = {v(D,(D))}. The equilibrium inequalities tell us whether i(D) > supDj, infDj < i(D) < supDj, or i(D) < infDj, which directly gives the sign vector representing v(D,(D)). Consider the case of V(D1, ^(D1)) = {v(D1,(f>(D1))} in the example system. 4>{Dl) was shown to be (Ka/ja, Kb/jb)'. From the threshold and equilibrium inequalities in figure 9, we know that Ka/ja > Q\ and Kb/lb > Q\- As a consequence, v(Dl, (D'))i represent the ith component of the sign vector v(D,(D')). Now, if X{ is a non-switching variable, then V(D^(D))i is given by V{DMD))i = in e {-1,0,1} I DmmD)v(D,cl>(D% < rt < ^rnax^ v(D, c/>(D% }. (23) On the other hand, if X{ is a switching variable, V(D, ^(D))i is given by V(D, y(D))i = {ri G {-1,0,1} | mm v(D, {D% < rt < max «(£>,<£(£>'))*}• (24) D'eR(D) D'eR(D) As an example, consider the computation of V(D2^(D2)), where D2 is a switching domain and R(D2) = (D1,!)3} (figure 6). By means of the parameter inequalities in figure 9, we find that v(D2,cf>(D1)) = (1,1) and v(D2, (f>(D3)) = (0,1). Using (23)-(24), bearing in mind that Xb is a switching variable in D2, we derive F(D2, \I/(D2)) = {0,1} x {} = {}. A second example concerns the computation of F(D4, \I/(D4)), where D4 is a switching domain and i?(D4) = {D3, D5} (figure 6). By means of the parameter inequalities in figure 9, v(D4, ip(D3)) and v(D4, ip(D5)) are calculated to be (0,1) and (0, -1), respectively. In this case, we find a set F(D4, W(D4)) = {0} x {0} = {(0,0)}. The qualitative states associated to D2 and D4 are (D2,{}) and (D4, {(0,0)}), respectively. 6.2 Computation of state transitions A qualitative state QS = (D, V(D, ^(D))) describes a domain and the position of the target equilibrium set with respect to the domain. Since the trajectories in D tend towards the target equilibrium set, this information can be exploited to determine which domains in the boundary of D, or domains that have D in their boundary, are accessible from D. Since each of these domains is itself associated with a qualitative state, this amounts to computing the possible transitions between qualitative states of the system. The possible transitions are defined by two transition rules. The relative position of the domains D and D' is given by V(D^D'). As can be easily verified, V(D^D') always consists of a single sign vector, that is, V(D^D') = {w}. The domains D and D' are associated with qualitative states QS and QS', calculated to be (D,V(D,^(D))} and {D',V(D',^(D'))}, respectively. INRIA Qualitative Simulation oj Genetic Regulatory Networks 17 Transition rule 1 Let D' G A(D). There is a transition from QS to QS', if 1. F(D,W(D)) + {}, and 2. if X{ is a switching variable in D', but not in D, then there is some v G F(D,\I/(D)), such that = 1- Transition rule 2 Let D € A(D'). There is a transition from QS to if 1. V(D'^(D')) + {}, and 2. if is a switching variable in D, but not in D', then there is some v' G V^D', \I/(D')), such that ^ -1. Intuitively, the first transition rule says that, in order to enter a switching domain D' in the boundary of D, some trajectories must tend towards D' (condition 2). If D is a switching domain, then there must exist sliding mode trajectories in D (condition 1). The second transition rule says that, in order to enter a domain D' from a switching domain D in the boundary of D', the trajectories in D' must not tend towards D (condition 2). If D' is a switching domain, then there must exist sliding mode trajectories in D' (condition 1). We will illustrate the rules by means of the two examples in figure 11, again derived from the model in figure 9. All qualitative states mentioned below have been computed in the way explained in section 6.1. Consider the possible transitions from the qualitative state QS3 associated with regulatory domain D3 to qualitative states associated with the boundary domains A(D3) = {-D2, D4, D7, D8, D9} in (a). We have to verify whether the conditions 1 and 2 of the first transition rule are verified. V(D3,^(D3)) is calculated to be {(0,1)}, while V(D3,D4) equals {(0,1)}. WTith xh a switching variable in D4, but not in D3, we find that conditions 1 and 2 are satisfied. Consequently, there exists a transition from QS3 to QS4. Transitions from QS3 to the other candidate successor states are ruled out, because they violate condition 2. Because D3 is not in the boundary of any domain, the second transition rule cannot be applied. In figure 11(b), the transitions from the qualitative state QS4 = (D4, {(0,0)}) to the states associated with the boundary domain A(D4) = {D9} are investigated. The transition to QS9 = (D9, {}) is excluded, because the condition 2 of the first transition rule is not satisfied. In addition, we consider transitions to QS4 from qualitative states associated with domains that have D4 in their boundary. The second rule is valid for these cases. As can be verified in the figure, D4 G A(D3) and D4 G A(D5). D3 and D5 are regulatory domains, so the condition 1 is trivially satisfied. However, with xy, being a switching variable in D4, but not in D3 and D5, condition 2 is satisfied in neither case. We therefore conclude that there are no transitions from QS4. 6.3 Simulation algorithm The qualitative simulation algorithm can be summarized as follows. Given an initial domain D°, describing the initial protein concentrations Xq , the simulation algorithm computes the initial qualitative state QS°, and then determines all possible transitions from QS° to successor qualitative states by means of the rules above. The generation of successor states is repeated in a recursive manner until all qualitative states reachable from the initial qualitative state have been found. The simulation algorithm is more precisely defined below. The candidate successor states of a qualitative state QS associated with a domain D are the qualitative states associated with domains D', such that D' is in the boundary of D (D' G A(D)) or D is in the boundary of D' (D G A(D')). An actual successor state of QS is a candidate successor state QS' that satisfies the conditions specified in the transition rules. RR n° 4407 18 made de Jong et, al. Xb maxf, 01 n4 j |9 U L D3 L p.2 T 8 D L qs4 qsi ' qs'j • qs' qs2 Xb maxb 01 01 |9 U L D3 qs" qs1 o X - • qs9 qs (b) Figure 11: Possible and impossible (x) transitions from the qualitative states (a) QS3 and (b) QS4, as determined by the transition rules. determine qualitative state QS associated with D°; push(sfae£;, QS°); while not stack is empty do current qualitative state QS <- pop(stack); determine candidate successor states of QS; for every candidate successor QS' do if QS' is actual successor state and not QS' reached before then mark QS' as reached; push(sfae£;, QS') end end end The simulation algorithm generates a directed graph of qualitative states and transitions between qualitative states, a so-called transition graph. Let QS = (D, V(D, ^(D))) be a qualitative state in the graph. If 0 € V(D^(D)), then QS is called a qualitative equilibrium state of the system. Qualitative equilibrium states may correspond to sinks of the graph, although this is not necessarily the case. A set INRIA Qualitative Simulation oj Genetic Regulatory Networks 19 of qualitative states forming a cycle in the graph is called a qualitative cycle of the system. Since the number of possible qualitative states is finite, every path in the transition graph will reach a qualitative equilibrium state or a qualitative cycle at some point. Each path starting at the initial qualitative state QS° and ending at a qualitative equilibrium state or a qualitative cycle forms a possible qualitative behavior predicted by the simulation algorithm. The set of qualitative states included in any path leading to a qualitative equilibrium state or a qualitative cycle is called the attraction set of that state or cycle. Figure 12(a) shows the transition graph for a qualitative simulation of the example system, starting in the regulatory domain D1, where both xa and Xf, lie below their first threshold. As can be seen, the simulation results in five qualitative behaviors leading to different qualitative equilibrium states. In QSlfi, protein A is present at a high concentration (xa = 62), whereas protein B is present at a low concentration (0 < Xf, < 0\). In QS4, protein A is present at a low concentration (0 < xa < 0\) and protein B at a high concentration (xf, = 02). In QS1, protein A and protein B are present at intermediate concentrations (xa = 9\ and x\, = 9\). By simulating the regulatory system for every initial domain, and connecting the resulting transition graphs, the complete phase space behavior of the system can be analyzed. The results for the example system are shown in figure 12(b). They reveal that the system has three qualitative equilibrium states. The attraction set of the qualitative equilibrium state QS4 is indicated in the figure. The simulation results confirm the switch-like character of the network, established by earlier mathematical studies of genetic regulatory networks with the same or a similar structure (e.g., [6, 26, 31, 32, 62, 66]). The qualitative equilibrium states QS4 and QSW correspond to stable equilibria of the system, located at (0,62)' and (62,0)', respectively. The qualitative equilibrium state QS7 corresponds to an unstable equilibrium located at (9\, 0^)'. This equilibrium lies on a separatrix dividing the initial regulatory domain D° into two parts. On one side of the separatrix, the trajectories move towards (92,0)', whereas on the other side of the separatrix the trajectories move towards (0,02)'. Only for the particular case that they start on the separatrix, will the trajectories end up in (9l,9\)'. The three cases are described by the qualitative behaviors (QS1, QS2, QSA, QS4), (QS1, QS6, QS11, QSW), and (QS1, QS7), respectively. The stable equilibria present the two functional states of the system: (1) gene a on and gene b off, (2) gene a off and gene b on. The qualitative behaviors in a transition graph describe how protein concentrations evolve over time. In figure 13, one of the qualitative behaviors of the example system is explored in more detail. The figure shows how the bounds on the concentration of proteins A and B evolve over time, thus displaying the qualitative dynamics of the regulatory system. The simulation method has been implemented in Java 1.3, in a program called GNA (Genetic Network Analyzer) [9].5 The program reads and parses input files specifying the model of the system (state equations, threshold and equilibrium inequalities) and the initial domain. From this information it produces a transition graph. GNA is accessible through a graphical user-interface, which allows the network of interactions between genes to be displayed, as well as the transition graph resulting from the simulation. In addition, the user can analyze the qualitative equilibrium states and qualitative cycles with their attraction sets, and focus on selected qualitative behaviors in order to study the temporal evolution of protein concentrations in more detail. 7 Properties of qualitative simulation Given a qualitative PL model and an initial regulatory domain D°, what can be said about the correctness of the behaviors produced by qualitative simulation? We denote by X the set of solutions x(t) on [0, t[ of all quantitative PL models corresponding to the qualitative model, such that x(0) = Xq € D°. In section 6, the aim of qualitative simulation has been described as finding the set of 5GNA is available for non-profit academic research purposes at http://www-helix.inrialpes.fr/gna. RR n° 4407 20 made de Jong et, al. Xb maxb 4 u LJ L D1 I 6 £,11 i ,16 qs" qs qs Xb rnaxi, 01 |10 jylo £ f9 jjlA £ 20 £)25 19 £>24 U L D3 L P'> T 8 £)13 £ }7 r)12 r IS £)23 17 £>22 D1 I u ±_ 6 £,11 i ,16 £,21 qs:' qs10 qsv' qs20 qs2* -•---•--o (b) Figure 12: Phase spaces and transition graphs for the model describing the network in figure 9. Qualitative states associated with regulatory domains and switching domains are indicated by unfilled and filled dots, respectively. Qualitative equilibrium states are circled in addition, (a) Phase space and transition graph resulting from a simulation of the example system starting in the domain D1. (b) Phase space and transition graph showing all qualitative states of the system and their mutual transitions. The attraction set of the qualitative equilibrium state QS4 is marked by the dotted border. INRIA Qualitative Simulation oj Genetic Regulatory Networks 21 max a 61 j-r maxb n si QS1 QS2 QS3 QS states ■ QS1 QS2 QS3 QS states■ Figure 13: Detailed description of the qualitative behavior (QS1, QS2, QS3, QS4) in figure 12. all qualitative behaviors that abstract from some x G X. More precisely, we would like that (1) for every x G X, the transition graph contains a qualitative behavior QB, such that QB = QB(x,0,t) (soundness), and (2) for every qualitative behavior QB contained in the transition graph, there is some x G X, such that QB = QB(x,0,t) (completeness). 7.1 Soundness Theorem 1 The qualitative simulation algorithm is sound. Proof. Let a; be a solution in X. On [0, r[, x(t) traverses a sequence of domains D°,..., Dm, where Xq G D°. Like in the soundness proof of the qualitative simulation algorithm QSIM [35], we will show by induction that the qualitative behavior QB(x,0,t) is generated by the algorithm. The qualitative state (D°, V(D°, ^(D0))) is generated as described in section 6.1. Since ^(D°) includes #(D°), V(D°, ^(D0)) is guaranteed to include F(D°,#(D0)). A necessary condition for solution trajectories in D° to enter D1 is either that D1 is included in the boundary of D°, i.e. Dl G A(D°), or that D° is included in the boundary of D1, i.e. D° G A(Dl). (1) If D1 G A(D°), then the transition from the qualitative state associated with D° to the qualitative state associated with D1 is generated by means of transition rule 1. For, in order to reach the switching domain D1, x must remain in D° for some time. That is, D° must be a regulatory domain or x must be a sliding mode solution in D°. Moreover, the solution trajectory in D° must tend towards D1. Together with F(D°,W(D0)) D F(D°,$(D0)), this implies that the two conditions of transition rule 1 are satisfied. (2) If D° G A(Dl), then the transition from the qualitative state associated with D° to the qualitative state associated with D1 is generated by means of transition rule 2. For, in order to enter D1 from D°, x must remain in D1 for some time. That is, D1 must be a regulatory domain or x must be a sliding mode solution in D1. Moreover, the solution trajectory in D1 must point away from D°. As a consequence, with V(D°, ^(D0)) D V(D°, #(D0)), it follows that the two conditions of transition rule 2 are satisfied. We conclude that the transition from the qualitative state associated with D° to the qualitative state associated with D1 is generated by the simulation algorithm. Now suppose that the qualitative state associated with the domain DJ, 0 < j < m, has been generated. By an argument parallel to that given above, it can be shown that the qualitative state associated with the domain DJ+1 will be generated by the algorithm as well. Consequently, the qualitative behavior QB(x,0,t) is generated step by step. □ The soundness of the qualitative simulation algorithm in section 6 has important consequences. It implies that, whatever the exact numerical values for the parameters may be, if these values are consistent with the threshold and equilibrium inequalities specified in the qualitative PL model, the solutions of the quantitative PL model are covered by the transition graph generated by the algorithm. More in particular, an equilibrium attained by some solution of a quantitative PL model consistent with the parameter inequalities necessarily corresponds to some qualitative equilibrium state in the transition graph. For, by proposition 2, we know that, if domain D contains an equilibrium, then 0 G RR n° 4407 22 made de Jong et, al. V(D,$>(D)). Soundness of the simulation algorithm guarantees that the qualitative state associated with D is actually generated, and that 0 € V(D,^(D)). 7.2 Incompleteness On the other hand, the qualitative simulation algorithm is not complete, as will be illustrated by means of a counterexample. The network in figure 14 consists of two genes, a and b, which encode the proteins A and B, respectively. The proteins form a heterodimer A-B repressing the expression of both genes. The qualitative PL model is shown in figure 15. It is assumed that the heterodimer represses the two genes above the same threshold concentration A-B -O A r b i5 nS U L D1 1 u ,4 D7 maxa Xn, rnaxb 01 4,(1) QS-* QS2 QS1 ^D ) Ola maxa QS" *® QSB QS:' • —-o QS4 QS7 (a) (b) (c) Figure 16: (a) Phase space box and transition graph for the genetic regulatory network described by the model in figure 15. (b) Determination of \I/(D6). (c) Transition graph resulting from a simulation of the system starting in D1. In order to explain this result, consider the computation of the qualitative state associated with the switching domain D6. R(D6) = {D3,D9}, and 0(D3) = (Ka/ja,Kh/jhy and 4>{D9) = (0,0)'. Using the equilibrium inequalities in figure 15, we compute V(D%, \I/(D6)) as described in section 6.1. We find V(D%, \I/(D6)) = {(0, —1), (0,0)}, where \I/(D6) is the intersection of the smallest rectangle including 4>{DZ) and (D6), as can be seen in figure 17. The smallest closed convex set co({(D9)}) is not red ({(D9)}), but rather the linear segment connecting the two target equilibria. As a consequence, #(D6) is a single point in fi, located at the intersection of the linear segment connecting <^>(-D3) and (D6)) = {(0,0)}, in (b) V(D6, 3>(D6)) = {(0,-1)}, and in (c) V(D6, 3>(D6)) = {(0,-1)}. In other words, the parameter inequalities in figure 15 are consistent with two different qualitative states for QS6, both of which are covered by the qualitative state actually determined by the algorithm. Each configuration in figure 17(a)-(c) implies an additional constraint on the parameter values. In particular, we have (a) ^ > % (b) ^ = % and (c) ^ < % Kalb 6k Kalb #a Ka lb &k With the above constraints, we should have obtained the three different transition graphs shown in figure 17(d)-(f). Each of these is a subgraph of the transition graph actually generated by the simulation algorithm. Notice that the transition graphs in figure 17 contain a single qualitative equilibrium state, which is a sink of the graph. The example illustrates that the transition graph resulting from a simulation may contain more qualitative behaviors than necessary to cover the trajectories corresponding to solutions in X. As a consequence, the algorithm is incomplete. Among other things this implies that, although we are sure to never miss a qualitative equilibrium state, not all qualitative equilibrium states in the transition graph need to correspond to some equilibrium of a quantitative PL model corresponding to the qualitative PL model. The incompleteness of the simulation algorithm raises issues for further research that will be addressed in the next section. RR n° 4407 24 made de Jong et, al. qs4 qs7 qs4 qs7 qs4 (d) (e) (f) Figure 17: (a)-(c) Depending on the exact position of {DZ). 8 Discussion We have presented a method for the qualitative simulation of genetic regulatory networks described by a class of piecewise-linear (PL) differential equations that has been well-studied in mathematical biology. The method allows the behavior emerging from large and complex networks of genetic regulatory interactions to be predicted in a qualitative manner. In the accompanying paper, we describe a model of the network underlying the initiation of sporulation in B. subtilis, and we compare predictions obtained through simulation with observations of the behavior of wild type and mutant bacteria [8]. The application of the qualitative simulation method is supported by a computer tool, called Genetic Network Analyzer (GNA) [9]. The PL models employed in this paper are based on step function approximations of the regulatory interactions involved in the synthesis and degradation of proteins. The step functions provide a succinct description of the regulatory logic, while abstracting from the details of molecular interactions. The biological validity of the step function expressions derives from experimental evidence that the activation of a gene, as a function of the concentration of a regulatory protein, often follows a steep sigmoidal curve [52, 67]. That is, below a certain threshold concentration of the protein, the gene will be hardly expressed at all, whereas above this threshold its expression rapidly saturates. Recent experimental studies have shown that some aspects of the qualitative dynamics of genetic regulatory networks synthesized in vivo correspond well with the predictions obtained from mathematical models based on switch-like approximations of regulatory interactions [4, 13, 17]. The use of step functions gives rise to discontinuities in the right-hand side of the differential equations, which may lead to nontrivial mathematical problems, as illustrated in figure 4. Several ways INRIA Qualitative Simulation oj Genetic Regulatory Networks 25 to deal with the step function discontinuities have been proposed in the literature. The application of the PL models can be restricted to systems without autoregulation, which excludes situations of the type described in figure 4(b) [20]. Alternatively, when a trajectory arrives at a switching domain from which it cannot be continued, it may simply be stipulated to come to a dead stop [48]. Another solution, based on an idea of Plahte et al. [44, 48], consists of avoiding the discontinuities altogether, by replacing the step functions s+(xj,6j) by so-called logoid functions l+(xj,6j,S) that monotonically increase from 0 to 1 in a 5-interval around the threshold 6j. The logoid functions approach their step function homologues as S —» 0. Each of the above solutions is unsatisfactory in some way. In the first place, autoregulation is an ubiquitous feature in genetic regulatory networks [60]. (In fact, the network shown in figure 1 is a simplified version of the molecular switch determining the response of E. coli to phage A infection [52].) Ignoring trajectories that cannot be continued in a switching domain will cause important behavioral properties of regulatory systems to be missed, like the singular steady states in the example. The use of logoid functions instead of step functions is attractive at first sight, but leads to nonlinear differential equation models that are difficult to treat in a qualitative way (see an earlier version of the simulation method in [10]). Here, we have adopted another solution, based on an approach to deal with differential equations with discontinuous right-hand sides originally proposed by Filippov [15]. This approach, recently applied to PL models of the form (2), has the advantage of putting no restrictions on the class of genetic regulatory networks that can be handled, while explicitly defining the behavior of the system in the threshold planes by means of simple-to-analyze PL models [27]. The simulation method described in this paper allows genetic regulatory networks to be simulated in the absence of numerical values for the kinetic parameters and molecular concentrations. The possible qualitative behaviors of the system are determined by means of qualitative constraints on the threshold and rate parameters. Whereas exact numerical values for the parameters are usually not available, the choice of appropriate threshold and equilibrium inequalities can be based on biological data, or is at least strongly constrained by the latter. If the choice of parameter inequalities is not unambiguously determined by the data, the consequences of opting for one combination of inequalities rather than another can be explored by simulating the system for each of the alternatives. Notice that in this case, instead of continuous ranges of numerical values, only a limited and often restricted number of alternative parameter inequalities need to be considered. The underdetermination of initial conditions by available biological data can be handled in the same way. The predictions obtained through qualitative simulation have the form of a graph of qualitative states and transitions between qualitative states. Each qualitative state in the graph corresponds to a domain in the phase space where the system behaves in a qualitatively homogeneous manner, while a transition between two qualitative states corresponds to solution trajectories that leave one domain and enter another. Qualitative simulation makes predictions about the qualitative shape of solutions, illustrated in figure 13, as opposed to the quantitative predictions obtained through numerical simulation. The qualitative nature of the predictions is well-adapted to current measurement techniques in genomics, which have limited quantitative precision, but are able to detect qualitative changes in gene expression over time. The soundness property of qualitative simulation guarantees that no solution of a quantitative PL model consistent with the qualitative PL model is omitted. In fact, to each such solution there corresponds a qualitative behavior in the transition graph. The soundness of qualitative simulation provides for a straightforward check of the robustness of simulation results to changes in parameter values. If a certain behavior is not observed for given parameter inequalities, one can be sure that it will not occur for any of the parameter values consistent with the threshold and equilibrium inequalities. That is, the behavioral variety of the system, given the constraints expressed in the parameter inequalities, is completely covered by the transition graph. This is a powerful result that may guide computationally-expensive simulations to test the sensitivity of the results to changes in parameter RR n° 4407 26 made de Jong et, al. values, as for example carried out in studies of the segment polarity network in Drosophila by von Dassow et al. [65] and simple metabolic pathways by Alves and Savageau [2, 3]. The lack of quantitative information on kinetic parameters and molecular concentrations has stimulated an interest in methods for modeling and simulation developed in the field of qualitative reasoning (QR), most notably QSIM [35] and QPT [16]. QR methods have been applied to the regulation of tryptophan synthesis [29] and A phage growth [28] in E. colt, and to the regulation of the transcription factor families ap-1 and NF-kB in different classes of animals [64]. A major problem with existing QR methods is their lack of upscalability, which causes the applicability of the methods to be limited to small regulatory systems of modest complexity. As its application to the sporulation example in the accompanying paper shows, the qualitative simulation method presented here is able to deal with large and complex networks. Upscaling of the method has been achieved by using PL models that strongly constrain the local dynamics of the system. Moreover, the representation of qualitative states and the transition rules have been tailored to this class of models, in order to maximally exploit their favorable mathematical properties. Thus, in comparison with existing qualitative simulation methods like QSIM and QPT, the expressivity and generality of the modeling formalism have been traded for upscalability. This is particularly clear for the description of qualitative states, which is achieved on a higher level of abstraction in our method. Domains containing a hyperplane x{ = i(D) = jvf give rise to a single qualitative state in our case, whereas in a method like QSIM several qualitative states must be defined, corresponding to the part of the domain located above, below, and inside the hyperplane [35]. Qualitative methods for the analysis of genetic regulatory systems have been developed in mathematical biology as well, the best-known example being Boolean networks [30, 58]. Simulation of Boolean networks rests on the assumption that a gene is either active or inactive, and that genes change their activation state synchronously. For the purpose of modeling actual genetic regulatory networks, these assumptions are usually too strong. Thomas and colleagues [62, 63] have proposed a generalized logical method that permits multivalued activation states and asynchronic transitions. On the formal level, the method of Thomas and colleagues is related to the approach presented in this paper. In fact, Snoussi has demonstrated that the logical equations can be interpreted as an abstraction of a special case of (1), where in the production term fi(x) = Y^iel Kn ^u(x) ^ holds that either bu(x) = s+(xj,6j) or bu(x) = s~(xj,6j), while in the degradation term gi(x) = Y^ieLlubui®) it holds that bu(x) = 1 [56]. Although some ideas of the generalized logical method have been retained in the method presented here, in particular the parameter constraints of section 4, which are related to the logical parameters in [62, 63], we have opted for differential equation models. We believe that the latter formalism is intuitively clear and of large generality. In particular, it allows for a transparent description of the behavior of the system in the threshold planes. Although for the class of PL models covered by the generalized logical method, certain patterns of logical states can be interpreted as indicating singular steady states of the system [57], a general description of the behavior of the system in the threshold planes is currently missing. The differential equation formalism has the additional advantage of facilitating the integration of any quantitative data becoming available through improvements of current measurement techniques. The simulation of genetic regulatory networks may lead to unwieldy transition graphs that are difficult to interpret. The availability of computer tools, allowing one to identify and to focus upon interesting parts of the transition graph, forms a prerequisite for the simulation of large and complex regulatory systems [9]. In addition, a weaker notion of soundness might be chosen. We call an algorithm almost sound, if for every x € Y C X, the transition graph contains a qualitative behavior QB, such that QB = QB(x,0,t), while for every x € X/Y, the transition graph does not contain a qualitative behavior QB, such that QB = QB(x,0,t). Moreover, the set X/Y is of measure 0 with respect to the set X. That is, the size of the set X/Y of trajectories not covered by a qualitative behavior in the transition graph is negligible in comparison with the size of the total set of trajectories X. INRIA Qualitative Simulation oj Genetic Regulatory Networks 27 The simulation algorithm could be made almost sound by adapting the transition rules in section 6.2. For example, transitions in which the order of the domain of the qualitative states changes by more than one might be excluded. In fact, these transitions correspond to sets of trajectories lying on a separatrix in a domain, like in the transition from QSl to QS1 in the example of figure 12. Although it still needs to be proven that the resulting algorithm is almost sound, there is some intuitive support for omitting these transitions. Biologically speaking, they represent the rare event that two or more protein concentrations reach a threshold at exactly the same moment [62]. Excluding these transitions would have the advantage of reducing the connectivity of the transition graph, and thus increasing their legibility, while retaining the transitions that are essential. The simulation algorithm is sound but incomplete, which means that the transition graph may contain qualitative behaviors that do not correspond to any solution of a quantitative PL model satisfying the parameter constraints in the qualitative PL model. An illustration of this point was given in section 7.2, when discussing the results of the simulation of the genetic regulatory network in figure 14. The example reveals that the transition graph generated by the algorithm should be refined into three subgraphs, corresponding to alternative implicit additional constraints on the parameter values. This conclusion required graphical and algebraic considerations that are far beyond the scope of the current simulation algorithm. An issue for further research would be to exploit methods for the symbolic analysis of sets of algebraic inequalities, in order to obtain more precise estimates of the position of a target equilibrium set #(-D) with respect to its domain D. Proposition 2 expresses a relation between the existence of an equilibrium in a domain D and the qualitative state associated with D. Previous work on PL models of the form (2) has led to theorems that can be reformulated as stating a relation between the existence of a limit cycle in fi and a certain type of qualitative cycle in the transition graph [24, 25, 45, 56]. Given the complexity of a general analysis of the occurrence of limit cycles for equations of the form (2), the relation is much weaker than its counterpart for equilibrium points. An interesting issue for further work would be to try to generalize the relation between limit cycles and qualitative cycles, using the literature on periodic solutions in differential equations with discontinuous right-hand sides (see [37] for a review). The qualitative simulation method provides predictions of the possible qualitative behaviors of a genetic regulatory network. The interest of these predictions is that they can be directly compared with gene expression profiles obtained by means of quantitative RT-PCR or DNA microarrays. The use of predicted qualitative behaviors in combination with observed gene expression profiles allows hypothesized models of regulatory networks to be rapidly tested, even when only imprecise data is available. Along these lines, we are currently working on extensions of the method to validate and identify models of genetic regulatory networks using gene expression data. Incorporation of these extensions in the computer tool described in section 6.2 would allow the simulation method to evolve into a more general approach towards the computer-supported analysis of genetic regulatory networks. Acknowledgments The authors would like to thank Gregory Batt for comments on previous versions of this report. References [1] R. Alur, C. Belta, F. Ivančic, V. Kumar, M. Mintz, G.J. Pappas, H. Rubin, and J. Schlug. Hybrid modeling and simulation of biomolecular networks. In M.D. Di Benedetto and A. Sangiovanni-Vincentelli, editors, Hybrid Systems: Computation and Control (HSCC 2001), volume 2034 of Lecture Notes in Computer Science, pages 19 32. Springer-Verlag, Berlin, 2001. [2] R. Alves and M.A. Savageau. Comparing systemic properties of ensembles of biological networks by graphical and statistical methods. Bioinformatics, 16(6):527^533, 2000. RR n° 4407 28 made de Jong et, al. [3] R. Alves and M.A. Savageau. Systemic properties of ensembles of metabolic networks: Application of graphical and statistical methods to simple unbranched pathways. Bioinformatics, 16(6):534^ 547, 2000. [4] A. Becskei and L. Serrano. Engineering stability in gene networks by autoregulation. Nature, 405:590-591, 2000. [5] A. Ben-Hur and H.T. Siegelmann. Computation in gene networks. In M. Margenstern and Y. Rogozhin, editors, Machines, Computations, and Universality: Third International Conference (MCU 2001), volume 2055 of Lecture Notes in Computer Science, pages 11 21. Springer-Verlag, Berlin, 2001. [6] J.L. Cherry and F.R. Adler. How to make a biological switch. Journal of Theoretical Biology, 203:117-133, 2000. [7] H. de Jong. Modeling and simulation of genetic regulatory systems: A literature review. Journal of Computational Biology, 9(1):69-105, 2001. [8] H. de Jong, J. Geiselmann, G. Batt, C. Hernandez, and M. Page. Qualitative simulation of the initiation of sporulation in B. subtilis. in preparation, 2002. [9] H. de Jong, J. Geiselmann, C. Hernandez, and M. Page. Genetic Network Analyzer: A tool for the qualitative simulation of genetic regulatory networks. Technical Report RR-4262, INRIA Rhone-Alpes, Montbonnot Saint-Martin, 2001. [10] H. de Jong, M. Page, C. Hernandez, and J. Geiselmann. Qualitative simulation of genetic regulatory networks: Method and application. In B. Nebel, editor, Proceedings of the Seventeenth International Joint Conference on Artificial Intelligence, IJCAI-01, pages 67^73, San Mateo, CA, 2001. Morgan Kaufmann. [11] R. Edwards and L. Glass. Combinatorial explosion in model gene networks. Chaos, 10(3):691-704, 2000. [12] R. Edwards, H.T. Siegelmann, K. Aziza, and L. Glass. Symbolic dynamics and computation in model gene networks. Chaos, 11(1):160-169, 2001. [13] M.B. Elowitz and S. Leibler. A synthetic oscillatory network of transcriptional regulators. Nature, 403:335^338, 2000. [14] D. Endy and R. Brent. Modelling cellular behavior. Nature, 409:391^395, 2001. [15] A.F. Filippov. Differential Equations with Discontinuous Righthand Sides. Kluwer Academic Publishers, 1988. [16] K.D. Forbus. Qualitative process theory. Artificial Intelligence, 24:85^168, 1984. [17] T.S. Gardner, C.R. Cantor, and J.J. Collins. Construction of a genetic toggle switch in Escherichia coli. Nature, 103:339 312. 2000. [18] R. Ghosh and C.J. Tomlin. Lateral inhibition through Delta-Notch signaling: A piecewise affine hybrid model. In M.D. Di Benedetto and A. Sangiovanni-Vincentelli, editors, Hybrid Systems: Computation and Control (HSCC 2001), volume 2034 of Lecture Notes in Computer Science, pages 232 2 10. Springer-Verlag, Berlin, 2001. [19] L. Glass. Classification of biological networks by their qualitative dynamics. Journal of Theoretical Biology, 54:85^107, 1975. INRIA Qualitative Simulation of Genetic Regulatory Networks 29 39 L. Glass. Combinatorial and topological methods in nonlinear chemical kinetics. Journal of Chemical Physics, 63(4):1325^1335, 1975. L. Glass. Global analysis of nonlinear chemical kinetics. In B. Berne, editor, Statistical Mechanics, Part B: Time Dependent Processes, pages 311 319. Plenum Press, New York, 1977. L. Glass and S.A. Kauffman. Co-operative components, spatial localization and oscillatory cellular dynamics. Journal of Theoretical Biology, 3 1:219 237. 1972. L. Glass and S.A. Kauffman. The logical analysis of continuous non-linear biochemical control networks. Journal of Theoretical Biology, 39:103^129, 1973. L. Glass and J.S. Pasternack. Prediction of limit cycles in mathematical models of biological oscillations. Bulletin of Mathematical Biology, 10:27 11. 1978. L. Glass and J.S. Pasternack. Stable oscillations in mathematical models of biological control systems. Journal of Mathematical Biology, 6:207^223, 1978. B.C. Goodwin. Temporal Organization in Cells. Academic Press, New York, N.Y., 1963. J.-L. Gouze and T. Sari. A class of piecewise linear differential equations arising in biological models. Technical Report RR-4207, INRIA Sophia-Antipolis, Sophia-Antipolis, 2001. K.R. Heidtke and S. Schulze-Kremer. Design and implementation of a qualitative simulation model of A phage infection. Bioinformatics, 14(1):81—91, 1998. P.D. Karp. Design methods for scientific hypothesis formation and their application to molecular biology. Machine Learning, 12:89-116, 1993. S.A. Kauffman. The Origins of Order: Self-Organization and Selection in Evolution. Oxford University Press, New York, 1993. A.D. Keller. Specifying epigenetic states with autoregulatory transcription factors. Journal of Theoretical Biology, 170:175-181, 1994. A. D. Keller. Model genetic circuits encoding autoregulatory transcription factors. Journal of Theoretical Biology, 172:169^185, 1995. J.L. Kelley. General Topology. Van Nostrand, New York, 1969. K.W. Kohn. Molecular interaction maps as information organizers and simulation guides. Chaos, 11(1):1-14, 2001. B. Kuipers. Qualitative Reasoning: Modeling and Simulation with Incomplete Knowledge. MIT Press, Cambridge, MA, 1994. B.J. Kuipers. Qualitative reasoning: Modeling and simulation with incomplete knowledge. Auto-matica, 25(4):571^585, 1989. R.I. Leine, D.H. van Campen, and B.L. van de Vrande. Bifurcations in nonlinear discontinuous systems. Nonlinear Dynamics, 23:105^164, 2000. J.E. Lewis and L. Glass. Steady states, limit cycles, and chaos in models of complex biological networks. International Journal of Bifurcation and Chaos, 1 (2): 177 1*3. 1991. J.E. Lewis and L. Glass. Nonlinear dynamics and symbolic dynamics of neural networks. Neural Computation, 1:621 6 12. 1992. RR n° 4407 30 tlidde de Jong et, al. [40] D.J. Lockhart and E.A. Winzeler. Genomics, gene expression and DNA arrays. Nature, 405:827^ 836, 2000. [41] H.H. McAdams and A. Arkin. Simulation of prokaryotic genetic circuits. Annual Review of Biophysics and Biomolecular Structure, 27:109 22 I. 1998. [42] H.H. McAdams and L. Shapiro. Circuit simulation of genetic networks. Science, 269:650^656, 1995. [43] L. Mendoza, D. Thieffry, and E.R. Alvarez-Buylla. Genetic control of flower morphogenesis in Arabidopsis thaliana: A logical analysis. Bioinformatics, 15(7-8) :593^606, 1999. [44] T. Mestl, E. Plahte, and S.W. Omholt. A mathematical framework for describing and analysing gene regulatory networks. Journal of Theoretical Biology, 176:291^300, 1995. [45] T. Mestl, E. Plahte, and S.W. Omholt. Periodic solutions in systems of piecewise-linear differential equations. Dynamics and Stability of Systems, 10(2):179—193, 1995. [46] S.WT. Omholt, X. Kefang, 0. Andersen, and E. Plahte. Description and analysis of switchlike regulatory networks exemplified by a model of cellular iron homeostasis. Journal of Theoretical Biology, 195:339^350, 1998. [47] A. Pandey and M. Mann. Proteomics to study genes and genomes. Nature, 405:837^846, 2000. [48] E. Plahte, T. Mestl, and S.WT. Omholt. Global analysis of steady points for systems of differential equations with sigmoid interactions. Dynamics and Stability of Systems, 9( 1):275 291. 1994. [49] E. Plahte, T. Mestl, and S.WT. Omholt. Stationary states in food web models with threshold relationships. Journal of Biological Systems, 3(2):569 577. 1995. [50] E. Plahte, T. Mestl, and S.WT. Omholt. A methodological basis for description and analysis of systems with complex switch-like interactions. Journal of Mathematical Biology, 36:321 3 IK. 1998. [51] E.I. Prokudina, R.Y. Valeev, and R.N. Tchuraev. A new method for the analysis of the dynamics of the molecular genetic control systems. II. Application of the method of generalized threshold models in the investigation of concrete genetic systems. Journal of Theoretical Biology, 151:89-110, 1991. [52] M. Ptashne. A Genetic Switch: Phage A and Higher Organisms. Cell Press k, Blackwell Science, Cambridge, MA, 2nd edition, 1992. [53] B. Ren, F. Robert, J.J. Wyrick, O. Aparicio, E.G. Jennings I. Simon, J. Zeitlinger, J. Schreiber, N. Hannett, E. Kanin, T.L. Volkert, C.J. Wilson, S.P. Bell, and R.A. Young. Genome-wide location and function of DNA binding proteins. Science, 290:2306^2309, 2000. [54] L. Sanchez and D. Thieffry. A logical analysis of the Drosophila gap genes. Journal of Theoretical Biology, 211:115 1 11. 2001. [55] M. Santillan and M.C. Mackey. Dynamic regulation of the tryptophan operon: A modeling study and comparison with experimental data. Proceedings of the National Academy of Sciences of the USA, 98(4):1364^1369, 2001. [56] E.H. Snoussi. Qualitative dynamics of piecewise-linear differential equations: A discrete mapping approach. Dynamics and Stability of Systems, 4(3-4):189^207, 1989. [57] E.H. Snoussi and R. Thomas. Logical identification of all steady states: The concept of feedback loop characteristic states. Bulletin of Mathematical Biology, 55(5):973 991. 1993. INRIA Qualitative Simulation of Genetic Regulatory Networks 31 [58] R. Somogyi and C.A. Sniegoski. Modeling the complexity of genetic networks: Understanding multigenic and pleiotropic regulation. Complexity, 1(6):45^63, 1996. [59] R.N. Tchuraev and V.A. Ratner. A continuous approach with threshold characteristics for simulation of gene expression. In K. Bellmann, editor, Molecular Genetic Information Systems: Modelling and Simulation, pages 64^80. Akademie-Verlag, Berlin, 1983. [60] D. Thieffry, A.M. Huerta, E. Perez-Rueda, and J. Collado-Vides. From specific gene regulation to genomic networks: A global analysis of transcriptional regulation in Escherichia coli. BioEssays, 20:433^440, 1998. [61] D. Thieffry and R. Thomas. Dynamical behaviour of biological networks: II. Immunity control in bacteriophage lambda. Bulletin of Mathematical Biology, 57(2):277 297. 1995. [62] R. Thomas and R. d'Ari. Biological Feedback. CRC Press, Boca Raton, FL, 1990. [63] R. Thomas, D. Thieffry, and M. Kaufman. Dynamical behaviour of biological regulatory networks: I. Biological role of feedback loops and practical use of the concept of the loop-characteristic state. Bulletin of Mathematical Biology, 57(2):217 276. 1995. [64] R.B. Trelease, R.A. Henderson, and J.B. Park. A qualitative process system for modeling NF-kB and AP-1 gene regulation in immune cell biology research. Artifical Intelligence in Medicine, 17:303-321, 1999. [65] G. von Dassow, E. Meir, E.M. Munro, and G.M. Odell. The segment polarity network is a robust developmental module. Nature, 406:188^192, 2000. [66] D.M. Wolf and F.H. Eeckman. On the relationship between genomic regulatory element organization and gene regulatory dynamics. Journal of Theoretical Biology, 195:167^186, 1998. [67] G. Yagil and E. Yagil. On the relation between effector concentration and the rate of induced enzyme synthesis. Biophysical Journal, 11:11 27. 1971. [68] H. Zhu and M. Snyder. Protein arrays and microarrays. Current Opinion in Chemical Biology, 5:11) 15. 2001. RR n° 4407 Unite de recherche INRIA Rhone-Alpes 655, avenue de l'Europe - 38330 Montbonnot-St-Martin (France) Unite de recherche INRIA Lorraine : LORIA, Technopole de Nancy-Brabois - Campus scientifique 615, rue du Jardin Botanique - BP 101 - 54602 Villers-les-Nancy Cedex (France) Unite de recherche INRIA Rennes : IRISA, Campus universitaire de Beaulieu - 35042 Rennes Cedex (France) Unite de recherche INRIA Rocquencourt : Domaine de Voluceau - Rocquencourt - BP 105 - 78153 Le Chesnay Cedex (France) Unite de recherche INRIA Sophia Antipolis : 2004, route des Lucioles - BP 93 - 06902 Sophia Antipolis Cedex (France) Editeur INRIA - Domaine de Voluceau - Rocquencourt, BP 105 - 78153 Le Chesnay Cedex (France) http://www .i n ria .f r ISSN 0249-6399