Modeling of gene regulatory networks : A review

Gene regulatory networks play an important role the molecular mechanism underlying biological processes. Modeling of these networks is an important challenge to be addressed in the post genomic era. Several methods have been proposed for estimating gene networks from gene expression data. Computational methods for development of network models and analysis of their functionality have proved to be valuable tools in bioinformatics applications. In this paper we tried to review the different methods for reconstructing gene regulatory networks.


INTRODUCTION
A gene regulatory network or genetic regulatory network (GRN) is a collection of DNA segments in a cell which interact with each other indirectly (through their RNA and protein expression products) and with other substances in the cell, thereby governing the rates at which genes in the network are transcribed into mRNA.GRNs provide a systematic understanding of molecular mechanisms underlying biological processes [1][2][3][4][5][6][7].The groups of genes, regulatory proteins and their interactions are often referred to as regulatory networks, whereas the complete set of metabolites and the enzyme-driven reactions constitute the metabolic networks.The nodes of this network are genes and the edges between nodes represent gene interactions through which the products of one gene affect those of another.These interactions can be inductive (the arrowheads), with an increase in the expression of one leading to an increase in the other, or inhibitory (the filled circles), with an increase in one leading to a decrease in the other.A series of edges indicates a chain of such dependences, with cycles corresponding to feedback loops.
Gene regulatory networks play a vital role in organism development by controlling gene expression.Understanding the structure and behavior of gene regulatory network is a fundamental problem in biology.With the availability of gene expression data and complete genome sequences, several novel experimental and computational approaches have recently been developed which helps to comprehensively characterize these regulatory networks by enabling the identification of their genomic or regulatory state components.Accurate prediction of the behavior of regulatory networks will also accelerate biotechnological projects and such predictions are quicker and cheaper than lab experiments.
Creating accurate dynamic models of GRNs is gaining importance in biomedical research and development.Gene expression microarrays monitor the transcription activities of thousands of genes simultaneously, which provides great opportunities to explore large scale regulatory networks.Constructing a GRN from expression data, a process which is called reverse-engineering, is not a computationally simple problem because an enormous amount of time is needed even when a trivial approach is applied.Various computational models developed for regulatory network analysis can be roughly divided into four classes (Figure 1).The first class 1) logical models, describes regulatory networks qualitatively.They allow users to obtain a basic understanding of the different functionalities of a given network under different conditions.Their qualitative nature makes them flexible and easy to fit to biological phenomena, although they can only answer qualitative questions.To understand and manipulate behaviors that depend on finer timing and exact molecular concentrations, a second class of models was developed 2) continuous models.For example, to simulate the effects of dietary restriction on yeast cells under different nutrient concentrations, users must resort to the finer resolution of continuous models.A third class of models was introduced following the observation that the functionality of regulatory networks is often affected by noise.As the majority of these models account for interactions between individual molecules, they are referred to 3) single molecule level models.The fourth class includes 4) hybrid models combining different techniques like neural networks and fuzzy rules.
A complete gene regulatory network model incorporates experimental knowledge about the components and their interactions as well as the initial state of these components, and leads to the known final state or dynamical behavior of the network.Validated models then are able to investigate cases that cannot be explored experimentally, for example changes in the initial state, in the components or in the interactions, and they can lead to predictions and insights into the functioning of the system robust is the system under extreme conditions.In this article we review the various modeling techniques for reconstructing gene regulatory network.

MODELLING TECHNIQUES
Figure 1 illustrates various Gene Regulatory Network construction models that are discussed in following sections.

Logical Models
The most basic and simplest modeling methodology is discrete and logic-based, and was introduced by Kauffman and Thomas [8,9].The reconstruction of the regulatory network that controls the development of sea urchin embryos is a seminal example of the profound insights that qualitative examination of regulatory network models can provide.This work demonstrates how maternal cues initiate the activity of the regulatory network and how this network orchestrates the developmental process.Logical models represent the local state of each entity in the system (for example, genes, proteins and small molecules) at any time as a discrete level, and the temporal development of the system is often assumed to occur in synchronous, discrete time steps.Entity levels are updated at each time step according to regulation functions.Discrete modeling allows researchers to rely on purely qualitative knowledge.Such models can be analyzed using a broad range of well established mathematical and statistical methods.

Boolean Network
Boolean networks are a dynamic model of synchronous interactions between nodes in a network.They are the simplest network models that exhibit some of the biological and systemic properties of real gene networks [10,11].Because of the simplicity they are relatively easier to interpret biologically.
A Boolean network is a directed graph G(X, E), where the nodes, x i ∈ X, are Boolean variables.To each node, where the arguments are all and only the parent nodes of x i in G. Together, at any given time, the states (values) of all nodes represent the state of the network, given by the vector   . For gene networks the node variables correspond to levels of gene expression, discretized to either up or down [12][13][14].The Boolean functions at the nodes model the aggregated regulation effect of all their parent nodes.The states of all nodes are updated at the same time (i.e., synchronously) according to their respective Boolean functions: All states' transitions together correspond to a state transition of the network from S(t) to the new network state, S(t + 1).A sample network is shown in Figure 2.
LIMITATION: These models are ultimately limited by their definition: they are Boolean and synchronous.In reality, of course, the levels of gene expression do not have only two states but can assume virtually continuous values.Thus discretization of the original data becomes a critical step in the inference, and often reducing the values to two states may not suffice.In addition, the updates of the network states in this model are synchronous, whereas biological networks are typically asynchronous.Finally, despite their simplicity, only small nets can be reverse engineered with the current state-of-the-art algorithms.

Probabilistic Boolean Network
Often, due to insufficient experimental evidence or in-Figure 2.An example Boolean network and three possible ways to represent it.The one on the left is a gene network modeled as a Boolean network, in the middle is a wiring diagram obviating the transitions between network states, and on the right is a truth table of all possible state transitions.complete understanding of a system, several candidate regulatory functions may be possible for an entity.This raises the need to express uncertainty in the regulatory logic.Shmulevich et al., [15,16] addressed this idea by modifying the Boolean network model such that an entity can have several regulation functions, each of which is given a probability based on its compatibility with prior data.At each time step, every entity is subjected to a regulation function that is randomly selected according to the defined probabilities.Hence the model is stochastic and an initial global state can lead to many trajectories of different probabilities.The new model, the probabilistic Boolean network (PBN), generates a sequence of global states that constitutes a Markov chain.For example, a PBN was used to model a 15 gene sub network that was inferred from human glioma expression data [15,16].This analysis demonstrates that the stationary distributions of entities may indicate possible regulatory relationships among them: entities that have the same states in a significant proportion of the global states are likely to be related.As the number of global states in the gene sub network was prohibitively large, one study estimated the stationary distribution by sampling the global states.
LIMITATION: Even though it is stochastic the state space is discrete.

Bayesian Network
The basic of Bayesian Network is Bayes' Theorem.It can be described as follows.Let X be a data sample whose class label is unknown.Let H be a hypothesis that X belongs to class C. For classification problems, determine P(H/X): the probability that the hypothesis holds given the observed data sample X.It is called posteriori probability.
P(H): prior probability of hypothesis H (i.e., the initial probability before we observe any data, reflects the background knowledge).
P(X): probability that sample data is observed.P(X|H): probability of observing the sample X, given that the hypothesis holds.
Given training data X, posteriori probability of a hypothesis H, P(H|X) follows the Bayes theorem:

P X H P H P H X P X 
A simple Bayesian Classifier will work as follows: Let D be a training set of tuples and their associated class labels.As usual, each tuple is represented by an n-dimensional attribute vector, , depicting n measurements made on the tuple from n attributes, respectively, 1 2 , , , n A A A  .Suppose that there are m classes, .
Given a tuple, X, the classifier will predict that X belongs to the class having the highest posterior probability, conditioned on X.That is, the naïve Bayesian classifier predicts that tuple X belongs to the class C i if and only if Thus we maximize P(C i |X).The class C i for which Bayesian classifiers assume that the effect of an attribute value on a given class is independent of the values of the other attributes.This assumption is called class conditional independence.It is made to simplify the computations involved and, in this sense, is considered "naïve".Bayesian belief networks are graphical models, which unlike naïve Bayesian classifiers allow the representation of dependencies among subsets of attributes.
Bayesian networks are a class of graphical probabilistic models.Formally a Bayesian network [17,18] is a joint probability distribution over a set of random variables.They combine two very well developed mathematical areas: probability and graph theory.A Bayesian network consists of an annotated directed acyclic graph G(X, E), where the nodes, x i  X, are random variables representing genes' expressions and the edges indicate the dependencies between the nodes.The random variables are drawn from conditional probability distributions is the set of parents for each node.A Bayesian network implicitly encodes the Markov Assumption that given its parents, each variable is independent of its non-descendants.With this assumption each Bayesian network uniquely specifies a decomposition of the joint distribution over all variables down to the conditional distributions of the nodes: , , A belief network is defined by two components, a directed acyclic graph and a set of conditional probability tables [19].Each node in the directed acyclic graph represents a random variable.The variables may be discrete or continuous-valued.They may correspond to actual attributes given in the data or to "hidden variables" believed to form a relationship.If an arc is drawn from a node Y to a node Z, then Y is a parent or immediate predecessor of Z and Z is a descendant of Y.Each variable is conditionally independent of its non descendants in the graph, given its parents.
For example, let us consider the five variables in Fig- ure 3. Without using any independence assumptions, the joint probability distribution can be written as: , In contrast, using the independence assumptions implied by the network in Figure 3, the same distribution can be expressed as:

, , , P A B C D E P E P A P B A E P D A P C B 
. If the variables are all binary in this network, the former form requires 31 parameters, while the latter only needs 10 parameters.More generally, if G is defined over N binary variables and their maximal number of parents is bound by M, then instead of using 2N − 1 independent parameters to represent the full joint probability distribution, a Bayesian network model can represent the same joint distribution with at most 2MN parameters.
A node within the network can be selected as an "output" node, representing a class label attribute.There may be more than one output node.Various algorithms for learning can be applied to the network.Rather than returning a single class label, the classification process can return a probability distribution that gives the probability of each class.A major advantage of Bayesian network models is the ability to learn them from observed data.Bayesian networks can capture linear, non-linear, combinatorial, stochastic and other types of relationships among variables.They are suitable for modeling gene networks because of their ability to represent stochastic events, to describe locally interacting processes, to handle noisy or missing biological data in a principled statistical way and to possibly make causal inferences from the derived models [20,21].Hence, Bayesian networks, including their variants Dynamic Bayesian networks, Gaussian networks, Module networks, mixture Bayesian networks and state-space models (SSMs), etc., have become widely used tools for regulatory-network modeling.
LIMITATION: Although effective in dealing with noise, incompleteness and stochastic aspects of gene regulation, they fail to consider temporal dynamic aspects that are an important part of regulatory networks modeling.Dynamic Bayesian networks (DBN) evolved feedback loops to effectively deal with the temporal aspects of regulatory networks but their benefits are hindered by the high computational cost required for learning the conditional dependencies in the cases where large numbers of genes are involved.

Continuous Models
Biological experiments usually produce real, rather than discrete valued, measurements.Examples include reaction rates, cell mass [22][23][24][25], cell cycle length and gene expression intensities.Logical models require discretization of the real valued data, which reduces the accuracy of the data.Continuous models, using real valued parameters over a continuous timescale, allow a straightforward comparison of the global state and experimental data and can theoretically be more accurate.In practice, however, quantitative measurements are almost always partial (that is, they cover only a fraction of the system's entities).Therefore, some of the parameters of continuous models are usually based on estimations or inference.

Linear Model
The defining property of linear models is that each regulator contributes to the input of the regulation function independently of the other regulators, in an additive manner [10].In other words, the change in the level of each entity depends on a weighted linear sum of the levels of its regulators.This assumption allows a high level of abstraction and efficient inference of network structure and regulation functions.
A biological system can be considered to be a state machine, where the change in internal state of the system depends on the current internal state plus any external inputs.The mRNA levels form an important part of the internal state of a cell (ideally, we also want to measure protein levels, metabolites, etc.).As a first approximation, we fit the expression data with a purely linear model, where the change in expression level of each mRNA species is derived as a weighted sum of the expression levels of all other genes.Of course, a linear model can never be much more than a caricature of the real system, but perhaps we can still draw some interesting conclusions from it.
The basic linear model is of the form , where X i (t + Δt) is the expression level of gene i at time t + Δt, and W ij indicates how much the level of gene j inu-ences gene i.For each gene, we will also add an extra term indicating the influence of kainate, and a constant bias term to model the activation level of the gene in the absence of any other regulatory inputs.The differences in gene regulation due to tissue type will be modeled by a difference in bias.The final formula becomes: where kainate(t) is the kainate level at time t, K i is the influence of kainate on gene i, C i is a constant bias factor for each gene, and T i indicates the difference in bias between tissue types (T i = 0 when simulating spinal cord, so the total bias for spinal cord is C i , for hippocampus C i + T i ).LIMITATION: Linear additive regulation models revealed certain linear relations in regulatory systems but failed to capture nonlinear dynamics aspects of genes regulation.When higher sensitivity to detail is desired, more complex models are preferable.

Differential Equation Based Model
Differential equation models encode a gene network as a system of differential equations.Difference and differential equations allow more detailed descriptions of network dynamics, by explicitly modelling the concentration changes of molecules over time [26,27].
The basic difference equation model is of the form

  g t t g t w g t w g t t g t t g t w g t w g t t
where g i (t + Δt) is the expression level of gene i at time t + Δt, and w ij the weight indicating how much the level of gene i is influenced by gene j   , 1, , i j n   .Note that this model assumes a linear logic control model-the expression levels of genes at a time t + Δt, depends linearly on the expression levels of all genes at a time t.For each gene, one can add extra terms indicating the influence of additional substances.Differential equation models are similar to difference equation models, but follow concentration changes continuously, modelling the time difference between two time steps in infinitely small time increases, i.e.Δt is approaching 0.
Difference and differential models depend on numerical parameters, which are often difficult to measure experimentally.An important question for these models is stability-does the behaviour of the system depend on the exact values of these parameters and initial substance concentrations, or is it similar for different variations.It seems unlikely that an unstable system represents a biologically realistic model, while on the other hand, if the system is stable, the exact values of some parameters may not be essential.
The rate of change in concentration of a particular transcript is given by an influence function of other RNA concentrations.The non-linear differential equations describe the mutual activating and repressing influences of genes in a GRN at a high-level of abstraction.In particular, it is assumed that the rate of gene expression depends exclusively on the concentration of gene products arising from the nodes (genes) of the GRN.This means that the influence of other molecules (e.g., transcription factors) and cellular processes (translation) is not taken into account directly.Even with these limitations, dynamic GRN models of this kind can be useful in deciphering basic aspects of gene-regulatory interactions.
One major advantage of all three methods described below lies in their simple homogeneous structures, as this allows the settings of parameter discovering software to be easily customized for these structures.The three methods describe dynamic GRN models by means of a system (or set) of ordinary differential equations.For a GRN comprising N genes, N differential equations are used to describe the dynamics of N gene product concentrations, X i with 1, , i N   .In all three methods, the expression rate dX i /dt of a gene product concentration may depend on the expression level of one or more gene products of the genes X j , with .Thus, the gene product concentration X i may be governed by a self-regulatory mechanism (when i = j), or it may be regulated by products of other genes in the GRN.The three modeling methods differ in the way they represent and calculate expression rates.
1, , j  N 2.2.2.1.The Artificial Neural Network (ANN) Method Vohradsky [28] introduced ANNs as a modeling method capable of describing the dynamic behavior of GRNs.The way this method represents and calculates expression rates depends on the weighted sum of multiple regulatory inputs.This additive input processing is capable of representing logical disjunctions.The expression rate is restricted to a certain interval where a sigmoidal transformation maps the regulatory input to the expression interval.ANNs provide an additional external input which has an influence on this transformation in that it can regulate the sensitivity to the summed regulatory input.Finally, the ANN method defines the degradation of a gene product on the basis of standard mass-action kinetics.Formally, the ANN method is defined as: The parameters of the ANN method have the following biological interpretations: N: Number of genes in the GRN to be modeled.The genes of the GRN are indexed by i and j, where , 1, , i j N   .ϑ i : Influence of external input on gene i, which modulates the gene's sensitivity of response to activating or repressing influences.
f: Represents a non-linear sigmoid transfer function modifying the influence of gene expression products Xj and external input ϑ i to keep the activation from growing without bounds.
k i : Degradation of the i-th gene expression product.
The mathematical properties of the ANN method have been well studied because it is a special case of a recurrent neural network.In particular, the symmetry of the matrix of connection weights w ij influences whether the network dynamics are oscillatory or whether they converge on a steady (or even chaotic) state.High positive or negative values of the external input, ϑ i , reduce the effect of the connection weights.This is explored in Case D where ϑ i has been interpreted as a delay to the reaction kinetics of the transcriptional machinery.

The S-System (SS) Method
Savageau [29] proposed the synergistic system or S-system (SS) as a method to model molecular networks.When modeling GRNs with the SS method, the expression rates are described by the difference of two products of power-law functions, where the first represents the activation term and the second degradation term of a gene product X i .This multiplicative input processing can be used to define logical conjunctions for both the regulation of gene expression processes and for the regulation of degradation processes.The SS method has no restrictions in the gene expression rates and thus does not implicitly describe saturation.Formally, the SS method is defined as:

R
The parameters of the SS method have the following biological interpretations: N: Number of genes in the GRN to be modeled.The genes of the GRN are indexed by i and j, where ., 1, , i j N   α i : Rate constant of activation term; in SS GRN models, all activation (up-regulation) processes of a gene i are aggregated into a single activation term.
β i : Rate constant of degradation term; in SS GRN models, all degradation processes of a gene i are aggregated into a single degradation term.
g ij ,h ij : Exponential parameters called kinetic order.These parameters describe the interactive influences of gene j on gene i.Positive values of g ij indicate an activating influence on the expression of gene i, whereas inhibiting influences are represented by negative values.
Similarly, positive values of h ij indicate increasing degradation of the gene product X i , whereas decreasing degradation is represented by negative values.The parameters used in SS models have a clear physical meaning and can be measured experimentally, yet they describe phenomenological influences, as opposed to stoichiometric rate constants in general mass action (GMA) systems.The SS method generalizes mass-action kinetics by aggregating all individual processes into a single activition and a single degradation term (per gene).In contrast, the GMA system defines all individual processes k with 1, , k  with the sum of power-law functions according to: The parameters of the GMA system have the following biological interpretations: α i : Rate constant of activation process k. β ik : Rate constant of degradation process k. g ijk : Exponential parameter called kinetic order describing the interactive influence of X j on gene i of process k.
h ijk : Exponential parameter called kinetic order describing the interactive influence of X j on gene i of process k.

The General Rate Law of Transcription (GRLOT)
Method The GRLOT method has been used to generate benchmark time-series data sets to facilitate the evaluation of different reverse-engineering approaches.GRLOT models multiply individual regulatory inputs.Activation and inhibition are represented by different functional expressions that are similar to Hill kinetics, which allow the inclusion of cooperative binding events.Identical to the ANN, the degradation of gene products is defined via mass-action kinetics.Formally, the GRLOT method is defined as: The parameters of the GRLOT method have the following biological interpretations: Failed to capture nonlinear dynamics aspects of genes regulation.
Not sufficient if higher sensitivity to detail is desired.

Differential Equation Based Model
Simple homogeneous structures: this allows the settings of parameter discovering software to be easily customized for these structures.
Involve a large number of parameters-O(d 2 ) parameters where d is the number of genes modeled.

Single Molecule Level Model
The most detailed, can capture stochasticity.computationally expensive

Hybrid Model
In the real world systems both continuous aspects and discrete aspects are present.Hybrid models helps in modeling both together.
Computationally expensive number of activators A can be related to the total number of genes by I + A ≤ N. K ij : Concentration at which the influence of inhibitor j is half of its saturation value.K ak : Concentration at which the influence of activator k is half of its saturation value.
n j , n k : Regulate the sigmoidicity of the interaction behavior in the same way as Hill coefficients in enzyme kinetics.
k i : Degradation of the i-th gene expression product.LIMITATIONS: Unless they are restricted to simple function forms, differential equation models involve a large number of parameters-O(d 2 ) parameters where d is the number of genes modeled.Moreover, differential equation models require time-series data to learn the parameters

Single Molecule Level Model
Every biological network is composed of stochastic components, and therefore it may manifest different behaviours, even starting from the same initial conditions [30,31].When the number of involved molecules of each species is large, the law of mass action can be used to accurately calculate the change in concentrations, and little or no stochastic effect is observable.However, when the number of molecules is small, significant stochastic effects may be seen.This is particularly true for regulatory networks, in which the number of regulatory molecules is often low [32][33][34][35].Recently, single cell experimental assays demonstrated the stochastic behaviour of the processes of transcription and translation [36].

Hybrid Model
In the real world systems both continuous aspects and discrete aspects are present.In general, concentrations are expressed as continuous values, whereas the binding of a transcription factor to DNA is expressed as a discrete event (bound or unbound).However, the boundaries between the discrete and continuous aspects depend on the level of detail that our model is designed for.For instance, on single cell level the concentrations may have to be expressed by molecule counts and become discrete, whereas if we use thermodynamic equilibrium to model the protein-DNA binding, the variable describing the binding state becomes continuous.Hybrid models have been developed in an attempt to describe both, discrete and continuous aspects in one model.
An example of a hybrid model [37,38] is a multi-layer evolutionary trained neuro-fuzzy recurrent network (ENFRN) applied to the problem of GRN reconstruction, which addresses the major drawbacks of currently existing computational methods.This choice was driven by the benefits, in terms of computational power, that neural network based methods provide.The self-organized nature of ENFRN algorithm is able to produce an adaptive number of temporal fuzzy rules that describe the relationships between the input (regulating) genes and the output (regulated) gene.Related to that, another advantage of this approach is that it overcomes the need of prior data discretization, a characteristic of many computational methods which often leads to information loss.The dynamic mapping capabilities emerging from the recurrent structure of ENFRN and the incorporation of fuzzy logic drive the construction of easily interpretable fuzzy rules of the form: "IF gene x is highly expressed at time t THEN its dependent/target gene y will be lowly expressed at time t + 1".The evolutionary training, based on the PSO framework, tries to avoid the drawbacks of classical neural networks training algorithms [39].Additionally, we are approaching the under-determinism problem by selecting the most suitable set of regulatory genes via a time-effective procedure embedded in the construction phase of ENFRN.Also, besides determining the regulatory relations among genes, this method can determine the type of the regulation (activation or repression) and at the same time assign a score, which might be used as a measure of confidence in the retrieved regulation.
Comparison of different models discussed in this paper is given in Table 1.

CONCLUSION
In this paper we have reviewed the different modeling methods for reconstructing gene networks from gene expression data.All methods mentioned above are for reverse engineering of GRNs from gene expression data.The Boolean network models have the limitation of discrete apace and in reality, of course, the levels of gene expression do not have only two states but can assume virtually continuous values.The probabilistic methods have the flexibility of assuming different probability of expression for gene at a particular point of time and are closely related to real time situations.Also we discussed continuous models like linear and differential models using non-discrete values.Single molecule based models consider stochastic behavior of biological network and hybrid models combines different concepts for GRN reconstruction.
v i : Maximal expression rate of gene i.w ij : The connection weight or strength of control of gene j on gene i.Positive values of w ij indicate activating influences while negative values define repressing influences.

Table 1 .
Maximal expression rate of gene i.I Advantages and disadvantages of the different algorithms for gene network construction.
j : Inhibitor (repressor) j.A k : Activator k; the number of inhibitors I, and the OPEN ACCESS