A computational framework for qualitative simulation of nonlinear dynamical models of gene-regulatory networks

Background Due to the huge amount of information at genomic level made recently available by high-throughput experimental technologies, networks of regulatory interactions between genes and gene products, the so-called gene-regulatory networks, can be uncovered. Most networks of interest are quite intricate because of both the high dimension of interacting elements and the complexity of the kinds of interactions between them. Then, mathematical and computational modeling frameworks are a must to predict the network behavior in response to environmental stimuli. A specific class of Ordinary Differential Equations (ODE) has shown to be adequate to describe the essential features of the dynamics of gene-regulatory networks. But, deriving quantitative predictions of the network dynamics through the numerical simulation of such models is mostly impracticable as they are currently characterized by incomplete knowledge of biochemical reactions underlying regulatory interactions, and of numeric values of kinetic parameters. Results This paper presents a computational framework for qualitative simulation of a class of ODE models, based on the assumption that gene regulation is threshold-dependent, i.e. only effective above or below a certain threshold. The simulation algorithm we propose assumes that threshold-dependent regulation mechanisms are modeled by continuous steep sigmoid functions, unlike other simulation tools that considerably simplifies the problem by approximating threshold-regulated response functions by step functions discontinuous in the thresholds. The algorithm results from the interplay between methods to deal with incomplete knowledge and to study phenomena that occur at different time-scales. Conclusion The work herein presented establishes the computational groundwork for a sound and a complete algorithm capable to capture the dynamical properties that depend only on the network structure and are invariant for ranges of values of kinetic parameters. At the current state of knowledge, the exploitation of such a tool is rather appropriate and useful to understand how specific activity patterns derive from given network structures, and what different types of dynamical behaviors are possible.


Background
Although, up to now, there is no model that efficiently and accurately represents the gene interactions underlying regulatory mechanisms in their whole complexity, recent experimental evidence seems to confirm the adequacy of a specific class of ODEs to describe their essential dynamical features. These models assume that genes are controlled by transcriptor factors, and that the effect of a transcription factor on the transcription rate of a gene, the response function, is threshold-dependent. Such switch-like behaviors across variable thresholds are mathematically well represented by steep sigmoidal functions.
Although these models provide detailed description of gene regulatory molecular mechanisms [1][2][3], their predictive usefulness, at quantitative level, is rather limited even when the network at hand is very well studied. In fact, the exploitation of classical numerical approaches is mostly impracticable as precise and quantitative information on (i) the biochemical reaction mechanisms underlying regulatory interactions, and (ii) kinetic parameters and threshold concentrations are currently unknown and not identifiable from available data. However, at the current state of knowledge, qualitative predictions of the dynamical properties are not make-shift solutions but rather appropriate to get insight into the functioning of systems not completely understood as molecular interaction networks are.
The application of generic qualitative simulation approaches, as originally proposed in the Artificial Intelligence research framework [4] under the label QR, is not the proper solution. The mathematical tools they are grounded on are too simple to compensate for knowledge incompleteness. This results in a number of drawbacks, e.g. their inability to upscalability, the exponential growth of the generated behaviors, and the generation of spurious behaviors, that seriously limit the range of applicability of such methods to predict the nonlinear dynamics of regulatory networks. A qualitative study of GRNs dynamics could, in theory, be performed by analytical methods [5][6][7][8] based on the classical theory of qualitative analysis of dynamical systems, and properly adapted to the specific class of models. But, in practice, the network complexity makes quite hard or even unfeasible traditional analysis.
Pioneering work towards automated qualitative analysis and simulation of GRNs has resulted in a computational tool, called GNA [9]. Although GNA has been successfully applied to study real world networks, namely the initiation of sporulation in Bacillus subtilis [10], and the response to nutritional stress and carbon starvation in Escherichia coli [11,12], it is not flawless. GNA circumvents the hard problem of dealing with sigmoidal nonlinear response functions by approximating them with step functions, discontinuous in the threshold hyperplanes. Such an assumption considerably simplifies the analysis as the model results in piecewise-linear equations, but it raises the problem to find a proper continuous solution across the threshold hyperplanes, or, in other words, to seek for generalized solutions of ODEs with discontinuous right-side terms. The solution to this problem is not straightforward as (i) there exist in the literature several definitions of generalized solutions, (ii) it is not yet completely understood what are the relationships between different definitions, and then, (iii) it is not clear how to choose the "right" definition for a particular task [13]. GNA adopts the Filippov approach [14] that results particularly convenient to deal with control problems but it may present drawbacks when applied to approximate the limit solutions of a continuous ODE model: it might find "too many" solutions, and fail to reach all stable ones. Thus, GNA suffers from the same problems, that in addition to those raised by a further approximation introduced to deal with computational issues, might compromise its soundness and completeness (Dordan O, Ironi L, Panzeri L, Some critical remarks on GNA, in preparation).
The qualitative simulation algorithm we propose works under the assumptions that (i) threshold-dependent regulation mechanisms are modeled by continuous steep sigmoid functions, and (ii) any two genes are never regulated at the same threshold by a certain variable. The sigmoidal-nonlinearities make the problem quite hard to be tackled. But, the assumption that all sigmoids have very high steepness allows us to apply a systematic way of analysis. Let us observe that, due to the switch-like behavior of the response functions around the thresholds, the GRN dynamics occurs at different time-scales. To be able to deal with both slow and fast nonlinear dynamics we theoretically base our algorithm on a classical singular perturbation analysis method properly adapted to the assumed class of ODEs [8,15]. Such a method suitably combined with QR key concepts computationally drives, starting from initial conditions, the construction of all possible state transitions, and calculates the sets of symbolic inequalities on parameter values that hold when specific transitions occur.
A class of models of GRN dynamics As phenomenological model of the complex dynamics of GRNs made up of n components, we consider the following generic equations: where the dot denotes time derivative, x i is the concentration of the i-th gene product, g i > 0 is the decay rate of x i , Z is a vector with Z jk as components, and Z jk = S(x j , θ jk , q) is a sigmoid function with threshold θ jk . The response, or regulatory, function S : R + [0, 1] is a continuous monotonic S-shaped map depending on the parameter q (0 <q ≪ 1), that determines the steepness of S around the threshold value θ jk , such that for q 0 we have S(x j , θ jk , q) = 0 (respectively 1) when the value of x j is smaller (greater) than θ jk (Fig. 1).
Each x i , with domain Ω i ⊂ R + , is associated with n i thresholds ordered according to θ ij <θ ik if j <k. The state equations describe the balance between the production process f i (Z) and the degradation one, herein supposed to be linear. The functions f i are multilinear polynomials in the variables Z jk , and are frequently composed by algebraic equivalents of Boolean functions. More precisely, where il are real values that denote kinetic rate parameters, L i is the possibly empty number of interactions that synthesize x i , and, in accordance with the network structure, a jkl assumes value either equal to 1 when Z jk takes part in the l-th interaction or equal to 0 otherwise.
Let us represent the dynamics of the n state variables x i , modeled by the Eq. (1) and associated with appropriate initial conditions, in the phase space. The ordered sets Θ i of the n i threshold values θ ij associated with each x i naturally induce a partition of the phase space into qualitatively distinct domains. In the set Δ of all the domains identified by the partition, we can distinguish the set Δ s ⊂ Δ of switching domains from the set Δ r ⊂ Δ of regular domains, such that Δ = Δ s ∪ Δ r (Fig. 2).
A regular domain D r (also called a box) is an open rectangular region between adjacent threshold hyperplanes in which the values of all response functions Z jk Figure 1 Response function shape. The sigmoid response function S(x 1 , θ 1 , q) describes the relationship between the concentration of a gene product x 1 and the relative production rate of the regulated gene. The parameter q determines its steepness around the threshold θ 1 .

Figure 2
An example of phase space partition. Partition of the phase space associated to a dynamical system with two state variables, where each variable x i [0, x i ], i = 1, 2 is associated with two ordered thresholds θ ij , j = 1, 2. The white rectangles denote regular domains. The gray rectangles denote switching domains, whose width, δ > 0 is a monotonic function of the parameter q with δ(q) 0 for q 0. In the light gray regions one variable is switching and the other one is regular, whereas in the dark gray ones both variables are switching. BMC Bioinformatics 2009, 10(Suppl 12):S14 http://www.biomedcentral.com/1471-2105/10/S12/S14 equal either 0 or 1. A switching domain D s is a narrow region, whose width depends on the parameter q, surrounding either a section of a threshold hyperplane or an intersection of threshold hyperplanes. In a switching domain D s we distinguish s (D s ) switching variables, x s , from ns (D s ) regular ones x r . The values of a variables x s lies in the neighborhood of one of its thresholds and Z sk takes values in (0,1), while the values of a variable x r lies between two adjacent thresholds. The network dynamics in each domain D Δ is described by different models. The motion in D r Δ r is described by linear ODEs, and its analysis is straightforward. The motion in D s Δ s , or equivalently around the thresholds, is described by nonlinear equations, and occurs at different time-scales. Therefore, to capture the salient features of the nonlinear dynamics in a switching domain, and to determine how the trajectories cross it to move towards other domains, a specific mathematical method is required.

Methods
Singular perturbation analysis in the switching domains Singular Perturbation Analysis (SPA) is a classical approach to study phenomena that occur at different time-scales [33]. The dynamics of such phenomena are described by ODEs, associated with appropriate initial conditions, in which a small parameter (0 <q ≪ 1) multiplies either one of the derivatives or higher order derivative. Let us indicate this initial value problem by ℳ q , and by ℳ 0 the same problem where q = 0. As q 0, the solution of ℳ q identifies a "small" region, called boundary-layer region, of non-uniform convergence to the solution of the reduced system ℳ 0 . The region of uniformconvergence of ℳ q to ℳ 0 is called outer region. Taken together, the outer and boundary-layer solution approximate the solution of ℳ q for small nonzero values of q.
Such a general approach is not directly applicable to our problem, but it must be tailored to it [8]. In outline, (i) the Eqs. (1) related to the x s variables are rewritten into the standard form of SPA through a change of coordinate system, (ii) the boundary-layer and outer solutions are calculated in the new coordinates, and (iii) they are converted back into the usual frame of reference.
Let Σ: Ω ↦ [0, 1] n be the coordinate transformation that converts the x s coordinates into the Z s ones. As under our assumptions, where Z S , Z R are the vectors of switching and regular variables Z s and Z r , respectively.
The fast dynamics in the boundary-layer is studied by suitably scaling the time variable, namely τ = t/q. In the limit, the fast dynamics is obtained by solving the system: associated with appropriate initial conditions, where the prime denotes the derivative with respect to τ. As Z R is constant in any D s Δ s , we focus on the switching variables Z s only. The system (4) has a manifold of stationary points, called slow-manifold, given by the solutions, for all s, of the stationary equations ′ Z s = 0. We call exit point set (EP) the set of points in the slowmanifold that are stable and satisfy the Tikhonov-Levinson theorem, and we call Z-cube Ƶ(D s ) = [ , ] ( ) 0 1 s D s the frame of reference where we search for exit points.
The motion in the original time t along the exit points is described by the reduced system ℳ 0 . Thus, under the hypothesis that at least one exit point % Z S exists, the motion equations of regular variables, or equivalently the outer solution, is represented by: The problem (5) is linear, and then, given the initial conditions, the outer solution, that determines how the trajectories move along the x r directions, is easily calculated.
The calculation of the slow-manifold generally leads to an heavy, nonlinear computational problem as it consists in the solution of a set of polynomial equations.
In the current implementation we adopt a further assumption that sounds realistic: Assumption A . Every gene product only regulates one gene at each of its thresholds Such an assumption considerably simplifies the calculation of the slow-manifold. Mathematically, it implies that each Z jk only occurs in one equation, and thus the terms f s (Z R , Z S ) are linear.

Remark 1
The location of each exit point is crucial in our analysis as it indicates the next adjacent domains the trajectories are moving towards along the x s directions.
BMC Bioinformatics 2009, 10(Suppl 12):S14 http://www.biomedcentral.com/1471-2105/10/S12/S14 Let us observe that stationary points always exist on the vertices of Ƶ(D s ) as they are the roots of d s (Z s , θ s ) = 0. The computational cost of the search for all the other exit points could be quite high, but it can be considerably reduced by checking first a necessary condition for the existence of a stationary point on the other elements of Ƶ (D s ). Let F be a face or the interior of Ƶ(D s ). In [15], it has been proved that necessary condition for the existence of a stationary point in F is that the Jacobian matrix ⎟ restricted to the switching variables in F has a complete loop. This holds if and only if there is a nonzero loop involving all variables in J F , and can be checked by using concepts from graph theory.
The stability of % Z is checked by analyzing the spectrum of the Jacobian matrix, and the condition (ii) is verified when, Remark 2 SPA works out in the limit q 0, but the calculated solution approximates the solution of Eq. (3) for sufficiently small q (0 <q ≪ 1). Moreover, it can be proved that the Jacobian matrix is stable for 0 <q < q ≪ 1. Thus, the exit points calculated in the limit are also valid for values of q < q .

Remark 3
Under the Assumption A , the reduced equations are always independent of the variables Z s in the Eq. (4), and the two sets of equations are mutually independent. Thus, the behavior of the variables x s in a Z-cube is completely independent of the values of variables x r , and the motion in a switching domain may be studied by first analyzing the former variables, and then the latter ones.

Qualitative reasoning concepts
Among the generic qualitative approaches proposed in the literature, QSIM provides both the most suitable formalism and algorithm to represent and simulate models qualitatively abstracted from ODEs [4]. Thus, we revise and ad hoc tailor its key concepts to our specific class of models.

Qualitative value
The real values of each state variable x i with domain Ω i = [0, x i ] are discretized into a finite ordered set of values qualitatively important from the biological point of view. In our context, the qualitative value space of each x i is defined by the ordered set Θ i made up of both its n i threshold symbolic values and the endpoints of Ω i . The partition of the whole system domain, induced by the sets θ i , i = 1,..., n identifies qualitatively distinct hyperrectangles D that define all possible system qualitative values.

Qualitative state
Let A(D) be the set of domains adjacent to D Δ. The qualitative state of D, QS(D), is defined by all of its adjacent domains D k towards which a transition from it is possible: Each transition from D identifies a domain next traversed by a system trajectory. More precisely, if we number by i the domain D traversed at time t i , each D k QS(D) will be traversed by different trajectories at time t i+1 .

State transition
Qualitative simulation of network dynamics is achieved by iteratively applying local transition strategies from one domain to its adjacent domains. The possible transitions from any D are determined by different strategies according to whether D Δ r or D Δ s .
In the case D Δ r , like in traditional QR methods and in GNA [9], possible transitions are determined by the signs of & x i . As & x i are defined by linear expressions, such signs are easily determined by exploiting the inequalities that define the parameter space domain.
In the case D Δ s , a sign-based strategy is not practicable as the expressions for & x i are nonlinear. A convenient way to proceed is given by SPA: transitions from D towards adjacent D k are determined by the locations of the exit points on the boundary elements of the associated Ƶ(D). Due to Assumption A , each boundary element of Ƶ(D) may contain at most one exit point. But, as, in a qualitative context, knowledge incompleteness on the parameter values is expressed by coarse constraints, the boundary element of Ƶ(D) where an exit point is located is not, in general, uniquely determined. Thus, unless the exit point is located in the interior of Ƶ(D) and a stable solution is reached, the successors of D are not uniquely determined. However, through symbolic computation procedures, it is possible to calculate the set of inequalities, I k i , on parameter values that hold when the transition from D i to D k occurs. Then, the 3-tuple 〈D i , D k , I k i 〉 clearly and uniquely identifies each path from D i to D k .

Qualitative behavior
A finite sequence of paths, where each path is clearly both linked and consistent with its predecessor and successor, defines a qualitative behavior: , where D 0 is the initial domain, and D F either contains a stable fixed point or identifies a cycle, i.e it is an already visited domain. I 0 is the initial set of inequalities that defines the parameter space domain, and I F the set of inequalities on parameter values associated with D F .

Qualitative simulation
Starting from an initial domain D 0 and a set, I 0 , of symbolic inequalities on parameter values, qualitative simulation generates all possible state transitions, and represents them by a directed tree rooted in D 0 , BT(D 0 ), where the vertices correspond to D i , and the arcs, labeled by the inequalities I j i , to the transitions from D i to D j . Each branch in BT(D 0 ) defines a qualitative trajectory QB from D 0 , that traverses specific domains and occurs when the values of parameters satisfy its related inequalities. Such a trajectory abstracts all those numerical solution of the ODE model with initial conditions x 0 varying in D 0 , and with kinetic parameter values fulfilling the inequalities.

Results
Given as input, (i) n symbolic state equations of the form (1); (ii) n qualitative value spaces Θ i = {θ ij } of the state variables; (iii) D 0 Δ; (iv) a set of symbolic inequalities I 0 on parameters defining a parameter space domain PSD 0 , the algorithm, through an iterative procedure initialized by 〈D 0 , I 0 〉, provides as output the entire range of possible network dynamics represented by BT(D 0 ). Its main steps are outlined in the following: 1. Define the qualitative values by partitioning the phase space into regular and switching domains. Step 2, and in particular the calculation of the conditions on parameters I k i , is the core of the overall algorithm. A transition from D i to D k actually occurs, and then D k QS(D i ), only if the set I k i is consistent with the set I 0 , i.e. when it defines a not empty parameter space domain PSD k i such that PSD k i ⊆ PSD 0 . The calculation of an inequality set I k i is performed by two distinct algorithms that implement a different transition strategy according to whether D i Δ r or D i Δ s .

Moving from a regular domain
The algorithm in charge of the construction of the possible paths from regular domains is, in principle, similar to that one proposed by GNA, but it is more informative as it calculates the I k i s. In the limit, we indicate D i Δ r by the (a) Symbolic state equations in D i Δ r . In each box D i , Z jk equals either 0 or 1 in the step function limit. This simplifies Eq. (1) as they reduces to linear equations: where μ j depends on D i , and is given by the sum of some jl s. From Eq. (6) we can easily find the focal point if v j = 1 or (c) Possible transitions from D i to D k . If the calculated inequality set defines a not empty parameter space domain BMC Bioinformatics 2009, 10(Suppl 12):S14 http://www.biomedcentral.com/1471-2105/10/S12/S14 PSD k i ⊆ PSD 0 then a transition towards D k is possible and the qualitative state QS(D i ) is updated accordingly.
To exemplify how the algorithm works we consider the ODE model: where the parameters are positive, and all the response functions are expressed by the Hill function commonly used in the literature.
Let us consider the domain D 11 in Fig. 2 as current D i , and let us define I 0 as follows: The conditions on parameters in Table 1 Table 1, I 12 11 is the only one in agreement with I 0 . Thus, QS(D 11 ) = {D 12 } and the only possible transition from D 11 occurs when I 12 11 ∧ I 0 holds.

Moving from a switching domain
In a domain D i Δ s , the nonlinear dynamics is characterized by fast and slow motions, respectively associated with x s and x r , that are independently calculated. Let us reindex the variables x j , Z j such that the switching variables come first, and proceed first with the construction of the fast motion by exploiting the SPA approach.

A -Fast motion
The study of the fast dynamics is performed in Ƶ(D i ) in the scaled time τ = t/q, and aims at localizing the set of exit points in Ƶ(D i ) rather than at detailing the dynamics within it. Such points clearly identify the next adjacent domains the trajectories are moving towards from D i along the x s directions.
To this end, the algorithm defines, in the limit, the mapping Σ D i : D i Ƶ(D i ), where D i = D i ∪ A(D i ), such that the interior of Ƶ(D i ), and its boundary are the images of D i , and A(D i ), respectively. More precisely, the domains D k A(D i ) are mapped into the faces of Ƶ(D i ) when D k Δ s or into its vertices, otherwise. Let ℱ be the set of both the faces and the interior of Ƶ(D i ), and F its generic element.
To exemplify, let us consider the switching domain D 7 in Fig. 2. The set ℱ has five elements, the four faces corresponding to D 2 , D 6 , D 8 , D 12 , and the interior of Ƶ (D 7 ) that correspond to D 7 . Moreover, the vertices of Ƶ (D 7 ) are the images, through the mapping Σ, of the adjacent regular domains D 1 , D 11 , D 3 , and D 13 (Fig. 3).
(a) 1. Symbolic state equations in D i OE Δ s The algorithm symbolically calculates the boundarylayer equations (4) in the Z variables associated with the current switching domain.
As an example, let us consider the model (7), and the domain D 7 , where both variables x 1 , x 2 are switching in it, respectively, around the thresholds θ 11 , θ 21 . In D 7 , characterized by fast motion only, the boundary-layer system is given by: ) ( ) . q k g q (9)  . Finally, the exit point candidate set is updated with the points % Z 2 and % Z 12 (Fig. 3(a)). 2. imposes conditions I k, l on the sign of ′ Z l , ∀ l L F , in a neighborhood of % Z k so that the motion of Z l heads towards % Z k ; 3. for those % Z k located on elements of ℱ, imposes the further condition, I k,(0, 1) , 0 < % Z s k < 1 for each % Z s k that does not take value 0 or 1. holds.
Conditions 2. and 3. are easily checked, while the stability condition is checked by using concepts from graph theory, and the usual definition of stability based on the sign of the eigenvalues of J F . It follows from Assumption A that the matrix J F is block-structured, where each block is a permutation matrix associated with a sub-loop. It can be proved that % Z k is stable if: (i) J F k has no blocks with dimension strictly greater than 2; (ii) in 1-dimensional blocks the elements are negative; (iii) in 2-dimensional blocks the loop products are positive.

Remark 4
Let us consider the general case when D i Δ s is characterized by both switching and regular variables. The motion from an exit point located in the interior of D i occurs in a sliding mode along a stable point in the slow-manifold of the boundary-layer system, and is described, in the normal time, by the reduced system. In the example, both J F 2 and J F 12 are 1-dimensional block matrices with negative elements. Then, both exit points % Z 2 and % Z 12 fulfill the stability condition 1. without further constraints imposed on the parameters. The condition 2. on variable Z l , l = 2 imposes: The inequality I 12,2 defined by (11) is compatible with (8), but the inequality I 2,2 is not. Then, % Z 2 is removed from the exit point set. To be an exit point % Z 12 must satisfy the condition 3.: Finally, % Z 12 is an exit point if I 12 7 , defined by I 0 ∧ I 12,2 ∧ I 12,(0,1) , holds.
As for vertices in the example, the conditions 1. and 2. are fulfilled in the point % Z 11 = (0, 1) that corresponds to the vertex defined as image of D 11 by the map Σ.

B -Slow motion
The slow motion of regular variables x r along the exit points is studied in the normal time t in the usual frame of reference, and it is reconstructed from the reduced system (5) through the same symbolic procedure given for regular domains.
Qualitative simulation outcomes Through the described iterative procedure, the algorithm builds (i) a behavior tree, rooted in the initial domain, and (ii) sets of inequalities on parameter values that characterize specific state transitions. The behavior tree captures the entire spectrum of network dynamical behaviors that are invariant for ranges of parameter values defined by the calculated inequalities.
To illustrate the type of output the algorithm produces, let us consider the results of the simulation of the ODE model (7) starting from D 1 with the parameter space defined by the inequalities I 0 (8). Through the described iterative procedure, the algorithm builds the behavior tree, BT(D 1 ), showed in Fig. 4, and calculates the inequalities on the parameters, listed in Table 2, that are associated with each path in BT(D 1 ). As in the example n = 2, the trajectories described by the tree are also easily summarized in the phase plane (Fig. 5). Three reachable stable states, located in D 11 , D 12 and D 5 , are identified by the final leaf of each branch in BT. As D 12 Δ s , one of them is a singular stable point whereas the others are regular stable points. These stable states are reached by different predicted qualitative behaviors, each of them occurring under specific constraints on parameters. For example, the trajectory QB 13 starting from D 1 , crossing D 6 , and reaching a stable point in D 11 is allowed when the inequalities I I , and I 11 , consistent with I 0 , hold.
In general, for n > 2, due to graphical difficulties, each branch in the tree is given a time representation. More precisely, a specific qualitative behavior, e.g. QB 4 , is described by the temporal evolution of each variable (Fig.  6). The behavior in Fig. 6 abstracts all numerical trajectories of the ODE model (7) obtained with different initial values in D 1 , and different sets of parameter values that fulfill the same inequalities associated with QB 4 (Fig. 7).

Discussion
The algorithm is currently under implementation. As for symbolic calculus, the implementation requires to tackle complex tasks, such as: (i) update an inequality set with Behavior tree. Behavior tree rooted in D 1 , BT(D 1 ), obtained by the simulation of the model Eqs. (7) with initial condition 〈D 1 , I 0 〉, where I 0 is defined by Eq. (8). Each branch in the tree defines a qualitative trajectory that occurs when the inequalities that characterize each path in the branch hold.  BMC Bioinformatics 2009, 10(Suppl 12):S14 http://www.biomedcentral.com/1471-2105/10/S12/S14 an another one; (ii) check the consistency of two sets of inequalities I 1 and I 2 ; (iii) solve systems of equations; (iv) find cycles in the Jacobian matrix. As for (iii), the original equations are multilinear in Z s , but due to Assumption A they can be straightforward solved and analyzed for stability. Also the solution of problems (i) and (ii) benefits from Assumption A as the inequalities are always linear. Then, thanks to the Assumption A , and to algorithms proposed both by the literature and common symbolic computation package, such as Mathematica [34], the tasks (i)-(iii) are simplified and feasible. As for the task (iv), it is performed by using cycle-detection algorithms and matrix graph theory tools [35].

Soundness and completeness Soundness
The algorithm guarantees that the behavior tree captures all of the sound behaviors for values of q sufficiently small. A closed-form expression of a symbolic upper bound of q, q , in terms of model parameters guarantees the Jacobian matrix stability. Then, the generated behaviors are sound for values of q < q . However, for a complete formal proof of soundness we need to prove that assuming, in the limit, stability instead of asymptotic stability does not affect the results for values of q < q .

Figure 5
Representation of the qualitative trajectories in the phase space. Phase space representation of trajectories defined by BT(D 1 ) in Fig. 4, filtered of the spurious behaviors QB 2 and QB 11 . Stable states are denoted by ·. The labels I k i denote the sets of inequalities on parameter values that hold when transitions from D i to D k occur.

Completeness
At the current stage, the algorithm may generate spurious behaviors. There is a twofold explanation for that. First, we have not yet performed a thorough analysis with respect to entrance-exit transition, or in other words, we have not yet solved the problem of identifying the only admissible connections between entrance and exit points. Moreover, singular perturbation analysis is a local procedure that works quite well in a quantitative context but that needs, in a qualitative context, to be supported by a global criterion when local paths are connected to produce a specific trajectory. For example, the behavior QB 2 is spurious as I 12 7 is not consistent with I 11 12 . Similarly, QB 11 is spurious.
The characterization of the paths from one domain to the next ones by sets of inequalities constraining the model parameters is quite new in the field of qualitative simulation, as for both general-purpose and specifically tailored algorithms. Such a strategy may reveal quite useful in the definition of a "global criterion" that allows us to distinguish sound behaviors from spurious ones by requiring that the sets of inequalities that label the local paths in a specific trajectory are consistent with each other. Both the definition of such a criterion and its implementation are not a trivial task, especially from a computational point of view. As for algorithm completeness, another essential methodological and computational issue to be deepened deals with the definition of the transition map that properly connects the entrance points to the exit points associated with a switching domain.
Comparison between the algorithm and GNA Until now, to the best of our knowledge, GNA is the only computational counterpart of the algorithm herein proposed. Unlike GNA, that simplifies the simulation problem by approximating the sigmoid response functions by step functions discontinuous in the thresholds, our qualitative simulation algorithm tackles the problem in all of its complexity, and works for models of GNAs with continuous sigmoid response functions. As a matter of fact, it has been designed to both provide a tool for a more realistic modeling framework and overcome the limits of GNA. Furthermore, in addition to more solid mathematical grounds for soundness and completeness, our algorithm differs from GNA as far as the required input information and the output outcomes are concerned. In GNA, the trajectories are constructed under the assumption that the locations of the focal points associated with the regular domains, symbolically expressed by a set of inequalities on parameters, are specified as input. But, as this precise knowledge is unlikely available, a thorough study of the dynamics of the system at hand might require several runs of the algorithm to simulate the GNA dynamics under possible different assumptions on focal point locations. Then, the study could result quite heavy due to the high number of possible scenarios, and the need to check, to avoid further causes of generation of spurious behaviors, the consistency of all the parameter inequalities that precisely localize a specific focal point.
Our algorithm does not need precise information on the location of focal points but, to be fired, it just requires the parameter space domain, that is inequalities on parameters that univocally determine the sign of direction of change of each state variable in the initial domain. Afterwards, it is up to the algorithm to calculate and refine the parameter inequalities that distinguish each simulated behavior. Due to more relaxed constraints, a single run of our algorithm leads to a richer and more informative set of behaviors than that one obtained by GNA as it is highlighted by the comparison of Fig. 5 with Fig. 8. Let us observe that any other possible definition of parameter inequalities would not have led to the set of behaviors in Fig. 5 but to one of its subsets.
A complete and fair performance evaluation of our approach when compared with GNA would not disregard its application to those real world networks successfully simulated within the GNA framework, namely the initiation of sporulation in Bacillus subtilis [10], and the response to nutritional stress and carbon starvation in Escherichia coli [11,12].

Application to Systems/Synthetic Biology
The qualitative simulation algorithm presented provides powerful computational grounds for dry experiments to be performed in both Systems and Synthetic Biology contexts. It is an economic tool for hypothesis testing and knowledge discovery as it allows us to understand how specific activation patterns are derived from models of fixed network structures, and which different dynamical behaviors are possible.
In a Systems Biology context, its exploitation is useful in the formulation phase of new hypotheses and theories that explain, in both physiological and pathological situations, experimental data either previously unobserved or not interpretable by the well-established biological knowledge. The result of the comparison of the simulated behaviors with the observed ones allows us either to confirm or to refuse the model that represents the new formulated biological hypotheses. Moreover, as qualitative simulation derives the entire range of network dynamics consistent with the initial conditions, the simulation outcomes might capture behaviors never observed. Such behaviors need to be experimentally tested by ad hoc designed experiments to possibly further confirm or suggest how to revise the model, i.e. the underlying biological hypotheses. The effective applicative potential of the proposed algorithm within a model-based reasoning cycle for the deep comprehension of complex real world biological systems is currently under study. More precisely, it will be applied (1) to explain and interpret the dynamics of the complex regulatory network that controls motility in Bacillus subtilis [36], and (2) to explore the molecular mechanisms that regulate the transition of a normal cell phenotype to an invasive-metastatic phenotype [37].
The model-based reasoning cycle results particularly convenient and useful also when applied to Synthetic Biology problems in the design phase. The construction of a synthetic biological network in the cell, able to perform specific tasks or to produce desired behaviors of a biological process in response to external stimuli, could benefit from a preliminary dry design of the network. To this end, models of different network structures can be simulated in correspondence to different inputs, i.e signals either exogenous or cellular, till the desired stable behaviors are obtained. Let us remind that, in our simulation framework, each specific simulated behavior is characterized by a chain of parameter inequalities. Thus, the actual construction of a synthetic network could benefit from such a piece of information as it provides the admissible parameter space domain for the target behaviors. The synthetic network IRMA [38] built in the yeast Saccharomyces cerevisiae is made up of a small number of components but rather complex in their interconnections as it includes multiple feedback loops generated by the combination of transcriptional activators and repressors. It offers a suitable benchmark for in vivo testing our computational framework.

Conclusion
The assumption of continuous response functions makes the simulation problem hard to be tackled but it is crucial in view of the realization of tools that can be gradually extended to tackle more and more realistic models. Our algorithm is grounded on a set of symbolic computation algorithms that carry out the integration of qualitative reasoning techniques with singular analysis perturbation methods: the former techniques allow us to cope with uncertain and incomplete knowledge whereas the latter ones lay the mathematical groundwork for a sound and complete algorithm capable to deal with regulation processes that occur at different time-scales.
The modeling framework and the simulation algorithm proposed can be applied to predict the possible qualitative behaviors of GRNs of any size and complexity, both natural and synthetic. The range of applicability of such a computational approach is quite large. First, it makes possible the validation of hypothesized models of GRNs by matching the simulated predictions against observed gene expression profiles, the derivation of the most plausible model, and the identification of parameter inequalities that account for the observations. Then, it greatly facilitates the investigation of the effects on the dynamics of a specific GRN in response to external stimuli. Finally, it allows us to build a total envisionment of the model through the generation of all possible state transitions, and to search, within them, for the parameter constraints and initial conditions that allow us to achieve a desired behavior, e.g. a stable solution. In term of tasks, besides its undoubted contribution to the comprehension of complex networks of genes and interactions, the proposed computational approach could be fruitfully used in the design of either new drugs or synthetic regulatory networks.