Volume 10 Supplement 12

## Bioinformatics Methods for Biomedical Complex System Applications

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

- Liliana Ironi
^{1}Email author and - Luigi Panzeri
^{1}

**10(Suppl 12)**:S14

https://doi.org/10.1186/1471-2105-10-S12-S14

© 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 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.

## Keywords

## 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

*n*components, we consider the following generic equations:

*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).

*x*

_{ i }, with domain Ω

_{ i }⊂ R

^{+}, is associated with

*n*

_{ i }thresholds ordered according to

*θ*

_{ ij }<

*θ*

_{ ik }if

*j*<

*k*. The state equations describe the balance between the production process

*f*

_{ i }(

**Z**) and the degradation one, herein supposed to be linear. The functions

*f*

_{ i }are multilinear polynomials in the variables

*Z*

_{ jk }, and are frequently composed by algebraic equivalents of Boolean functions. More precisely,

where *κ*_{
il
}are real values that denote kinetic rate parameters, *L*_{
i
}is the possibly empty number of interactions that synthesize *x*_{
i
}, and, in accordance with the network structure, *α*_{
jkl
}assumes value either equal to 1 when *Z*_{
jk
}takes part in the *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.

*n*state variables

*x*

_{ i }, modeled by the Eq. (1) and associated with appropriate initial conditions, in the phase space. The ordered sets Θ

_{ i }of the

*n*

_{ i }threshold values

*θ*

_{ ij }associated with each

*x*

_{ i }naturally induce a partition of the phase space into qualitatively distinct domains. In the set Δ of all the domains identified by the partition, we can distinguish the set Δ

_{ s }⊂ Δ of

*switching domains*from the set Δ

_{ r }⊂ Δ of

*regular domains*, such that Δ = Δ

_{ s }∪ Δ

_{ r }(Fig. 2).

A regular domain *D*_{
r
}(also called a *box*) is an open rectangular region between adjacent threshold hyperplanes in which the values of all response functions *Z*_{
jk
}equal either 0 or 1. A switching domain *D*_{
s
}is a narrow region, whose width depends on the parameter *q*, surrounding either a section of a threshold hyperplane or an intersection of threshold hyperplanes. In a switching domain *D*_{
s
}we distinguish *σ* (*D*_{
s
}) switching variables, *x*_{
s
}, from *n* - *σ* (*D*_{
s
}) regular ones *x*_{
r
}. The values of a variables *x*_{
s
}lies in the neighborhood of one of its thresholds and *Z*_{
sk
}takes values in (0,1), while the values of a variable *x*_{
r
}lies between two adjacent thresholds. The network dynamics in each domain *D* ∈ Δ is described by different models. The motion in *D*_{
r
}∈ Δ_{
r
}is described by linear ODE s, and its analysis is straightforward. The motion in *D*_{
s
}∈ Δ_{
s
}, or equivalently around the thresholds, is described by nonlinear equations, and occurs at different 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 ℳ_{0} the same problem where *q* = 0. As *q* → 0, the solution of ℳ_{
q
}identifies a "small" region, called *boundary*-*layer region*, of non-uniform convergence to the solution of the *reduced system* ℳ_{0}. The region of uniform-convergence of ℳ_{
q
}to ℳ_{0} is called *outer region*. Taken together, the outer and boundary-layer solution approximate the solution of ℳ_{
q
}for small nonzero values of *q*.

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

^{ n }be the coordinate transformation that converts the

*x*

_{ s }coordinates into the

*Z*

_{ s }ones. As under our assumptions, , where

*d*

_{ s }is a continuous and limited function, we can rewrite Eq. (1) as in the following:

where **Z**_{
S
}, **Z**_{
R
}are the vectors of switching and regular variables *Z*_{
s
}and *Z*_{
r
}, respectively.

*τ*=

*t*/

*q*. In the limit, the fast dynamics is obtained by solving the system:

associated with appropriate initial conditions, where the prime denotes the derivative with respect to *τ*. As **Z**_{
R
} is constant in any *D*_{
s
}∈ Δ_{
s
}, we focus on the switching variables *Z*_{
s
}only. The system (4) has a manifold of stationary points, called slow-manifold, given by the solutions, for all *s*, of the stationary equations
= 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.

*t*along the exit points is described by the reduced system ℳ

_{0}. Thus, under the hypothesis that at least one exit point exists, the motion equations of regular variables, or equivalently the outer solution, is represented by:

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

The calculation of the 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

*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*_{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

*n*symbolic state equations of the form (1); (ii)

*n*qualitative value spaces Θ

_{ i }= {

*θ*

_{ ij }} of the state variables; (iii)

*D*

_{0}∈ Δ; (iv) a set of symbolic inequalities

*I*

_{0}on parameters defining a parameter space domain

*PSD*

_{0}, the algorithm, through an iterative procedure initialized by ⟨

*D*

_{0},

*I*

_{0}⟩, provides as output the entire range of possible network dynamics represented by BT(

*D*

_{0}). Its main steps are outlined in the following:

- 1.
*Define*the qualitative values by partitioning the phase space into regular and switching domains. - 2.
*Calculate*the qualitative state*QS*(*D*_{ i }) of the current domain*D*_{ i }. - (a)
Calculate the

*symbolic state equations in D*_{ i }, and the*set A*(*D*_{ i }); - (b)
Determine

*constraints**on parameters*that need to be fulfilled to have a transition from*D*_{ i }to*D*_{ k }∈*A*(*D*_{ i }); - (c)
- 3.
- 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

*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:

*μ*

_{ j }depends on

*D*

_{ i }, and is given by the sum of some

*κ*

_{ jl }s. From Eq. (6) we can easily find the

*focal point*the trajectories are heading towards. Herein, we assume that focal points do not belong to switching domains. If

**x*** belongs to

*D*

_{ i }, there is a stable point in it. Such a stable point exists in

*D*

_{ i }, i.e.

*D*

_{ i }∈

*QS*(

*D*

_{ i }), if the set of inequalities ∀

*j*∈ {1,...,

*n*} exists and is consistent with

*I*

_{0}. Otherwise, the trajectories move towards a switching domain

*D*

_{ k }adjacent to

*D*

_{ i }.

- (b)
*Constraints**on parameters*. ∀*D*_{ k }∈*A*(*D*_{ i }), the algorithm calculates the set of inequalities on parameters that need to be fulfilled to have a transition from*D*_{ i }to*D*_{ k }. As in*D*_{ i }all the equations (6) are linear, and all the trajectories head towards a focal point**x***in a regular domain, such inequalities are calculated by imposing that the signs of state variable rates match the relative position of*D*_{ k }with respect to*D*_{ i }. We define the relative position of*D*_{ k }with respect to*D*_{ i }, indicated by where*v*_{ j }∈ {-1, 0, 1}, through the comparison of the intervals defining*D*_{ i }and*D*_{ k }. , initialized to*I*_{0}, is updated, ∀*j*∈ {1,...,*n*}, with either the inequality if*v*_{ j }= 1 or if*v*_{ j }= -1. - (c)
*Possible transitions from D*_{ i }*to D*_{ k }. If the calculated inequality set defines a not empty parameter space domain ⊆*PSD*_{0}then a transition towards*D*_{ k }is possible and the qualitative state*QS*(*D*_{ i }) is updated accordingly.

where the parameters are positive, and all the response functions are expressed by the Hill function commonly used in the literature.

*I*

_{0}identify the possible transitions from

*D*

_{11}towards adjacent domains in

*A*(

*D*

_{11}) = {D

_{6},

*D*

_{7},

*D*

_{12},

*D*

_{16},

*D*

_{17}}. Among the inequalities in Table 1, is the only one in agreement with

*I*

_{0}. Thus,

*QS*(

*D*

_{11}) = {

*D*

_{12}} and the only possible transition from

*D*

_{11}occurs when ∧

*I*

_{0}holds.

### Moving from a switching domain

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

#### A – Fast motion

The study of the fast dynamics is performed in *Ƶ*(*D*_{
i
}) in the scaled time *τ* = *t/q*, and aims at localizing the set of exit points in *Ƶ*(*D*_{
i
}) rather than at detailing the dynamics within it. Such points clearly identify the next adjacent domains the trajectories are moving towards from *D*_{
i
}along the *x*_{
s
}directions.

To this end, the algorithm defines, in the limit, the mapping
:
→ *Ƶ*(*D*_{
i
}), where
= *D*_{
i
}∪ *A*(*D*_{
i
}), such that the interior of *Ƶ*(*D*_{
i
}), and its boundary are the images of *D*_{
i
}, and *A*(*D*_{
i
}), respectively. More precisely, the domains *D*_{
k
}∈ *A*(*D*_{
i
}) are mapped into the faces of *Ƶ*(*D*_{
i
}) when *D*_{
k
}∈ Δ_{
s
}or into its vertices, otherwise. Let ℱ be the set of both the faces and the interior of *Ƶ*(*D*_{
i
}), and *F* its generic element.

*D*

_{7}in Fig. 2. The set ℱ has five elements, the four faces corresponding to

*D*

_{2},

*D*

_{6},

*D*

_{8},

*D*

_{12}, and the interior of

*Ƶ*(

*D*

_{7}) that correspond to

*D*

_{7}. Moreover, the vertices of

*Ƶ*(

*D*

_{7}) are the images, through the mapping Σ, of the adjacent regular domains

*D*

_{1},

*D*

_{11},

*D*

_{3}, and

*D*

_{13}(Fig. 3).

#### (a) 1. Symbolic state equations in D_{
i
}∈ Δ_{
s
}

The algorithm symbolically calculates the boundary-layer 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 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*_{2} and *F*_{12}:
and
. Finally, the exit point candidate set is updated with the points
and
(Fig. 3(a)).

#### (b) Constraints on parameters

*I*

_{0}, is calculated for each candidate exit point by requiring that each point is stable and fulfills other motion conditions on parameters. More precisely, the algorithm:

- 1.
analyzes the spectrum of the Jacobian matrix

**J**_{ F }, and imposes possible conditions on parameters to ensure stability of ; - 2.
imposes conditions

*I*_{k, l}on the sign of , ∀*l*∈ ℒ_{ F }, in a neighborhood of so that the motion of*Z*_{ l }heads towards ; - 3.
for those located on elements of ℱ, imposes the further condition,

*I*_{k,(0, 1)}, 0 < < 1 for each that does not take value 0 or 1.

Finally, is an exit point if holds.

Conditions 2. and 3. are easily checked, while the stability condition is checked by using concepts from graph theory, and the usual definition of stability based on the sign of the eigenvalues of **J**_{
F
}. It follows from *Assumption*
that the matrix **J**_{
F
}is 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
}.

*Z*

_{ l },

*l*= 2 imposes:

*I*

_{12,2}defined by (11) is compatible with (8), but the inequality

*I*

_{2,2}is not. Then, is removed from the exit point set. To be an exit point must satisfy the condition 3.:

Finally,
is an exit point if
, defined by *I*_{0} ∧ *I*_{12,2} ∧ *I*_{12,(0,1)}, holds.

As for vertices in the example, the conditions 1. and 2. are fulfilled in the point
= (0, 1) that corresponds to the vertex defined as image of *D*_{11} by the map Σ.

Finally, the only exit domains are *D*_{12}, *D*_{11}, and then *QS*(*D*_{7}) = {*D*_{12}, *D*_{11}} (Fig. 3(b)).

#### B – Slow motion

The slow motion of regular variables *x*_{
r
}along the exit points is studied in the normal time *t* in the usual frame of reference, and it is reconstructed from the reduced system (5) through the same symbolic procedure given for regular domains.

### Qualitative simulation outcomes

Through the described iterative procedure, the algorithm builds (i) a behavior tree, rooted in the initial domain, and (ii) sets of inequalities on parameter values that characterize specific state transitions. The behavior tree captures the entire spectrum of network dynamical behaviors that are invariant for ranges of parameter values defined by the calculated inequalities.

*D*

_{1}with the parameter space defined by the inequalities

*I*

_{0}(8). Through the described iterative procedure, the algorithm builds the behavior tree, BT(

*D*

_{1}), showed in Fig. 4, and calculates the inequalities on the parameters, listed in Table 2, that are associated with each path in BT(

*D*

_{1}). As in the example

*n*= 2, the trajectories described by the tree are also easily summarized in the phase plane (Fig. 5). Three reachable stable states, located in

*D*

_{11},

*D*

_{12}and

*D*

_{5}, are identified by the final leaf of each branch in BT. As

*D*

_{12}∈ Δ

_{ s }, one of them is a singular stable point whereas the others are regular stable points. These stable states are reached by different predicted qualitative behaviors, each of them occurring under specific constraints on parameters. For example, the trajectory

*QB*

_{13}starting from

*D*

_{1}, crossing

*D*

_{6}, and reaching a stable point in

*D*

_{11}is allowed when the inequalities and

*I*

_{11}, consistent with

*I*

_{0}, hold.

*n*> 2, due to graphical difficulties, each branch in the tree is given a time representation. More precisely, a specific qualitative behavior, e.g.

*QB*

_{4}, is described by the temporal evolution of each variable (Fig. 6). The behavior in Fig. 6 abstracts all numerical trajectories of the ODE model (7) obtained with different initial values in

*D*

_{1}, and different sets of parameter values that fulfill the same inequalities associated with

*QB*

_{4}(Fig. 7).

## Discussion

The algorithm is currently under implementation. As for symbolic calculus, the implementation requires to tackle complex tasks, such as: (i) update an inequality set with an another one; (ii) check the consistency of two sets of inequalities *I*_{1} and *I*_{2}; (iii) solve systems of equations; (iv) find cycles in the Jacobian matrix. As for (iii), the original equations are multilinear in *Z*_{
s
}, but due to *Assumption*
they can be straightforward solved and analyzed for stability. Also the solution of problems (i) and (ii) benefits from *Assumption*
as the inequalities are always linear. Then, thanks to the *Assumption*
, and to algorithms proposed both by the literature and common symbolic computation package, such as Mathematica [34], the tasks (i)-(iii) are simplified and feasible. As for the task (iv), it is performed by using 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*_{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 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.

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.

## Declarations

### Acknowledgements

This research is carried out within the Interdepartmental CNR-BIOINFORMATICS 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/1471-2105/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/0022-5193(73)90208-7View ArticlePubMedGoogle Scholar - Plahte E, Mestl T, Omholt SW:
**A methodological basis for description and analysis of systems with complex switch-like 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/S1097-2765(00)00070-8View 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 Single-Cell 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/ip-syb: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.00455-08PubMed CentralView ArticlePubMedGoogle Scholar - Gentile A, D'Alessandro L, Lazzari L, Martinoglio B, Bertotti A, Mira A, Lanzetti L, Comoglio PM, Medico E:
**Met-driven 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 Reverse-Engineering 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.