 Research
 Open Access
A computational framework for qualitative simulation of nonlinear dynamical models of generegulatory networks
 Liliana Ironi^{1}Email author and
 Luigi Panzeri^{1}
https://doi.org/10.1186/1471210510S12S14
© Ironi and Panzeri; licensee BioMed Central Ltd. 2009
 Published: 15 October 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.
Keywords
 Ordinary Differential Equation
 Exit Point
 Switching Domain
 Regular Domain
 Adjacent Domain
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
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.
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.
where Z_{ S }, Z_{ R }are the vectors of switching and regular variables Z_{ s }and Z_{ r }, respectively.
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 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
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
 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
 (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:
 (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.
where the parameters are positive, and all the response functions are expressed by the Hill function commonly used in the literature.
Example parameter inequalities
Inequality Constraints  

 Sign Constraints  Parameter Inequalities 
 < 0 

 > 0 

 > 0, < 0 

 > 0 

 > 0, > 0 

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.
(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.
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
 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 }.
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.
Inequalities associated with BT(D_{1})
, I_{12}  I_{0} ∧ 

 I_{0} ∧ 
 I_{0} ∧ 
 I_{0} ∧ 
 I_{0} ∧ 
I _{5}  I_{0} ∧ 
 I _{0} 
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.
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.
Declarations
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.
Authors’ Affiliations
References
 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)902087View ArticlePubMedGoogle Scholar
 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/s002850050103View ArticlePubMedGoogle Scholar
 de Jong H: Modeling and Simulation of Genetic Regulatory Systems: A Literature Review. Journal of Computational Biology 2002, 9: 67–103. 10.1089/10665270252833208View ArticlePubMedGoogle Scholar
 Kuipers B: Qualitative Reasoning: Modeling and Simulation with Incomplete Knowledge. Cambridge, MA: MIT Press; 1994.Google Scholar
 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.View ArticleGoogle Scholar
 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/35066056View ArticlePubMedGoogle Scholar
 Gouzè JL, Sari T: A class of piecewise linear differential equations arising in biological models. Dynamical systems 2003, 17: 299–316.View ArticleGoogle Scholar
 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.014View ArticleGoogle Scholar
 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.010View ArticlePubMedGoogle Scholar
 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.009View ArticlePubMedGoogle Scholar
 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.Google Scholar
 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.005View ArticlePubMedGoogle Scholar
 Bacciotti A: Some remarks on generalized solutions of discontinuous differential equations. Int Journal of Pure and Applied Mathematics 2003, 10(3):257–266.Google Scholar
 Filippov AF: Differential equations with discontinuous right hand sides. Dordrecht: Kluwer Academic Publishers Group; 1988.View ArticleGoogle Scholar
 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.036View ArticlePubMedGoogle Scholar
 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.6840PubMed CentralView ArticlePubMedGoogle Scholar
 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.1896View ArticlePubMedGoogle Scholar
 Gardner T, Cantor C, Collins J: Construction of a genetic toggle switch in Escherichia coli . Nature 2000, 403(6767):339–342. 10.1038/35002131View ArticlePubMedGoogle Scholar
 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)000708View ArticlePubMedGoogle Scholar
 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.3167PubMed CentralView ArticlePubMedGoogle Scholar
 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.PubMedGoogle Scholar
 Guet CC, Elowitz MB, Hsing WH, Leibler S: Combinatorial synthesis of genetic networks. Science 2002, 296(5572):1466–1470. 10.1126/science.1067407View ArticlePubMedGoogle Scholar
 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.1230759100PubMed CentralView ArticlePubMedGoogle Scholar
 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.20142View ArticlePubMedGoogle Scholar
 Istrail S, Davidson EH: Logic functions of the genomic cis regulatory code. PNAS 2005, 102(14):4954–4959. 10.1073/pnas.0409624102PubMed CentralView ArticlePubMedGoogle Scholar
 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.1106914View ArticlePubMedGoogle Scholar
 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:20050024View ArticleGoogle Scholar
 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.0040045PubMed CentralView ArticlePubMedGoogle Scholar
 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.0790View ArticlePubMedGoogle Scholar
 Buchler NE, Gerland U, Hwa T: Nonlinear protein degradation and the function of genetic circuits. PNAS 2005, 102(27):9559–9564. 10.1073/pnas.0409553102PubMed CentralView ArticlePubMedGoogle Scholar
 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.007View ArticleGoogle Scholar
 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.006View ArticleGoogle Scholar
 Holmes M: Introduction to Perturbation Methods. Berlin: Springer; 1995.View ArticleGoogle Scholar
 Wolfram S: The Mathematica Book. Wolfram Media; 2003.Google Scholar
 Gross J, Yellen J: Graph Theory and its Applications. New York: Chapman & Hall/CRC Press; 2006.Google Scholar
 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.0045508PubMed CentralView ArticlePubMedGoogle Scholar
 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.173View ArticlePubMedGoogle Scholar
 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.055View ArticlePubMedGoogle Scholar
Copyright
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.