 Research
 Open Access
 Published:
A computational framework for qualitative simulation of nonlinear dynamical models of generegulatory networks
BMC Bioinformatics volume 10, Article number: S14 (2009)
Abstract
Background
Due to the huge amount of information at genomic level made recently available by highthroughput experimental technologies, networks of regulatory interactions between genes and gene products, the socalled generegulatory 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 generegulatory 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 thresholddependent, i.e. only effective above or below a certain threshold. The simulation algorithm we propose assumes that thresholddependent regulation mechanisms are modeled by continuous steep sigmoid functions, unlike other simulation tools that considerably simplifies the problem by approximating thresholdregulated 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 timescales.
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 ODE s 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 thresholddependent. Such switchlike behaviors across variable thresholds are mathematically well represented by steep sigmoidal functions.
Although these models provide detailed description of gene regulatory molecular mechanisms [1–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 makeshift 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 GRN s dynamics could, in theory, be performed by analytical methods [5–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 GRN s 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 piecewiselinear 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 ODE s with discontinuous rightside 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) thresholddependent 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 sigmoidalnonlinearities 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 switchlike behavior of the response functions around the thresholds, the GRN dynamics occurs at different timescales. 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 ODE s [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 GRN s made up of n components, we consider the following generic equations:
where the dot denotes time derivative, x_{ i }is the concentration of the ith gene product, γ_{ 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 Sshaped 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, α_{ jkl }assumes value either equal to 1 when Z_{ jk }takes part in the lth interaction or equal to 0 otherwise.
The validity of the biological assumptions underlying the model described by Eq. (1) has been confirmed by recent experimental evidence [16–28] and theoretical studies [29–32]. Thus, we can reasonably assume that such a model is suitable to give a phenomenological description of a wide range of regulatory systems in which the combined effects of a series of genetic processes, e.g. transcriptional and translational regulation, proteinprotein interactions, metabolic processes, etc., can be properly described by thresholddependent response functions.
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 }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 σ (D_{ s }) switching variables, x_{ s }, from n  σ (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 ODE s, 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 timescales. 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 timescales [33]. The dynamics of such phenomena are described by ODE s, 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 boundarylayer region, of nonuniform 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 boundarylayer 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 boundarylayer 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 d_{ s }is a continuous and limited function, we can rewrite Eq. (1) as in the following:
where Z_{ S }, Z_{ R }are the vectors of switching and regular variables Z_{ s }and Z_{ r }, respectively.
The fast dynamics in the boundarylayer 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 slowmanifold, given by the solutions, for all s, of the stationary equations = 0. We call exit point set (EP) the set of points in the slowmanifold that are stable and satisfy the TikhonovLevinson theorem, and we call Zcube Ƶ(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 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 slowmanifold 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 . Every gene product only regulates one gene at each of its thresholds
Such an assumption considerably simplifies the calculation of the slowmanifold. 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.
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.
Let be a stationary point located on ∈ Ƶ(D_{ s }), and = {l : l ∈ {l,..., σ (D_{ s })}, ∈ {0, 1}}. is an exit point if (i) it is stable. If is on the boundary of Ƶ(D_{ s }), it is required, in addition to (i), that (ii) Z_{ l }, ∀ l ∈ , heads towards . The stability of is checked by analyzing the spectrum of the Jacobian matrix, and the condition (ii) is verified when, ∀ l ∈ , the sign of , in a neighborhood of is in agreement with the motion of Z_{ l }towards . More precisely, in a neighborhood of , > 0 and = 1 or < 0 and = 0.
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 < ≪ 1. Thus, the exit points calculated in the limit are also valid for values of q <.
Remark 3
Under the Assumption , 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 Zcube 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 ODE s [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, ] 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 . As 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 signbased strategy is not practicable as the expressions for 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 , 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, , on parameter values that hold when the transition from D_{ i }to D_{ k }occurs. Then, the 3tuple ⟨D_{ i }, D_{ k }, ⟩ 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 , 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.

2.
Calculate the qualitative state QS(D_{ i }) of the current domain D_{ i }.

(a)
Calculate the symbolic state equations in D_{ i }, and the set A(D_{ i });

(b)
Determine constraints on parameters that need to be fulfilled to have a transition from D_{ i }to D_{ k }∈ A(D_{ i });

(c)
Check consistency of with I_{0}, and build the path e_{ ik }= D_{ i }→ D_{ k };

3.
Append ⟨D_{ i }, D_{ k }, ⟩ to BT(D_{0}), and mark D_{ i }as visited domain.

4.
Repeat from step 2 for each not visited D_{ k }.
Step 2, and in particular the calculation of the conditions on parameters , 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 is consistent with the set I_{0}, i.e. when it defines a not empty parameter space domain such that ⊆ PSD_{0}. The calculation of an inequality set 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 s. In the limit, we indicate D_{ i }∈ Δ_{ r }by the product where (θ_{ ji }, θ_{j(i+1)}) denotes the interval of x_{ j }in D_{ i }.

(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:
(6)
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 the trajectories are heading towards. Herein, we assume that focal points do not belong to switching domains. If x* belongs to D_{ i }, there is a stable point in it. Such a stable point exists in D_{ i }, i.e. D_{ i }∈ QS(D_{ i }), if the set of inequalities ∀j ∈ {1,..., n} exists and is consistent with I_{0}. Otherwise, the trajectories move towards a switching domain D_{ k }adjacent to D_{ i }.

(b)
Constraints on parameters. ∀ D_{ k }∈ A(D_{ i }), the algorithm calculates the set of inequalities on parameters that need to be fulfilled to have a transition from D_{ i }to D_{ k }. As in D_{ i }all the equations (6) are linear, and all the trajectories head towards a focal point x* in a regular domain, such inequalities are calculated by imposing that the signs of state variable rates match the relative position of D_{ k }with respect to D_{ i }. We define the relative position of D_{ k }with respect to D_{ i }, indicated by where v_{ j }∈ {1, 0, 1}, through the comparison of the intervals defining D_{ i }and D_{ k }. , initialized to I_{0}, is updated, ∀ j ∈ {1,..., n}, with either the inequality if v_{ j }= 1 or if v_{ j }= 1.

(c)
Possible transitions from D_{ i }to D_{ k }. If the calculated inequality set defines a not empty parameter space domain ⊆ 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 that are consistent with I_{0} identify the possible transitions from D_{11} towards adjacent domains in A(D_{11}) = {D_{6}, D_{7}, D_{12}, D_{16}, D_{17}}. Among the inequalities in Table 1, 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_{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 }), where = 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 }∈ Δ_{ 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 boundarylayer system is given by:
2. Identification of the candidate next traversed domains D_{ k }
Among all the domains D_{ k }∈ A(D_{ i }) only those that correspond to an element on the boundary of Ƶ(D_{ i }) where a stationary point is located are possible domains next traversed by the trajectories. Then, the algorithm first builds the subset of elements of A(D_{ i }) that are candidate successors of D_{ i }.
Let us denote by EP the set of stationary points, initially made up of the vertices of Ƶ(D_{ i }). Under the Assumption A, each element F ∈ ℱ, may contain at most one stationary point. A necessary condition for the existence of a stationary point on F is the presence of a nonzero loop in J_{ F }, the Jacobian matrix calculated on F. Then, ∀ F ∈ ℱ, the algorithm calculates J_{ F }and searches for a nonzero loop involving all variables in J_{ F }. In case, it symbolically calculates the stationary point on F, and updates accordingly the set of candidate exit points EP. Let us observe that one internal stationary point exists if all variables in Ƶ(D_{ i }) are involved in a loop.
In the example, only and have a nonzero loop. Then, the algorithm looks for the stationary state on F_{2} and F_{12}: and . Finally, the exit point candidate set is updated with the points and (Fig. 3(a)).
(b) Constraints on parameters
The inequality set , initialized to I_{0}, is calculated for each candidate exit point by requiring that each point is stable and fulfills other motion conditions on parameters. More precisely, the algorithm:

1.
analyzes the spectrum of the Jacobian matrix J_{ F }, and imposes possible conditions on parameters to ensure stability of ;

2.
imposes conditions I_{k, l}on the sign of , ∀ l ∈ ℒ_{ F }, in a neighborhood of so that the motion of Z_{ l }heads towards ;

3.
for those located on elements of ℱ, imposes the further condition, I_{k,(0, 1)}, 0 < < 1 for each that does not take value 0 or 1.
Finally, is an exit point if 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 that the matrix J_{ F }is blockstructured, where each block is a permutation matrix associated with a subloop. It can be proved that is stable if: (i) has no blocks with dimension strictly greater than 2; (ii) in 1dimensional blocks the elements are negative; (iii) in 2dimensional 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 slowmanifold of the boundarylayer system, and is described, in the normal time, by the reduced system. Then, a stable stationary point exists in D_{ i }, i.e. D_{ i }∈ QS(D_{ i }), if a stable point exists in the interior of Ƶ(D_{ i }), and the regular variable rates are zero in a point inside D_{ i }. The set of inequalities that checks the latter condition are defined by ∀j ∈ {σ (D_{ i }) + 1,..., n}.
(c) Possible transition from D_{ i }to D_{ k }
The exit points located on Ƶ(D_{ i }) clearly identify the set of all possible exit domains, i.e. those domains towards which a transition from D_{ i }is possible. Such domains are easily calculated by applying the map Σ^{1} to each element of Ƶ(D_{ i }) that contains an exit point. Let us observe that the remaining unstable stationary points in EP are possible entrance points to D_{ i }.
In the example, both and are 1dimensional block matrices with negative elements. Then, both exit points and 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, is removed from the exit point set. To be an exit point must satisfy the condition 3.:
Finally, is an exit point if , 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 = (0, 1) that corresponds to the vertex defined as image of D_{11} by the map Σ.
Finally, the only exit domains are D_{12}, D_{11}, and then QS(D_{7}) = {D_{12}, D_{11}} (Fig. 3(b)).
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 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 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 they can be straightforward solved and analyzed for stability. Also the solution of problems (i) and (ii) benefits from Assumption as the inequalities are always linear. Then, thanks to the Assumption , 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 cycledetection 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 closedform expression of a symbolic upper bound of q, , in terms of model parameters guarantees the Jacobian matrix stability. Then, the generated behaviors are sound for values of 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 <.
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 entranceexit 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 is not consistent with . 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 generalpurpose 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 GNA s 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 wellestablished 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 modelbased 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 invasivemetastatic phenotype [37].
The modelbased 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 timescales.
The modeling framework and the simulation algorithm proposed can be applied to predict the possible qualitative behaviors of GRN s 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 GRN s 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.
The code will be freely available for nonprofit academic research by making a user licence request to ironi@imati.cnr.it.
References
 1.
Glass L, Kauffman SA: The logical analysis of continuous, nonlinear biochemical control networks. Journal of Theoretical Biology 1973, 39: 103–129. 10.1016/00225193(73)902087
 2.
Plahte E, Mestl T, Omholt SW: A methodological basis for description and analysis of systems with complex switchlike interactions. Journal of Mathematical Biology 1998, 36(4):321–348. 10.1007/s002850050103
 3.
de Jong H: Modeling and Simulation of Genetic Regulatory Systems: A Literature Review. Journal of Computational Biology 2002, 9: 67–103. 10.1089/10665270252833208
 4.
Kuipers B: Qualitative Reasoning: Modeling and Simulation with Incomplete Knowledge. Cambridge, MA: MIT Press; 1994.
 5.
Glass L: Global analysis of nonlinear chemical kinetics. In Statistical Mechanics, Part B: Time Dependent Processes. Edited by: Berne B. Plenum Press, New York; 1977:311–349.
 6.
Hasty J, McMillen D, Isaacs F, Collins JJ: Computational studies of gene regulatory networks: In numero molecular biology. Nat Rev Genet 2001, 2: 268–279. 10.1038/35066056
 7.
Gouzè JL, Sari T: A class of piecewise linear differential equations arising in biological models. Dynamical systems 2003, 17: 299–316.
 8.
Plahte E, Kjøglum S: Analysis and generic properties of gene regulatory networks with graded response functions. Physica D: Nonlinear Phenomena 2005, 201(1–2):150–176. 10.1016/j.physd.2004.11.014
 9.
de Jong H, Gouzé JL, Hernandez C, Page M, Sari T, Geiselmann J: Qualitative Simulations of Genetic Regulatory Networks Using Piecewise Linear Models. Bulletin of Mathematical Biology 2004, 66(2):301–340. 10.1016/j.bulm.2003.08.010
 10.
de Jong H, Geiselmann J, Batt G, Hernandez C, Page M: Qualitative simulation of the initiation of sporulation in Bacillus subtilis . Bulletin of Mathematical Biology 2004, 66(2):261–300. 10.1016/j.bulm.2003.08.009
 11.
Ropers D, de Jong H, Geiselmann J: Mathematical modeling of genetic regulatory networks: Stress response in Escherichia coli . In Systems and Synthetic Biology. Edited by: Fu P, Latterich M, Panke S. Wiley & Sons; in press.
 12.
Ropers D, de Jong H, Page M, Schneider D, Geiselmann J: Qualitative simulation of the carbon starvation response in Escherichia coli . BioSystems 2006, 84(2):124–152. 10.1016/j.biosystems.2005.10.005
 13.
Bacciotti A: Some remarks on generalized solutions of discontinuous differential equations. Int Journal of Pure and Applied Mathematics 2003, 10(3):257–266.
 14.
Filippov AF: Differential equations with discontinuous right hand sides. Dordrecht: Kluwer Academic Publishers Group; 1988.
 15.
Veflingstad SR, Plahte E: Analysis of gene regulatory network models with graded and binary transcriptional responses. Biosystems 2007, 90: 323–339. 10.1016/j.biosystems.2006.09.036
 16.
Goldbeter A, Koshland DE: An Amplified Sensitivity Arising from Covalent Modification in Biological Systems. PNAS 1981, 78(11):6840–6844. 10.1073/pnas.78.11.6840
 17.
Yuh CH, Bolouri H, Davidson EH: Genomic cis Regulatory Logic: Experimental and Computational Analysis of a Sea Urchin Gene. Science 1998, 279(5358):1896–1902. 10.1126/science.279.5358.1896
 18.
Gardner T, Cantor C, Collins J: Construction of a genetic toggle switch in Escherichia coli . Nature 2000, 403(6767):339–342. 10.1038/35002131
 19.
Rossi FMV, Kringstein AM, Spicher A, Guicherit OM, Blau HM: Transcriptional control: Rheostat converted to on/off switch. Molecular Cell 2000, 6(3):723–728. 10.1016/S10972765(00)000708
 20.
Biggar SR, Crabtree GR: Cell signaling can direct either binary or graded transcriptional responses. The EMBO Journal 2001, 20: 3167–3176. 10.1093/emboj/20.12.3167
 21.
Yuh CH, Bolouri H, Davidson EH: cis regulatory logic in the endo16 gene: switching from a specification to a differentiation mode of control. Development 2001, 128(5):617–629.
 22.
Guet CC, Elowitz MB, Hsing WH, Leibler S: Combinatorial synthesis of genetic networks. Science 2002, 296(5572):1466–1470. 10.1126/science.1067407
 23.
Setty Y, Mayo AE, Surette MG, Alon U: Detailed map of a cis regulatory input function. PNAS 2003, 100(13):7702–7707. 10.1073/pnas.1230759100
 24.
Kramer BP, Fischer C, Fussenegger M: BioLogic gates enable logical transcription control in mammalian cells. Biotechnol Bioeng 2004, 87(4):478–484. 10.1002/bit.20142
 25.
Istrail S, Davidson EH: Logic functions of the genomic cis regulatory code. PNAS 2005, 102(14):4954–4959. 10.1073/pnas.0409624102
 26.
Rosenfeld N, Young JW, Alon U, Swain PS, Elowitz MB: Gene Regulation at the SingleCell Level. Science 2005, 307(5717):1962–1965. 10.1126/science.1106914
 27.
Kim J, Bates DG, Postlethwaite I, Ma L, Iglesias PA: Robustness analysis of biochemical network models. Syst Biol (Stevange) 2006, 153(3):96–104. 10.1049/ipsyb:20050024
 28.
Mayo AE, Setty Y, Shavit S, Zaslaver A, Alon U: Plasticity of the cis Regulatory Input Function of a Gene. PLoS Biology 2006, 4(4):e45. 10.1371/journal.pbio.0040045
 29.
Wolf DM, Eeckman FH: On the Relationship Between Genomic Regulatory Element Organization and Gene Regulatory Dynamics. Journal of Theoretical Biology 1998, 195(2):167–186. 10.1006/jtbi.1998.0790
 30.
Buchler NE, Gerland U, Hwa T: Nonlinear protein degradation and the function of genetic circuits. PNAS 2005, 102(27):9559–9564. 10.1073/pnas.0409553102
 31.
Bintu L, Buchler NE, Garcia HG, Gerland U, Hwa T, Kondev J, Phillips R: Transcriptional regulation by the numbers: models. Current Opinion in Genetics & Development 2005, 15(2):116–124. 10.1016/j.gde.2005.02.007
 32.
Bintu L, Buchler NE, Garcia HG, Gerland U, Hwa T, Kondev J, Kuhlman T, Phillips R: Transcriptional regulation by the numbers: applications. Current Opinion in Genetics & Development 2005, 15(2):125–135. 10.1016/j.gde.2005.02.006
 33.
Holmes M: Introduction to Perturbation Methods. Berlin: Springer; 1995.
 34.
Wolfram S: The Mathematica Book. Wolfram Media; 2003.
 35.
Gross J, Yellen J: Graph Theory and its Applications. New York: Chapman & Hall/CRC Press; 2006.
 36.
Calvio C, Osera C, Amati G, Galizzi A: Autoregulation of swrAA and motility in Bacillus subtilis . J Bacteriol 2008, 190: 5720–5728. 10.1128/JB.0045508
 37.
Gentile A, D'Alessandro L, Lazzari L, Martinoglio B, Bertotti A, Mira A, Lanzetti L, Comoglio PM, Medico E: Metdriven invasive growth involves transcriptional regulation of Arhgap12. Oncogene 2008, 27: 5590–5598. 10.1038/onc.2008.173
 38.
Cantone I, Marucci L, Iorio F, Ricci M, Belcastro V, Bansal M, Santini S, di Bernardo M, di Bernardo D, Cosma M: A Yeast Synthetic Network for In Vivo Assessment of ReverseEngineering and Modeling Approaches. Cell 2009, 137: 172–181. 10.1016/j.cell.2009.01.055
Acknowledgements
This research is carried out within the Interdepartmental CNRBIOINFORMATICS Project. We gratefully acknowledge O. Dordan, E. Plahte, and V. Simoncini for the useful discussions on different methodological aspects and mathematical problems related to the presented work.
This article has been published as part of BMC Bioinformatics Volume 10 Supplement 12, 2009: Bioinformatics Methods for Biomedical Complex System Applications. The full contents of the supplement are available online at http://www.biomedcentral.com/14712105/10?issue=S12.
Author information
Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
LI devised the work, and drafted the manuscript. Both authors participated in the design and development of the algorithms. LP implemented the software. All authors read and approved the final manuscript.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
Ironi, L., Panzeri, L. A computational framework for qualitative simulation of nonlinear dynamical models of generegulatory networks. BMC Bioinformatics 10, S14 (2009). https://doi.org/10.1186/1471210510S12S14
Published:
DOI: https://doi.org/10.1186/1471210510S12S14
Keywords
 Ordinary Differential Equation
 Exit Point
 Switching Domain
 Regular Domain
 Adjacent Domain