Inferring Biologically Relevant Models: Nested Canalyzing Functions

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.


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 [21] 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 y to node x implies that changes in the concentration of y could change the concentration of x. The other part of the model is the dynamics of the network, which describes how exactly the concentration of x is affected by that of y. 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, finitespace dynamical systems. This includes Boolean networks [13], Logical models [24], Petri nets [22], and algebraic models [14]. The latter is a straightforward generalization of Boolean networks to multistate systems. Furthermore, in [25], 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 [14,10,8,25] for more details. Throughout this paper, we will be talking about gene regulatory networks, however, the methods apply for biochemical networks in general.
Suppose that the gene regulatory network that we want to infer has n genes and that we have a set D of r state transition pairs (s j , t j ), j = 1, ..., r. The input s j and the output t j are n-tuples of 0 and 1 encoding the state of genes x 1 , . . . , x n . Real time data points are not Boolean but could be discretized (and in particular, could be made Boolean) using different methods [4]. The goal is to find a model Notice that, since any function over a finite field is a polynomial, each f i is a polynomial. An algorithm that finds all models f is presented in [14]. This is done by identifying, for each gene i, the set of all possible functions for f i . This set can be represented as the coset f + I, where f is a particular such function and I ⊂ F 2 [x 1 , . . . , x n ] is the ideal of all Boolean polynomials that vanish on the input data set, that is, I = I({s 1 , . . . , s r }). The algorithm in [14] then proceeds to find a particular model from the model space f + I. The chosen model, which is the normal form of f in the ideal I, 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 [3], using the Gröbner fan of the ideal I, the authors developed a method that produces a probabilistic model using all possible normal forms. Other improvements on this algorithm can be found in [9,23].
Another approach toward improving the model selection process is by restricting the model space f + I by requiring not only that the chosen model fits the data but also satisfies some other conditions, such as its network being sparse or scale-free, the polynomials f i being monomials, the dynamics of the model having some desirable properties such as fixed points are the only limit cycles (that is, 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 f + I 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 [19], certain post classes of Boolean functions have been studied in [20], and chain functions in [5], to name few. Another class of Boolean functions, which was introduced by S. Kauffman et al. [12], 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 C. H. Waddington [26] 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 [6,17], and Boolean nested canalyzing networks to be robust and stable [11,12,17].
For the purpose of restricting the model space f + I 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 [10]. 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 F 2 . 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 f + I that fit the given data set.
In the next section we briefly recall some definitions and results from [10,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 [15].

Nested Canalyzing Functions: Background
We recall some of the definitions and the results from [10,8] that we need to make this presentation selfcontained. Throughout this paper, when we refer to a function on n variables, we mean that h depends on all n variables, that is, for i = 1, . . . , n, there exists (a 1 , . . . , a n ) ∈ F n 2 such that h(a 1 , . . . , a i−1 , a i , a i+1 , . . . , a n ) = h(a 1 , . . . , a i−1 , 1 + a i , a i+1 , . . . , a n ).

Definition 2.1.
Let h be a Boolean function on n variables, i.e., h : F n 2 → F 2 .
• The function h is a nested canalyzing function (NCF) with respect to a permutation σ on the n variables, canalyzing input value a i and canalyzed output value b i , for i = 1, . . . , n, if it can be represented in the form 1 and · · · and x σ (n−1) = a n−1 and x σ (n) = a n , b n if x σ (1) = a 1 and · · · and x σ (n) = a n . (1) • The function h is nested canalyzing if h is nested canalyzing with respect to some permutation σ , canalyzing input values a 1 , . . . , a n and canalyzed output values b 1 , . . . , b n , respectively.
Remark 2.2. The definition above has been generalized to multistate functions in [16], where it is also shown that the dynamics of these functions are similar to their Boolean counterparts. In [18], 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 Indexing monomials by the subsets of [n] := {1, . . . , n} corresponding to the variables appearing in the monomial, the elements of R can be written as The main result in [10] is the identification of the set of nested canalyzing functions in R with a subset V nc f of F 2 n 2 by imposing relations on the coordinates of its elements.
Corollary 2.5. The set of points in F 2 n 2 corresponding to the set of all nested canalyzing functions with respect to a permutation σ on [n], denoted by V nc f σ , is defined by It was shown in [8] that V nc f σ is an algebraic variety, and its ideal I(V nc f σ ) is a binomial prime ideal in the polynomial ring F 2 [{c S : S ⊆ [n]}], where F 2 is the algebraic closure of F 2 . Namely, Furthermore, the variety of all nested canalyzing functions is In the next section, we identify the set f + I with the rational points in an algebraic affine variety. This will allow us to identify all nested canalyzing functions in the model space f + I.

Nested Canalyzing Models
Recall that we are given the data set D = {(s 1 , t 1 ), . . . , (s r , t r )} ⊂ F n 2 × F n 2 . The model space could be presented by the set f + I where, f = ( f 1 , . . . , f n ) and, for i = 1, . . . , n, In particular, f i is a polynomial that interpolates the data for gene i and I is the ideal of points of {s 1 , . . . , s r }. Furthermore, the ideal I is a principal ideal in the ring R/J: (1 − (x e − s j,e ))) .
Now a polynomial h ∈ f i + I if and only if h = f i + g(x 1 , . . . , x n ) ∏ r j=1 (1− ∏ n e=1 (1− (x e − s j,e ))), for some polynomial g, say g = ∑ H⊆[n] b H ∏ i∈H x i . By expanding the right-hand side and collecting terms, we get Then ker(Φ) is the ideal of all polynomials that fit the data set D. In particular, the rational points in the variety V(ker(Φ)) is the set of all models that fit the data set D, namely f + I.
Since the ideal of all NCFs is I(V nc f ), the following corollary is straightforward.

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, i.e., a square matrix of dimension n, describing the influence relationships among the n genes in the network. For each variable x i , a  (10) For each variable x i do 1. Compute f i as in (5) 2. Let g = f i + h * p; its coefficients are the same as W S above 3. Compute a Gröbner basis G for the ideal generated by the coefficients of g − q using any elimination order to eliminate all b S from G 4. Concatenate generators of G and I(V nc f ) 5. Compute the primary decomposition of G + I(V nc f ) to obtain necessary and sufficient conditions on the coefficients of all NCFs fitting the data set D End An implementation of this algorithm is available as a Singular library [7,2].

Application: Inferring the Cell Cycle Network in Budding Yeast
Li et al. [15] constructed a Boolean threshold model of the cell cycle in budding yeast. The network of the model has the key known regulators of the cell cycle process and the known interactions among these regulators in the literature. The Boolean function at each node, however, is a threshold function, which is completely determined by the numbers of active activators and active inhibitors, and is not necessarily biologically motivated. However, 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 above, 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 network. It consists of 11 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 [15] 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 [15]. 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. [15] 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 G 1 , S (Synthesis), G 2 , and M (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

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.

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.

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. [15]. We then showed that the dynamics of almost any of these models is strikingly similar to that  Figure 1 interpolating the time course in Table 1. x-axis: size of basin of attraction for given trajectory, y-axis: number of networks observed, out of 2000.  Figure 1 interpolating the time course in Table 1. x-axis: size of basin of attraction for given trajectory, y-axis: number of networks observed, out of 2000.
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.