Inferring Biologically Relevant Models: Nested Canalyzing Functions

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 (𝐬_𝑗,𝐭_𝑗), 𝑗=1,…,𝑟. 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:𝑓=𝑓₁,𝑓₂,…,𝑓_𝑛∶𝔽^𝑛2⟶𝔽^𝑛2,(1) such that for 𝑗=1,…,𝑟:𝑓𝐬_𝑗=𝑓₁𝐬_𝑗,…,𝑓_𝑛𝐬_𝑗=𝐭_𝑗.(2) 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 𝑖=1,…,𝑛, there exists (𝑎₁,…,𝑎_𝑛)∈𝔽^𝑛2 such that ℎ(𝑎₁,…,𝑎_𝑖−1,𝑎_𝑖,𝑎_𝑖+1,…,𝑎_𝑛)≠ℎ(𝑎₁,…,𝑎_𝑖−1,1+𝑎_𝑖,𝑎_𝑖+1,…,𝑎_𝑛).

Definition 1. Let ℎ be a Boolean function on 𝑛 variables, that is, ℎ∶𝔽^𝑛2→𝔽₂:(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 𝑖=1,…,𝑛, if it can be represented in the form: ℎ𝑥₁,𝑥₂,…,𝑥_𝑛=⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩𝑏¹𝑏₂if𝑥_𝜎(1)=𝑎₁,𝑏₃if𝑥_𝜎(1)≠𝑎₁,𝑥_𝜎(2)=𝑎₂,if𝑥𝑥_𝜎(1)≠𝑎_1𝜎(3),𝑥_𝜎(2)≠𝑎₂,𝑏⋮⋮_𝑛=𝑎₃,i𝑥f𝑥_{𝜎(1)𝜎(𝑛)}≠𝑎₁,…,𝑥_{𝜎(𝑛−1)}≠𝑎_𝑛−1,=𝑎_𝑛,𝑏_𝑛if𝑥_𝜎(1)≠𝑎₁,…,𝑥_𝜎(𝑛)≠𝑎_𝑛;(3)(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 𝐽=⟨𝑥^2𝑖−𝑥_𝑖∶1≤𝑖≤𝑛⟩. Indexing monomials by the subsets of [𝑛]∶={1,…,𝑛} corresponding to the variables appearing in the monomial, the elements of 𝑅 can be written as𝑅=_{𝑆⊆[𝑛]}𝑐_𝑆_𝑖∈𝑆𝑥_𝑖∶𝑐_𝑆∈𝔽₂.(4)

As a vector space over 𝔽₂, 𝑅 is isomorphic to 𝔽^2𝑛₂ via the correspondence:𝑅∋_{𝑆⊆[𝑛]}𝑐_𝑆_𝑖∈𝑆𝑥_𝑖⟷𝑐_∅,…,𝑐_[𝑛]∈𝔽^2𝑛₂.(5)

The main result in [7] is the identification of the set of nested canalyzing functions in 𝑅 with a subset 𝑉^ncf of 𝔽^2𝑛₂ 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 [𝑟^𝜎𝑆]={𝜎(1),𝜎(2),…,𝜎(𝑟^𝜎𝑆)}.
Note that, if 𝜎 is the identity permutation, then the completion is [𝑟_𝑆]:= {1,2,…,𝑟_𝑆}, 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 𝑖=1,…,𝑛, if and only if 𝑐_[𝑛]=1 and, for any proper subset 𝑆⊆[𝑛]: 𝑐_𝑆=𝑐_[𝑟^𝜎𝑆_{]𝑟𝜎(}_{𝑖)∈𝜎𝑆⧵𝑆}𝑐_{[𝑛]⧵{𝜎(𝑖)}}.(6)

Corollary 5. The set of points in 𝔽^2𝑛₂ corresponding to the set of all nested canalyzing functions with respect to a permutation 𝜎 on [𝑛], denoted by 𝑉^ncf_𝜎, is defined by 𝑉^ncf_𝜎=⎧⎪⎨⎪⎩𝑐_∅,…,𝑐_[𝑛]∈𝔽^2𝑛₂∶𝑐_[𝑛]=1,𝑐_𝑆=𝑐_[𝑟^𝜎𝑆_]×_{𝑟𝜎(𝑖)∈𝜎𝑆⧵𝑆}𝑐_{[𝑛]⧵{𝜎(𝑖)}},for𝑆⊆[𝑛]⎫⎪⎬⎪⎭.(7)

It was shown in [8] that 𝑉^ncf_𝜎 is an algebraic variety, and its ideal 𝕀(𝑉^ncf_𝜎) is a binomial prime ideal in the polynomial ring 𝔽₂[{𝑐_𝑆∶𝑆⊆[𝑛]}], where 𝔽₂ is the algebraic closure of 𝔽₂. Namely,𝐼_𝜎==𝕀𝑐𝑉^ncf_𝜎[𝑛]−1,𝑐_𝑆−𝑐_[𝑟^𝜎𝑆_{]𝑟𝜎(}_{𝑖)∈𝜎𝑆⧵𝑆}𝑐_{[𝑛]⧵{𝜎(𝑖)}}∶𝑆⊂[𝑛].(8)

Furthermore, the variety of all nested canalyzing functions is𝑉^ncf=_𝜎𝑉^ncf_𝜎,(9)

and its ideal is𝕀𝑉^ncf=_𝜎𝐼_𝜎.(10)

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 𝐷={(𝐬₁,𝐭₁),…,(𝐬_𝑟,𝐭_𝑟)}⊂𝔽^𝑛2×𝔽^𝑛2. The model space could be presented by the set 𝑓+𝐼, where 𝑓=(𝑓₁,…,𝑓_𝑛) and, for 𝑖=1,…,𝑛,𝑓_𝑖𝑥₁,…,𝑥_𝑛=^𝑟_𝑗=1𝑡_{𝑛𝑗,𝑖}_𝑒=11−𝑥_𝑒−𝑠_𝑗,𝑒.(11) 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 𝑅/𝐽:𝐼==𝕀𝐬₁,…,𝐬_𝑟_𝑟_𝑗=1=𝕀𝐬_𝑗_𝑟_𝑗=1𝑥₁=_𝑟−𝑠_𝑗,1,…,𝑥_𝑛−𝑠_𝑗,𝑛_𝑗=11−^𝑛_𝑒=11−𝑥_𝑒=−𝑠_𝑗,𝑒_𝑟_𝑗=11−^𝑛_e=11−𝑥_𝑒−𝑠_𝑗,𝑒.(12) Now a polynomial ℎ∈𝑓_𝑖+𝐼 if and only if ℎ=𝑓_𝑖+𝑔(𝑥₁,…,𝑥_𝑛)∏^𝑟𝑗=1(1−∏^𝑛𝑒=1(1−(𝑥_𝑒−𝐬_𝑗,𝑒))), 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 𝑗=1,…,𝑟.

The proof of the following theorem follows directly from Theorem  2.4.2 in [24].

Theorem 6. Consider the ring homomorphism: Φ∶𝔽₂𝑐_𝑆∶𝑆⊆[𝑛]⟶𝔽₂𝑏_𝐻∶𝐻⊆[𝑛](13) given by, for 𝑆⊆[𝑛], 𝑐_𝑆⟼𝑊_𝑆𝑏_𝐻,𝐬_𝑗,𝐭_𝑗.(14) Then ker(Φ) is the ideal of all polynomials that fit the data set 𝐷. In particular, the rational points in the variety 𝕍(ker(Φ)) is the set of all models that fit the data set 𝐷, namely 𝑓+𝐼.

Since the ideal of all NCFs is 𝕀(𝑉^ncf), the following corollary is straightforward.

Corollary 7. The ideal of all nested canalyzing functions that fit the data set 𝐷 is 𝕀(𝑉^ncf)+ker(Φ).

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 2⁵=32 different monomials in 5 variables, and hence 2³²=4,294,967,296 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 (s_𝑗,𝑖1,…,s_{𝑗,𝑖𝑠},t_𝑗,𝑖), 𝑗={1,…,𝑟}, 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 𝕀(𝑉^ncf) 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 𝕀(𝑉^ncf);(5)compute the primary decomposition of 𝐺+𝕀(𝑉^ncf) 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.

Table 1 

The temporal evolution of the Boolean cell-cycle model in [21]; corresponding to the biological 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.

Figure 1 

The simplified cell cycle network in budding yeast, which is based on the model in [21].

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 330,559,488 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.

Table 2 

For each protein 𝑖, we list the number of inputs, the number of possible Boolean functions (the cardinality of 𝑓_𝑖+𝐼), the number of nested canalyzing functions with the given number of inputs, and finally the number of nested canalyzing functions in the model space 𝑓_𝑖+𝐼.

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.

Figure 2 

Nested canalyzing functions with the wiring diagram in Figure 1 interpolating the time course in Table 1. 𝑥-axis: size of basin of attraction for given trajectory; 𝑦-axis: number of networks observed, out of 2000.

Figure 3 

Not nested canalyzing functions with the wiring diagram in Figure 1 interpolating the time course in Table 1. 𝑥-axis: size of basin of attraction for given trajectory; 𝑦-axis: number of networks observed, out of 2000.

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.