Abstract
Inferring dynamic biochemical networks is one of the main challenges in systems biology. Given experimental data, the objective is to identify the rules of interaction among the different entities of the network. However, the number of possible models fitting the available data is huge, and identifying a biologically relevant model is of great interest. Nested canalyzing functions, where variables in a given order dominate the function, have recently been proposed as a framework for modeling gene regulatory networks. Previously, we described this class of functions as an algebraic toric variety. In this paper, we present an algorithm that identifies all nested canalyzing models that fit the given data. We demonstrate our methods using a well-known Boolean model of the cell cycle in budding yeast.
1. Introduction
Inferring dynamic biochemical networks is one of the main challenges in systems biology. Many mathematical and statistical methods, within different frameworks, have been developed to address this problem, see [1] for a review of some of these methods. Starting from experimental data and known biological properties only, the idea is to infer a “most likely” model that could be used to generate the experimental data. Here the model could have two parts. The first one is the static network which is a directed graph showing the influence relationships among the components of the network, where an edge from node to node implies that changes in the concentration of could change the concentration of . The other part of the model is the dynamics of the network, which describes how exactly the concentration of is affected by that of . Due to the fact that biological networks are not well understood and the available data about the network is usually limited, many models end up fitting the available information, and the criteria for choosing a particular model are usually not biologically motivated but rather a consequence of the modeling framework.
A framework that has long been used for modeling gene regulatory networks is time-discrete, finite-space dynamical systems. This includes Boolean networks [2], Logical models [3], Petri nets [4], and algebraic models [5]. The latter is a straightforward generalization of Boolean networks to multistate systems. Furthermore, in [6], it was shown that logical models as well as Petri nets could be viewed and analyzed as algebraic models. The inference methods we develop here are within the algebraic models framework. To be self-contained, we briefly describe this framework and state some of the known results that we need in this paper, see [5–8] for more details. Throughout this paper, we will be talking about gene regulatory networks; however, the methods apply to biochemical networks in general.
Suppose that the gene regulatory network that we want to infer has genes and that we have a set of state transition pairs , . The input and the output are -tuples of 0 and 1 encoding the state of genes . Real-time data points are not Boolean but could be discretized (and in particular, could be made Boolean) using different methods [9]. The goal is to find a model: such that for : Notice that, since any function over a finite field is a polynomial, each is a polynomial. An algorithm that finds all models is presented in [5]. This is done by identifying, for each gene , the set of all possible functions for . This set can be represented as the coset , where is a particular such function, and is the ideal of all Boolean polynomials that vanish on the input data set, that is, . The algorithm in [5] then proceeds to find a particular model from the model space . The chosen model, which is the normal form of in the ideal , depends on the term ordering used in the Gröbner bases computation. So different ordering of the variables (genes) might lead to the selection of different models. This presents a problem as the term ordering, which is a needed for computational reasons, clearly influence the model selection process.
Several modifications have since been presented to address this problem. For example, in [10], using the Gröbner fan of the ideal , the authors developed a method that produces a probabilistic model using all possible normal forms. Other improvements on this algorithm can be found in [11, 12].
Another approach toward improving the model selection process is by restricting the model space by requiring not only that the chosen model fits the data but also satisfing some other conditions, such as its network being sparse or scale-free, the polynomials being monomials, the dynamics of the model having some desirable properties such as fixed points are the only limit cycles (i.e., starting from any initialization, the model always reaches a steady state), or that the model is robust and stable which could roughly mean that the number of attractors in the phase space is small. In a nutshell, some but not all functions in the model space are biologically relevant, and hence restricting the space to only relevant models will improve the model selection process.
By desiring a particular property, several classes of functions have been proposed as biologically relevant functions such as biologically meaningful rules [13], certain postclasses of Boolean functions have been studied in [14], and chain functions in [15], to name few. Another class of Boolean functions, which was introduced by Kauffman et al. [16], is called (nested) canalyzing functions (NCF), where an input to a single variable exclusively determines the value of the function regardless of the values of all other variables. This is a natural characterization of “canalisation” which was introduced by geneticist Waddington [17] to represent the ability of a genotype to produce the same phenotype regardless of environmental variability. Indeed, known biological functions have been shown to be canalyzing [18, 19], and Boolean nested canalyzing networks to be robust and stable [16, 19, 20].
For the purpose of restricting the model space of all Boolean polynomial models to NCFs only, we previously studied nested canalyzing functions, gave necessary and sufficient conditions on the coefficients of a boolean polynomial function to be nested canalyzing, and showed that NCFs are nothing but unate cascade functions [7]. Furthermore, in [8], the class of all nested canalyzing functions is parameterized as the rational points of an affine algebraic variety over the algebraic closure of . This variety was shown to be toric, that is, defined by a collection of binomial polynomial equations. In this paper, we present an algorithm that restricts the model space to only nested canalyzing functions by identifying all NCFs from the model space that fit the given data set.
In the next section, we briefly recall some definitions and results from [7, 8]. Our algorithm is presented in Section 3, and its implementation in Singular is discussed in Section 4. Before we conclude this paper, we demonstrate the algorithm in Section 5, where we identify all nested canalyzing models for the cell cycle in budding yeast using time course data from the Boolean model in [21].
2. Nested Canalyzing Functions: Background
We recall some of the definitions and the results from [7, 8] that we need to make this presentation self-contained. Throughout this paper, when we refer to a function on variables, we mean that depends on all variables, that is, for , there exists such that .
Definition 1. Let be a Boolean function on variables, that is, :(i)the function is a nested canalyzing function (NCF) with respect to a permutation on the variables, canalyzing input value and canalyzed output value , for , if it can be represented in the form: (ii)the function is nested canalyzing if is nested canalyzing with respect to some permutation , canalyzing input values and canalyzed output values , respectively.
Remark 2. The definition above has been generalized to multistate functions in [22], where it is also shown that the dynamics of these functions are similar to their Boolean counterparts. In [23], the authors introduce what they called kinetic models with unate structure, which are continuous models having the canalization property, and they presented an algorithm for identifying such models.
Using the polynomial form of any Boolean function, the ring of Boolean functions is isomorphic to the quotient ring , where . Indexing monomials by the subsets of corresponding to the variables appearing in the monomial, the elements of can be written as
As a vector space over , is isomorphic to via the correspondence:
The main result in [7] is the identification of the set of nested canalyzing functions in with a subset of by imposing relations on the coordinates of its elements.
Definition 3. Let be a permutation of the elements of the set . We define a new order relation on the elements of as follows: if and only if . Let be the maximum element of a nonempty subset of with respect to the order relation . For any nonempty subset of , the completion of S with respect to the permutation , denoted by , is the set .
Note that, if is the identity permutation, then the completion is := , where is the largest element of .
Theorem 4. Let , and let be a permutation of the set . The polynomial is nested canalyzing with respect to , input value and corresponding output value , for , if and only if and, for any proper subset :
Corollary 5. The set of points in corresponding to the set of all nested canalyzing functions with respect to a permutation on , denoted by , is defined by
It was shown in [8] that is an algebraic variety, and its ideal is a binomial prime ideal in the polynomial ring , where is the algebraic closure of . Namely,
Furthermore, the variety of all nested canalyzing functions is
and its ideal is
In the next section, we identify the set with the rational points in an algebraic affine variety. This will allow us to identify all nested canalyzing functions in the model space .
3. Nested Canalyzing Models
Recall that we are given the data set . The model space could be presented by the set , where and, for , In particular, is a polynomial that interpolates the data for gene and is the ideal of points of . Furthermore, the ideal is a principal ideal in the ring : Now a polynomial if and only if , for some polynomial , say . By expanding the right-hand side and collecting terms, we get that , where for , the coefficient is determined by for all and .
The proof of the following theorem follows directly from Theorem 2.4.2 in [24].
Theorem 6. Consider the ring homomorphism: given by, for , Then is the ideal of all polynomials that fit the data set . In particular, the rational points in the variety is the set of all models that fit the data set , namely .
Since the ideal of all NCFs is , the following corollary is straightforward.
Corollary 7. The ideal of all nested canalyzing functions that fit the data set is .
Remark 8. It is clear that the model space of Boolean functions is huge, since the number of monomials grows exponentially in the number of variables. For example, if a function has 5 inputs, there are different monomials in 5 variables, and hence different Boolean functions. This clearly shows that a search for NCFs inside the model space is computationally not feasible, which justifies the need for algorithms like the one above.
4. Algorithm
In this section, we present an algorithm for identifying all nested canalyzing models from the model space of a given data set.
Input
A wiring diagram, that is, a square matrix of dimension , describing the influence relationships among the genes in the network. For each variable , a table consisting of the rows , , where are the indices of the genes that affect , as specified in the wiring diagram.
Output
For each variable, the complete list of all nested canalyzing functions interpolating the given data set on the given wiring diagram. A function is in the output if it is nested canalyzing in at least one variable order. If needed, the code can easily be modified to find only nested canalyzing functions of a particular variable order.
Algorithm 9. It is a well-known fact, that a Gröbner basis for the kernel of is a basis for intersected with the ring ([24], Theorem 2.4.2) Using a similar notation as above, the algorithm is outlined as follows:use ring ;define as ideal in ;define ; define ;compute the polynomial that generates as in (12);for each variable do (1)compute as in (11);(2)let ; its coefficients are the same as above;(3)compute a Gröbner basis for the ideal generated by the coefficients of using any elimination order to eliminate all from ;(4)concatenate generators of and ;(5)compute the primary decomposition of to obtain necessary and sufficient conditions on the coefficients of all NCFs fitting the data set ;End. An implementation of this algorithm is available as a Singular library [25, 26].
5. Application: Inferring the Cell Cycle Network in Budding Yeast
The interactions among proteins constitute complex molecular networks that regulate cell behavior such as the decision to undergo cell division. We use data points generated by a previously published model of the cell-cycle regulatory network of the budding yeast Saccharomyces cerevisiae to demonstrate the algorithm described in this paper [21]. The model is a Boolean model, where nodes represent proteins and edges describe protein interactions that are either activating or inhibiting. Proteins are updated according to a threshold rule, that is, a protein is activated, respectively deactivated, if the weighted sum of the activating input proteins is greater than, respectively less than, the weighted sum of the inhibiting input proteins. The model contains the key regulators of the cell cycle process and the known interactions among these regulators. This model captures the known features of the global dynamics of the cell cycle, it is robust and stable, and the trajectory of the known cell-cycle sequence is a stable and attracting trajectory as it has 1764 states out of the total number of 2048 states. The remaining states are distributed into 6 small trajectories.
In this section, we use the time course corresponding to the biological cell-cycle sequence, see Table 1, to infer nested canalyzing models of the cell cycle. That is, assuming the same wiring diagram as the threshold model in [21], we use our algorithm to identify, for each gene in the network, all nested canalyzing functions that fit the cell cycle sequence.
We start by describing the model. There are twelve nodes that represent eleven proteins and a start signal. The proteins are members of the following three classes: cyclins (Cln1,2, Cln3, Clb1,2, and Clb5,6), inhibitors, degraders, and competitors of cyclin complexes (Sic1, Cdh1, Cdc20, and Cdc14), and transcription factors (SBF, MBF, Swi5, and Mcm1/SFF). This simplified network (Figure 1) is almost identical to the network in [21], where the only difference is that we do not force self-degradation, as it was added to some nodes in the network because they did not have inhibitors, but without biological justification [21]. Furthermore, we do not impose activation or inhibition in the network. As we do not use threshold functions but more general boolean functions, a variable can both increase and decrease the concentration of another substrate, depending on the concentrations of other proteins.
Li et al. [21] use their model to generate a time course of temporal evolution of the cell-cycle network, shown in Table 1. This time course is in agreement with the behavior of the cell-division process cycling through the four distinct phases , (synthesis), , and (mitosis).
We used this time course along with the network as the only input to the algorithm to obtain all nested canalyzing functions that interpolate the time course. In the fifth column of Table 2, we list the number of NCFs for each protein. By requiring the Boolean function to be nested canalyzing, we have significantly reduced the number of possible functions for each protein as it is evident when comparing the numbers in the third and fifth columns of Table 2. However, even after this reduction, there are possible nested canalyzing models that fit the time course in Table 1. To reduce this number more, one needs to use additional time courses or request that the models incorporate additional biological information about yeast cell cycle.
5.1. Dynamics
To analyze the dynamics of the resulting nested canalyzing models, we randomly sampled 2000 models and analyzed them. The average number of basins of attraction (components) per network is 3.09, and the average size of the component containing the given trajectory is 1889. In only 6 models, the trajectory in Table 1 is not in the largest component; however, the average size of the component containing the trajectory is 833.5.
These results clearly show that nested canalyzing models for the cell cycle network are in agreement with the original threshold model of Li et al., and since such models are known to be robust and stable, any of these models could be used as a model for the cell cycle in budding yeast. Furthermore, especially when there is no evidence for choosing a particular type of functions, a nested canalyzing function has an advantage over other possible choices.
5.2. Comparison with Random Networks
To understand the effect of the network itself on its dynamics, we sampled 2000 models on the same network, where the local function of each gene in each one of these models is chosen randomly from all possible functions in the model space. We found that the given cell cycle trajectory has oftentimes much smaller basin of attraction, and hence random functions on the cell cycle network could not in general produce the desired dynamics. A comparison of the statistics from the sampled networks is shown in Figures 2 and 3.
6. Conclusion
In this paper, we have presented an algorithm for identifying all Boolean nested canalyzing models that fit a given time course or other input-output data sets. Our algorithm uses methods from computational algebra to present the model space as an algebraic variety. The intersection of this variety with the variety of all NCFs, which was parameterized in [8], gives us the set of all NCFs that fit the data. We demonstrate our algorithm by finding all nested canalyzing models of the cell cycle network form Li et al. [21]. We then showed that the dynamics of almost any of these models is strikingly similar to that of the original threshold model. Unless the chosen model is required to meet other conditions, and in that case the model space will be reduced further, any one of the models that our algorithm found is an acceptable model of the cell cycle process in the budding yeast.
One limitation of the current algorithm, which we left for future work, is that it does not distinguish between activation and inhibition in the network as we do not have a systematic method of knowing when a given variable in a (nested canalyzing) polynomial is an activator or inhibitor.
As our algorithm relies heavily on different Gröbner-based computations, the current implementation in Singular allows a given gene to have at most 5 regulators. This is due to the fact that the number of monomials then is 32 which is already a burden especially when the primary decomposition of an ideal is what we are after. We are working on a better implementation so that we can infer larger and denser networks.
Acknowledgments
F. Hinkelmann was supported by a grant from the U.S. Army Research Office. The authors thank the referees for their comments and suggestions.