# Inferring Biologically Relevant Models: Nested Canalyzing Functions

Authors: Franziska Hinkelmann; Abdul Salam Jarrah
Journal: ISRN Biomathematics
 (2012)

---

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


---

## Body

## 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 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 istime-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 hasn genes and that we have a set D of r state transition pairs (sj,tj), j=1,…,r. The input sj and the output tj are n-tuples of 0 and 1 encoding the state of genes x1,…,xn. 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:(1)f=(f1,f2,…,fn):F2n⟶F2n,
such that for j=1,…,r:(2)f(sj)=(f1(sj),…,fn(sj))=tj.
Notice that, since any function over a finite field is a polynomial, each fi is a polynomial. An algorithm that finds all models f is presented in [5]. This is done by identifying, for each gene i, the set of all possible functions for fi. This set can be represented as the coset f+I, where f is a particular such function, and I⊂𝔽2[x1,…,xn] is the ideal of all Boolean polynomials that vanish on the input data set, that is, I=𝕀({s1,…,sr}). The algorithm in [5] 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 [10], 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 [11, 12].Another approach toward improving the model selection process is by restricting the model spacef+I 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 fi 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 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 [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 spacef+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 [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 𝔽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 [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 n variables, we mean that h depends on all n variables, that is, for i=1,…,n, there exists (a1,…,an)∈𝔽2n such that h(a1,…,ai-1,ai,ai+1,…,an)≠h(a1,…,ai-1,1+ai,ai+1,…,an).Definition 1.
Leth be a Boolean function on n variables, that is, h:𝔽2n→𝔽2:(i)
the functionh is a nested canalyzing function (NCF) with respect to a permutation σ on the n variables, canalyzing input value ai and canalyzed output value bi, for i=1,…,n, if it can be represented in the form:
(3)h(x1,x2,…,xn)={b1ifxσ(1)=a1,b2ifxσ(1)≠a1,xσ(2)=a2,b3ifxσ(1)≠a1,xσ(2)≠a2,xσ(3)=a3,⋮⋮bnifxσ(1)≠a1,…,xσ(n-1)≠an-1,xσ(n)=an,bn¯ifxσ(1)≠a1,…,xσ(n)≠an;(ii)
the functionh is nested canalyzing if h is nested canalyzing with respect to some permutation σ, canalyzing input values a1,…,an and canalyzed output values b1,…,bn, 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 ringR=𝔽2[x1,…,xn]/J, where J=〈xi2-xi:1≤i≤n〉. 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(4)R={∑S⊆[n]cS∏i∈Sxi:cS∈F2}.As a vector space over𝔽2, R is isomorphic to 𝔽22n via the correspondence:(5)R∋∑S⊆[n]cS∏i∈Sxi⟷(c∅,…,c[n])∈F22n.The main result in [7] is the identification of the set of nested canalyzing functions in R with a subset Vncf of 𝔽22n by imposing relations on the coordinates of its elements.Definition 3.
Letσ be a permutation of the elements of the set [n]. We define a new order relation <σ on the elements of [n] as follows: σ(i)<σσ(j) if and only if i<j. Let rSσ be the maximum element of a nonempty subset S of [n] with respect to the order relation <σ. For any nonempty subset S of [n], the completion of S with respect to the permutationσ, denoted by [rSσ], is the set [rSσ]={σ(1),σ(2),…,σ(rSσ)}.
Note that, ifσ is the identity permutation, then the completion is [rS]:= {1,2,…,rS}, where rS is the largest element of S.Theorem 4.
Leth∈R, and let σ be a permutation of the set [n]. The polynomial h is nested canalyzing with respect to σ, input value ai and corresponding output value bi, for i=1,…,n, if and only if c[n]=1 and, for any proper subset S⊆[n]:
(6)cS=c[rSσ]∏σ(i)∈[rSσ]∖Sc[n]∖{σ(i)}.Corollary 5.
The set of points in𝔽22n corresponding to the set of all nested canalyzing functions with respect to a permutation σ on [n], denoted by Vσncf, is defined by
(7)Vσncf={∏σ(i)∈[rSσ]∖S(c∅,…,c[n])∈F22n:c[n]=1,cS=c[rSσ]×∏σ(i)∈[rSσ]∖Sc[n]∖{σ(i)},forS⊆[n]}.It was shown in [8] that Vσncf is an algebraic variety, and its ideal 𝕀(Vσncf) is a binomial prime ideal in the polynomial ring 𝔽¯2[{cS:S⊆[n]}], where 𝔽¯2 is the algebraic closure of 𝔽2. Namely,(8)Iσ=I(Vσncf)=〈c[n]-1,cS-c[rSσ]∏σ(i)∈[rSσ]∖Sc[n]∖{σ(i)}:S⊂[n]〉.Furthermore, the variety of all nested canalyzing functions is(9)Vncf=⋃σVσncf,and its ideal is(10)I(Vncf)=⋂σIσ.In the next section, we identify the setf+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.


## 3. Nested Canalyzing Models

Recall that we are given the data setD={(s1,t1),…,(sr,tr)}⊂𝔽2n×𝔽2n. The model space could be presented by the set f+I, where f=(f1,…,fn) and, for i=1,…,n,(11)fi(x1,…,xn)=∑j=1rtj,i∏e=1n(1-(xe-sj,e)).
In particular, fi is a polynomial that interpolates the data for gene i and I is the ideal of points of {s1,…,sr}. Furthermore, the ideal I is a principal ideal in the ring R/J:(12)I=I({s1,…,sr})=⋂j=1rI({sj})=⋂j=1r〈x1-sj,1,…,xn-sj,n〉=⋂j=1r〈1-∏e=1n(1-(xe-sj,e))〉=〈∏j=1r(1-∏e=1n(1-(xe-sj,e)))〉.
Now a polynomial h∈fi+I if and only if h=fi+g(x1,…,xn)∏j=1r(1-∏e=1n(1-(xe-sj,e))), for some polynomial g, say g=∑H⊆[n]bH∏i∈Hxi. By expanding the right-hand side and collecting terms, we get that h=∑S⊆[n]WS(bH,sj,tj)∏l∈Sxl, where for S⊆[n], the coefficient WS(bH,sj,tj) is determined by bH,sj,tj for all H⊆[n] and j=1,…,r.The proof of the following theorem follows directly from Theorem  2.4.2 in [24].Theorem 6.
Consider the ring homomorphism:(13)Φ:F¯2[{cS:S⊆[n]}]⟶F¯2[{bH:H⊆[n]}]
given by, for S⊆[n],
(14)cS⟼WS(bH,sj,tj).
Then ker(Φ) is the ideal of all polynomials that fit the data set D. In particular, the rational points in the variety 𝕍(ker(Φ)) is the set of all models that fit the data set D, namely f+I.Since the ideal of all NCFs is𝕀(Vncf), the following corollary is straightforward.Corollary 7.
The ideal of all nested canalyzing functions that fit the data setD is 𝕀(Vncf)+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 are25=32 different monomials in 5 variables, and hence 232=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 dimensionn, describing the influence relationships among the n genes in the network. For each variable xi, a table consisting of the rows (sj,i1,…,sj,is,tj,i), j={1,…,r}, where i1,…,is are the indices of the genes that affect xi, 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 〈{cS-WS:S⊆[n]}〉 intersected with the ring F¯2[{cS:S⊆[n]}] ([24], Theorem  2.4.2) Using a similar notation as above, the algorithm is outlined as follows:use ring𝔽2[{x1,…,xn,bS,cS:S⊆[n]}];define𝕀(Vncf) as ideal in 𝔽2[{cS:S⊆[n]}];defineh=∑H⊆[n]bH∏i∈Hxi;defineq=∑H⊆[n]cH∏i∈Hxi;compute the polynomialp that generates I as in (12);for each variablexi do
(1)
computefi as in (11);(2)
letg=fi+h*p; its coefficients are the same as WS above;(3)
compute a Gröbner basisG for the ideal generated by the coefficients of g-q using any elimination order to eliminate all bS from G;(4)
concatenate generators ofG and 𝕀(Vncf);(5)
compute the primary decomposition ofG+𝕀(Vncf) 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 [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 yeastSaccharomyces 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 Table1, 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.
TimeCln3MBFSBFCln1,2Cdh1Swi5Cdc14,20Clb5,6Sic1Clb1,2Mcm1/SFF1100010001002011010001003011110001004011100000005011100010006011100010117000100110118000001100119000001101111000000110101110000111010012000011001001300001000100We 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 (Figure1) 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 G1, S (synthesis), G2, 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 Table2, 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 proteini, we list the number of inputs, the number of possible Boolean functions (the cardinality of fi+I), the number of nested canalyzing functions with the given number of inputs, and finally the number of nested canalyzing functions in the model space fi+I.
Protein (i)Inputsfi+INCFsNCFs infi+ICln31121MBF38642SBF38642Cln1,21121Cdh14204873612Swi54204873614Cdc20 and Cdc1438644Clb5,638643Sic1522410,634336Clb1,2522410,63461Mcm1/SFF38642


### 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 Table1 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 Figures2 and 3.Figure 2
Nested canalyzing functions with the wiring diagram in Figure1 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 3
Not nested canalyzing functions with the wiring diagram in Figure1 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.


## 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 Table1 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 Figures2 and 3.Figure 2
Nested canalyzing functions with the wiring diagram in Figure1 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 3
Not nested canalyzing functions with the wiring diagram in Figure1 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.


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



---

*Source: 613174-2012-06-12.xml*