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

Abstract

Background

Due to the huge amount of information at genomic level made recently available by high-throughput experimental technologies, networks of regulatory interactions between genes and gene products, the so-called gene-regulatory networks, can be uncovered. Most networks of interest are quite intricate because of both the high dimension of interacting elements and the complexity of the kinds of interactions between them. Then, mathematical and computational modeling frameworks are a must to predict the network behavior in response to environmental stimuli. A specific class of Ordinary Differential Equations (ODE) has shown to be adequate to describe the essential features of the dynamics of gene-regulatory networks. But, deriving quantitative predictions of the network dynamics through the numerical simulation of such models is mostly impracticable as they are currently characterized by incomplete knowledge of biochemical reactions underlying regulatory interactions, and of numeric values of kinetic parameters.

Results

This paper presents a computational framework for qualitative simulation of a class of ODE models, based on the assumption that gene regulation is threshold-dependent, i.e. only effective above or below a certain threshold. The simulation algorithm we propose assumes that threshold-dependent regulation mechanisms are modeled by continuous steep sigmoid functions, unlike other simulation tools that considerably simplifies the problem by approximating threshold-regulated response functions by step functions discontinuous in the thresholds. The algorithm results from the interplay between methods to deal with incomplete knowledge and to study phenomena that occur at different time-scales.

Conclusion

The work herein presented establishes the computational groundwork for a sound and a complete algorithm capable to capture the dynamical properties that depend only on the network structure and are invariant for ranges of values of kinetic parameters. At the current state of knowledge, the exploitation of such a tool is rather appropriate and useful to understand how specific activity patterns derive from given network structures, and what different types of dynamical behaviors are possible.

Background

Although, up to now, there is no model that efficiently and accurately represents the gene interactions underlying regulatory mechanisms in their whole complexity, recent experimental evidence seems to confirm the adequacy of a specific class of 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 threshold-dependent. Such switch-like behaviors across variable thresholds are mathematically well represented by steep sigmoidal functions.

Although these models provide detailed description of gene regulatory molecular mechanisms [1–3], their predictive usefulness, at quantitative level, is rather limited even when the network at hand is very well studied. In fact, the exploitation of classical numerical approaches is mostly impracticable as precise and quantitative information on (i) the biochemical reaction mechanisms underlying regulatory interactions, and (ii) kinetic parameters and threshold concentrations are currently unknown and not identifiable from available data. However, at the current state of knowledge, qualitative predictions of the dynamical properties are not make-shift solutions but rather appropriate to get insight into the functioning of systems not completely understood as molecular interaction networks are.

The application of generic qualitative simulation approaches, as originally proposed in the Artificial Intelligence research framework [4] under the label QR, is not the proper solution. The mathematical tools they are grounded on are too simple to compensate for knowledge incompleteness. This results in a number of drawbacks, e.g. their inability to upscalability, the exponential growth of the generated behaviors, and the generation of spurious behaviors, that seriously limit the range of applicability of such methods to predict the nonlinear dynamics of regulatory networks. A qualitative study of 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 piecewise-linear equations, but it raises the problem to find a proper continuous solution across the threshold hyperplanes, or, in other words, to seek for generalized solutions of ODE s with discontinuous right-side terms. The solution to this problem is not straightforward as (i) there exist in the literature several definitions of generalized solutions, (ii) it is not yet completely understood what are the relationships between different definitions, and then, (iii) it is not clear how to choose the "right" definition for a particular task [13]. GNA adopts the Filippov approach [14] that results particularly convenient to deal with control problems but it may present drawbacks when applied to approximate the limit solutions of a continuous ODE model: it might find "too many" solutions, and fail to reach all stable ones. Thus, GNA suffers from the same problems, that in addition to those raised by a further approximation introduced to deal with computational issues, might compromise its soundness and completeness (Dordan O, Ironi L, Panzeri L, Some critical remarks on GNA, in preparation).

The qualitative simulation algorithm we propose works under the assumptions that (i) threshold-dependent regulation mechanisms are modeled by continuous steep sigmoid functions, and (ii) any two genes are never regulated at the same threshold by a certain variable. The sigmoidal-nonlinearities make the problem quite hard to be tackled. But, the assumption that all sigmoids have very high steepness allows us to apply a systematic way of analysis. Let us observe that, due to the switch-like behavior of the response functions around the thresholds, the GRN dynamics occurs at different time-scales. To be able to deal with both slow and fast nonlinear dynamics we theoretically base our algorithm on a classical singular perturbation analysis method properly adapted to the assumed class of 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:

(1)

where the dot denotes time derivative, x _i is the concentration of the i-th 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 S-shaped map depending on the parameter q (0 <q ≪ 1), that determines the steepness of S around the threshold value θ _jk , such that for q → 0 we have S(x _j , θ _jk , q) = 0 (respectively 1) when the value of x _j is smaller (greater) than θ _jk (Fig. 1).

Figure 1

Response function shape. The sigmoid response function S(x₁, θ₁, q) describes the relationship between the concentration of a gene product x₁ and the relative production rate of the regulated gene. The parameter q determines its steepness around the threshold θ₁.

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,

(2)

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 l-th 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, protein-protein interactions, metabolic processes, etc., can be properly described by threshold-dependent 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).

Figure 2

An example of phase space partition. Partition of the phase space associated to a dynamical system with two state variables, where each variable x _i ∈ [0, ], i = 1, 2 is associated with two ordered thresholds θ _ij , j = 1, 2. The white rectangles denote regular domains. The gray rectangles denote switching domains, whose width, δ > 0 is a monotonic function of the parameter q with δ(q) → 0 for q → 0. In the light gray regions one variable is switching and the other one is regular, whereas in the dark gray ones both variables are switching.

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 time-scales. Therefore, to capture the salient features of the nonlinear dynamics in a switching domain, and to determine how the trajectories cross it to move towards other domains, a specific mathematical method is required.

Methods

Singular perturbation analysis in the switching domains

Singular Perturbation Analysis (SPA) is a classical approach to study phenomena that occur at different time-scales [33]. The dynamics of such phenomena are described by 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 ℳ₀ the same problem where q = 0. As q → 0, the solution of ℳ _q identifies a "small" region, called boundary-layer region, of non-uniform convergence to the solution of the reduced system ℳ₀. The region of uniform-convergence of ℳ _q to ℳ₀ is called outer region. Taken together, the outer and boundary-layer solution approximate the solution of ℳ _q for small nonzero values of q.

Such a general approach is not directly applicable to our problem, but it must be tailored to it [8]. In outline, (i) the Eqs. (1) related to the x _s variables are rewritten into the standard form of SPA through a change of coordinate system, (ii) the boundary-layer and outer solutions are calculated in the new coordinates, and (iii) they are converted back into the usual frame of reference.

Let Σ: Ω ↦ [0, 1]ⁿ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:

(3)

where Z _S , Z _R are the vectors of switching and regular variables Z _s and Z _r , respectively.

The fast dynamics in the boundary-layer is studied by suitably scaling the time variable, namely τ = t/q. In the limit, the fast dynamics is obtained by solving the system:

(4)

associated with appropriate initial conditions, where the prime denotes the derivative with respect to τ. As Z _R is constant in any D _s ∈ Δ _s , we focus on the switching variables Z _s only. The system (4) has a manifold of stationary points, called slow-manifold, given by the solutions, for all s, of the stationary equations = 0. We call exit point set (EP) the set of points in the slow-manifold that are stable and satisfy the Tikhonov-Levinson theorem, and we call Z-cube Ƶ(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 ℳ₀. 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:

(5)

The problem (5) is linear, and then, given the initial conditions, the outer solution, that determines how the trajectories move along the x _r directions, is easily calculated.

The calculation of the slow-manifold generally leads to an heavy, nonlinear computational problem as it consists in the solution of a set of polynomial equations. In the current implementation we adopt a further assumption that sounds realistic:

Assumption . Every gene product only regulates one gene at each of its thresholds

Such an assumption considerably simplifies the calculation of the slow-manifold. Mathematically, it implies that each Z _jk only occurs in one equation, and thus the terms f _s (Z _R , Z _S ) are linear.

Remark 1

The location of each exit point is crucial in our analysis as it indicates the next adjacent domains the trajectories are moving towards along the x _s directions.

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 non-zero 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 Z-cube is completely independent of the values of variables x _r , and the motion in a switching domain may be studied by first analyzing the former variables, and then the latter ones.

Qualitative reasoning concepts

Among the generic qualitative approaches proposed in the literature, QSIM provides both the most suitable formalism and algorithm to represent and simulate models qualitatively abstracted from 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 hyper-rectangles 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 sign-based 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 3-tuple ⟨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₀ 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₀ 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₀ and a set, I₀, of symbolic inequalities on parameter values, qualitative simulation generates all possible state transitions, and represents them by a directed tree rooted in D₀, BT(D₀), 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₀) defines a qualitative trajectory QB from D₀, 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⁰ varying in D₀, 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₀ ∈ Δ; (iv) a set of symbolic inequalities I₀ on parameters defining a parameter space domain PSD₀, the algorithm, through an iterative procedure initialized by ⟨D₀, I₀⟩, provides as output the entire range of possible network dynamics represented by BT(D₀). 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₀, and build the path e _ik = D _i → D _k ;

• 3. Append ⟨D _i , D _k , ⟩ to BT(D₀), 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₀, i.e. when it defines a not empty parameter space domain such that ⊆ PSD₀. 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₀. 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₀, 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₀ 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:

(7)

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₁₁ in Fig. 2 as current D _i , and let us define I₀ as follows:

(8)

The conditions on parameters in Table 1 that are consistent with I₀ identify the possible transitions from D₁₁ towards adjacent domains in A(D₁₁) = {D₆, D₇, D₁₂, D₁₆, D₁₇}. Among the inequalities in Table 1, is the only one in agreement with I₀. Thus, QS(D₁₁) = {D₁₂} and the only possible transition from D₁₁ occurs when ∧ I₀ holds.

Table 1 Example parameter inequalities

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₇ in Fig. 2. The set ℱ has five elements, the four faces corresponding to D₂, D₆, D₈, D₁₂, and the interior of Ƶ(D₇) that correspond to D₇. Moreover, the vertices of Ƶ(D₇) are the images, through the mapping Σ, of the adjacent regular domains D₁, D₁₁, D₃, and D₁₃ (Fig. 3).

Figure 3

Mapping between D _i and Ƶ ( D _i ). (a) Mapping of D₇ and its adjacent domains in Fig. 2 into the elements of Ƶ(D₇). The empty circles on the boundary of Ƶ(D₇) denote the candidate exit points. (b) The filled circles denote the exit points from Ƶ(D₇), whereas the empty circles denote unstable stationary points that correspond to possible entrance points to Ƶ(D₇). The single-headed arrows denote the direction of change of Z _l , ∀l ∈ ℒ _F in a neighborhood of the exit points.

(a) 1. Symbolic state equations in D _i ∈ Δ _s

The algorithm symbolically calculates the boundary-layer 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₇, where both variables x₁, x₂ are switching in it, respectively, around the thresholds θ₁₁, θ₂₁. In D₇, characterized by fast motion only, the boundary-layer system is given by:

(9)

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 non-zero loop in J _F , the Jacobian matrix calculated on F. Then, ∀ F ∈ ℱ, the algorithm calculates J _F and searches for a non-zero 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 non-zero loop. Then, the algorithm looks for the stationary state on F₂ and F₁₂: 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₀, 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 block-structured, where each block is a permutation matrix associated with a sub-loop. It can be proved that is stable if: (i) has no blocks with dimension strictly greater than 2; (ii) in 1-dimensional blocks the elements are negative; (iii) in 2-dimensional blocks the loop products are positive.

Remark 4

Let us consider the general case when D _i ∈ Δ _s is characterized by both switching and regular variables. The motion from an exit point located in the interior of D _i occurs in a sliding mode along a stable point in the slow-manifold of the boundary-layer system, and is described, in the normal time, by the reduced system. 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 1-dimensional 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:

(10)

(11)

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₀ ∧ 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₁₁ by the map Σ.

Finally, the only exit domains are D₁₂, D₁₁, and then QS(D₇) = {D₁₂, D₁₁} (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₁ with the parameter space defined by the inequalities I₀ (8). Through the described iterative procedure, the algorithm builds the behavior tree, BT(D₁), showed in Fig. 4, and calculates the inequalities on the parameters, listed in Table 2, that are associated with each path in BT(D₁). 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₁₁, D₁₂ and D₅, are identified by the final leaf of each branch in BT. As D₁₂ ∈ Δ _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₁₃ starting from D₁, crossing D₆, and reaching a stable point in D₁₁ is allowed when the inequalities and I₁₁, consistent with I₀, hold.

Table 2 Inequalities associated with BT(D₁)

Figure 4

Behavior tree. Behavior tree rooted in D₁, BT(D₁), obtained by the simulation of the model Eqs.(7) with initial condition ⟨D₁, I₀⟩, where I₀ is defined by Eq. (8). Each branch in the tree defines a qualitative trajectory that occurs when the inequalities that characterize each path in the branch hold.

Figure 5

Representation of the qualitative trajectories in the phase space. Phase space representation of trajectories defined by BT(D₁) in Fig. 4, filtered of the spurious behaviors QB₂ and QB₁₁. Stable states are denoted by ·. The labels denote the sets of inequalities on parameter values that hold when transitions from D _i to D _k occur.

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₄, 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₁, and different sets of parameter values that fulfill the same inequalities associated with QB₄ (Fig. 7).

Figure 6

Temporal representation of a qualitative trajectory. Qualitative temporal representation of the behavior labeled QB₄ in Fig. 4, namely of the variables x₁ (a) and x₂ (b). On the horizontal axis, the discrete sequence of symbolic instants at which the domains are sequentially traversed. On the vertical axis, the qualitative value space of each variable is reported. The qualitative values (∘) abstract all real values in the related open and closed intervals, whose width is highlighed by a dashed line. The stable states are denoted by ·. The continuous line connecting the qualitative values has been introduced to give a hint of the qualitative profile.

Figure 7

Numerical simulation results. Time courses of the variables x₁ (a) and x₂ (b) obtained by numerical simulations of the model (7). The different plots result from simulations with different initial numeric values in D₁, and different sets of parameter numeric values but fulfilling the same inequalities that characterize the qualitative trajectory QB₄. The horizontal continuous lines in correspondence to threshold values should facilitate the comparison with Fig. 6. All plots related to each variable are qualitatively abstracted by the related plot in Fig. 6.

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₁ and I₂; (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 cycle-detection algorithms and matrix graph theory tools [35].

Soundness and completeness

Soundness

The algorithm guarantees that the behavior tree captures all of the sound behaviors for values of q sufficiently small. A closed-form expression of a symbolic upper bound of q, , 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 entrance-exit transition, or in other words, we have not yet solved the problem of identifying the only admissible connections between entrance and exit points. Moreover, singular perturbation analysis is a local procedure that works quite well in a quantitative context but that needs, in a qualitative context, to be supported by a global criterion when local paths are connected to produce a specific trajectory. For example, the behavior QB₂ is spurious as is not consistent with . Similarly, QB₁₁ is spurious.

The characterization of the paths from one domain to the next ones by sets of inequalities constraining the model parameters is quite new in the field of qualitative simulation, as for both general-purpose and specifically tailored algorithms. Such a strategy may reveal quite useful in the definition of a "global criterion" that allows us to distinguish sound behaviors from spurious ones by requiring that the sets of inequalities that label the local paths in a specific trajectory are consistent with each other. Both the definition of such a criterion and its implementation are not a trivial task, especially from a computational point of view. As for algorithm completeness, another essential methodological and computational issue to be deepened deals with the definition of the transition map that properly connects the entrance points to the exit points associated with a switching domain.

Comparison between the algorithm and GNA

Until now, to the best of our knowledge, GNA is the only computational counterpart of the algorithm herein proposed. Unlike GNA, that simplifies the simulation problem by approximating the sigmoid response functions by step functions discontinuous in the thresholds, our qualitative simulation algorithm tackles the problem in all of its complexity, and works for models of 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.

Figure 8

GNA simulation outcomes. Simulation outcomes of the example model obtained by the algorithm GNA given as input D₁ and the parameter inequalities, compatible with I₀, ; ; ; (panel (a)), ; ; ; (panel (b)). The states generated by GNA, and labeled S _n , have been located at their position in the phase space to make the comparison with the results in Fig. 5 possible. Both sets of trajectories in (a) and (b) are subsets of the set of trajectories generated by our algorithm.

A complete and fair performance evaluation of our approach when compared with GNA would not disregard its application to those real world networks successfully simulated within the GNA framework, namely the initiation of sporulation in Bacillus subtilis [10], and the response to nutritional stress and carbon starvation in Escherichia coli [11, 12].

Application to Systems/Synthetic Biology

The qualitative simulation algorithm presented provides powerful computational grounds for dry experiments to be performed in both Systems and Synthetic Biology contexts. It is an economic tool for hypothesis testing and knowledge discovery as it allows us to understand how specific activation patterns are derived from models of fixed network structures, and which different dynamical behaviors are possible.

In a Systems Biology context, its exploitation is useful in the formulation phase of new hypotheses and theories that explain, in both physiological and pathological situations, experimental data either previously unobserved or not interpretable by the well-established biological knowledge. The result of the comparison of the simulated behaviors with the observed ones allows us either to confirm or to refuse the model that represents the new formulated biological hypotheses. Moreover, as qualitative simulation derives the entire range of network dynamics consistent with the initial conditions, the simulation outcomes might capture behaviors never observed. Such behaviors need to be experimentally tested by ad hoc designed experiments to possibly further confirm or suggest how to revise the model, i.e. the underlying biological hypotheses. The effective applicative potential of the proposed algorithm within a model-based reasoning cycle for the deep comprehension of complex real world biological systems is currently under study. More precisely, it will be applied (1) to explain and interpret the dynamics of the complex regulatory network that controls motility in Bacillus subtilis [36], and (2) to explore the molecular mechanisms that regulate the transition of a normal cell phenotype to an invasive-metastatic phenotype [37].

The model-based reasoning cycle results particularly convenient and useful also when applied to Synthetic Biology problems in the design phase. The construction of a synthetic biological network in the cell, able to perform specific tasks or to produce desired behaviors of a biological process in response to external stimuli, could benefit from a preliminary dry design of the network. To this end, models of different network structures can be simulated in correspondence to different inputs, i.e signals either exogenous or cellular, till the desired stable behaviors are obtained. Let us remind that, in our simulation framework, each specific simulated behavior is characterized by a chain of parameter inequalities. Thus, the actual construction of a synthetic network could benefit from such a piece of information as it provides the admissible parameter space domain for the target behaviors. The synthetic network IRMA[38] built in the yeast Saccharomyces cerevisiae is made up of a small number of components but rather complex in their interconnections as it includes multiple feedback loops generated by the combination of transcriptional activators and repressors. It offers a suitable benchmark for in vivo testing our computational framework.

Conclusion

The assumption of continuous response functions makes the simulation problem hard to be tackled but it is crucial in view of the realization of tools that can be gradually extended to tackle more and more realistic models. Our algorithm is grounded on a set of symbolic computation algorithms that carry out the integration of qualitative reasoning techniques with singular analysis perturbation methods: the former techniques allow us to cope with uncertain and incomplete knowledge whereas the latter ones lay the mathematical groundwork for a sound and complete algorithm capable to deal with regulation processes that occur at different time-scales.

The modeling framework and the simulation algorithm proposed can be applied to predict the possible qualitative behaviors of 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 non-profit academic research by making a user licence request to ironi@imati.cnr.it.