Representing perturbed dynamics in biological network models

Representing perturbed dynamics in biological network models Gautier Stoll1,3 ,∗ Jacques Rougemont3 ,† and Felix Naef1,2,3‡ 1 arXiv:q-bio/0702059v1 [q-bio.MN] 28 Feb 2007 3 NCCR Molecular Oncology, ch. des Boveresses 155, 1066 Epalinges, Switzerland School of Life Sciences, ISREC, Ecole polytechnique Fédérale de Lausanne 1015 Lausanne, Switzerland and 3 Swiss Institute of Bioinformatics, Quartier Sorge-Genopode, 1015 Lausanne, Switzerland We study the dynamics of gene activities in relatively small size biological networks (up to a few tens of nodes), e.g. the activities of cell-cycle proteins during the mitotic cell-cycle progression. Using the framework of deterministic discrete dynamical models, we characterize the dynamical modifications in response to structural perturbations in the network connectivities. In particular, we focus on how perturbations affect the set of fixed points and sizes of the basins of attraction. Our approach uses two analytical measures: the basin entropy H and the perturbation size ∆, a quantity that reflects the distance between the set of fixed points of the perturbed network to that of the unperturbed network. Applying our approach to the yeast-cell cycle network introduced by Li et al. provides a low dimensional and informative fingerprint of network behavior under large classes of perturbations. We identify interactions that are crucial for proper network function, and also pinpoints functionally redundant network connections. Selected perturbations exemplify the breadth of dynamical responses in this cell-cycle model. PACS numbers: Keywords: I. INTRODUCTION Recent experimental developments in the fields of genomics, e.g. whole genome DNA sequencing or proteomics, are opening possibilities for systems level studies in biology [1, 2, 3, 4]. In particular, the notion that biological functions may rely on a large number of interconnected variables (for example genes) working in concert has stimulated general theoretical interest about properties of biological networks [5]. Studies of the statistical properties of large (typically thousands of nodes) biological networks have identified a number of functional building block, termed network motifs, that occur more frequently than random [7]. These findings support the idea that some systems are designed around a modular architecture, in which autonomous modules are wired together to generate versatile biological functions [1, 4, 8, 25]. While structural (or topological) properties are key for network characterization, functional properties are ultimately encoded in dynamical, or timedependent changes in the state variables of the nodes. The sizes of systems that can be modeled dynamically are typically much smaller (10-100 nodes). One common modeling approach, for example for the yeast cell-cycle [10], is to simulate the nonlinear system of chemical rate equations describing the putative biochemical processes. Modeling approaches have been applied to a number of systems, including the cell-cycle [10, 11], the lambdaphage switch in E. coli [9]. Although these models provide a detailed description, this approach suffers from the ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] caveat that most parameters are currently not accessible experimentally. In addition, the number of parameters is typically about five per reaction, resulting in a prohibitively large parameter space. This last point makes it difficult to grasp the full solution space of the model. Recent approaches based on sampling the parameter space in optimal regions have been developed [19]. At the opposite end of model complexity, dynamical rules based on boolean state variables have been useful for studying more global dynamical properties of topological classes of networks [20, 21]. In addition, boolean models have been successfully applied to the yeast cell-cycle [22, 26] and the body patterning in drosophila embryos [23, 24]. In this study, we develop a systematic approach to describe how the dynamical landscape of small (less than about 50 nodes) boolean networks is affected by perturbations in the network connectivity. In particular, we consider the basin entropy H, a quantity that considers the size distribution of the basins of attraction. We complement entropy with a measure of distance between the stable fixed points of a perturbed network and those in the unperturbed network. This combination gives a low-dimensional and compact representation of the patterns induced by a large number of perturbations. We illustrate our methods using the yeast cell-cycle network introduced in [22], and discuss examples of structural perturbations producing a range of modified basins of attraction. II. DEFINITIONS Following [22] a network of N nodes can be represented by a N × N adjacency matrix A, in which an activating link between node i and node j is represented by Aij = 1 and an inhibiting link by Aij = −1. The possibility 2 of self-inhibitory (or activating links) Aii = ±1 is not excluded. In the Boolean approximation, each node has two possible states, so that the global state of all nodes can be represented by a vector S, with Si = 1 when the node i is on and Si = 0 if the node is off. The full phase space containing 2N states is denoted by Λ. A. Boolean dynamics A simple dynamical rule that characterizes the temporal evolution of the state variable can be defined following [22], which is closely related to update rules applied in perceptron models. If the network is in the state S(t) at time t, the state at the next time-step S(t + 1) is given by: P 1 if Pj Aij Sj (t) > 0 Si (t + 1) = Si (t) if Pj Aij Sj (t) = 0  0 if j Aij Sj (t) < 0   (1) For a given network, we apply this rule to every possible initial condition in Λ. This defines orbits (trajectories) that must end in a limit cycle (periodic attractor) since we are dealing with a dynamical system on a finite space. A fixed point is a cycle of length one. Accordingly, Λ can be decomposed into a disjoint union SK of K basins of attraction Bk of size dk : Λ = k=1 Bk . In a biological network, the attractors correspond to functional endpoints, and it is important that the states in the attractors are consistent with observed data. For example, by far the largest endpoint in the cell-cycle network of Li et al. (see appendix) corresponds to the stationary G1 phase in the cycle. Other systems are more switch-like, for instance in signal transduction, where a cell might change its state from growth to differentiation according to an external trigger. To characterize these attractors, we introduce the following definitions: • We compute the number of attractors K: an attractor is a limit cycle or a fixed point. An attractor A has a basin of attraction B which is the set of all initial conditions whose orbit converges to A. • The basin entropy H is defined as follows: let pk = 2−N dk be the probability that an initial state belongs to basin Bk . Then, the entropy reads H := − K X pk log (pk ) (2) k=1 H is maximum (H = log(K)) if each state is its own basin of size one, and minimum (H = 0) when there is one single basin. H is a natural measure for characterizing basin structures [27]. Because it takes into account the relative basin sizes, it is quite insensitive to appearance of small and biologically irrelevant basins. • The perturbation size ∆ measures the distance between attractors of a perturbed and a reference network: from every initial conditions, the Hamming distance between the fixed points is computed, and the average over all initial conditions is taken. More precisely, if FPG (S) is the fixed point of the trajectory starting at S and generated by the network G, then 1 X ∆G,G′ := N (3) HAM(FPG (S), FPG′ (S)) 2 S where HAM(·, ·) is the Hamming distance between two boolean states, namely 1 X |Si − Ti | (4) HAM(S, T) := N i The value ∆ has the following interpretation: it is the average probability (taken over all nodes) that, for a random initial condition, the final state of a node differs. In this study, the reference network G will be the cell-cycle network of Li et al., which has one very a large basin of attraction and several smaller ones. If some trajectories in the perturbed networks G′ end in a limit cycle, ∆ is defined as the average of the Hamming distance along the cycle. B. Network models and perturbations Our goal is to assess how network dynamics is affected by several types of perturbations. We consider two classes: one which randomizes the adjacency matrix while keeping a number of topological characteristics from the original network invariant. The second class mimics biological perturbations, as would occur for example through mutations in the interaction partners that constitute the network links.The two classes are defined as follows: • Shuffle (class I): all activating and inhibiting arrows are cut in half and re-wired randomly. This ensure that the connectivity at each node is conserved. As compared to the Li et al. [22] study, we generate random networks that are more constrained, since the connectivity at each node is forced to remain unchanged after randomization. Such perturbations are applied in the studies of network motifs [4, 7]. • Remove (class II): the arrows are simply suppressed. We extend this class of perturbations beyond single link removal. III. RESULTS We study the yeast cell-cycle network of Li et al. [22] (the Yeast cell-cycle network or YCC ), in which a 3 boolean model reproducing the different phases of the cycle is constructed (see appendix). This model has a main fixed point attracting 86% of the intial conditions. Biologically this state corresponds to the G1 stationary phase of the cell-cycle, as reflected by the activities of the respective nodes. Using computer simulations, the authors further showed that the cell-cycle dynamics had certain robustness properties when challenged with perturbations. In particular, it was shown that in a majority of cases, removal of one link or addition of a link at random did not change much the size of the largest basin of attraction. Finally, the studied network had unusual trajectory channeling properties, when compared to random networks with equal number of nodes and links. Here we extend the characterization of this model by introducing a combination of measures to characterize the structure of basins of attraction as they are modified by structural perturbations. In particular we investigate the consequences of combined mutations and show that they can lead to cancellation effect. A. Study of shuffled networks (Class I perturbations) This type of perturbation allows to study the dynamical characteristics of a biological network in comparison with random networks belonging to a topological class. Figure 1 shows the Number of attractors (K) and the Entropy (H) of the YCC and randomly shuffled (Class I) versions thereof. The location of the reference network in the H − K plane respective to the scatter of the perturbed networks allows us to asses how typical a network behaves with respect to a class. Accordingly, the YCC is atypical, as seen by its marginal location in the lower left corner. Indeed, this network has lower entropy and fewer basins than most networks, consistent with [22]. B. Study of mutated networks (Class II perturbations) The previous discussion shows how entropy characterizes the system of attractors. However, H contains only information about the relative weights of the attractors, irrespective of their biological relevance. For example a perturbation can decrease the entropy while shifting the fixed point away from from that in the unperturbed, biologically relevant state. For this reason we introduced a second quantity, ∆ (Equation 3), a probabilistic measure of the change in the fixed point after perturbation. Therefore, ∆ reflects the change in the biological relevance of the basin structure. We first repeat Figure 1 for class II perturbations which shows that networks with few perturbations cluster around the wild-type model (Figure 2A), while the sread for networks with four perturbations resembles the FIG. 1: Entropy vs. number of attractors after class I perturbation (shuffled arrows). The range of possible H values is indicated by the dashed gray lines. The open red circle represents the reference network, the other points show the perturbed networks. shuffled models (Figure 1). Turning to the measure of ∆, we find that ∆-distribution (Figure 2B) is bimodal, showing two distinct populations of perturbations: (∆ . 0.2 and ∆ & 0.2). In the second case, the perturbed model does not reproduce the biologically correct cell-cycle progression. But if ∆ is small, then the system of attractors of the perturbed network is still consistent with the biology and entropy allows to discriminate between networks with a larger or smaller main basin of attraction. For this reason, the entropy and ∆ are complementary for describing the dynamical landscape (Figure 2C). The two different modes in the ∆-histogram are clearly reflected on this 2D representation. Noticeably, the ∆ values span a broad range for any number of removed arrows, on the other hand higher entropies are more frequent for larger number (> 2) of removed arrows. Qualitatively, the spread of points in the H − ∆ plane conveys a measure of robustness. Accordingly, the ∆ measure appears more fragile than the entropy property, especially when few arrows are removed. We now interpret the different locations in the H − ∆ plane: 1. If ∆ is large (∆ & 0.2), the model does have attractor states which coincide with the gene activities of the different cell-cycles phases. Such perturbations are specially interesting if the number of removed arrows is small (dark colors). Such links are then essential for the model, as their removal disrupts the cell-cycle very 4 C. Examples of mutations In the first example (Table II), the dynamics has a large main basin of attraction like in the unperturbed model (Table I). However, the fixed point is significantly different from wild-type as the system is blocked in the a state of the M-phase and cannot finish properly the cellcycle (see Appendix for the recapitulation of the wildtype mode from [22]). Cln3 MBF SBF Cln1,2 Cdh1 Swi5 Cdc20,14 Clb5,6 Sic1 Clb1,2 Mcm1 % 0 0 0 0 1 0 0 0 1 0 0 0.8613 0 0 1 1 0 0 0 0 0 0 0 0.0737 0 1 0 0 1 0 0 0 1 0 0 0.0532 0 0 0 0 0 0 0 0 1 0 0 0.0043 0 0 0 0 0 0 0 0 0 0 0 0.0034 0 1 0 0 0 0 0 0 1 0 0 0.0034 0 0 0 0 1 0 0 0 0 0 0 0.0004 TABLE I: Basins of attraction with their respective probabilities in (%) for the original YCC network. The largest basin ends at the G1 stationary state. Entropy H = 0.543, Number of attractors K = 7. In the second example (Table III), the dynamics has the same main fixed point as the wild-type, but with a smaller basin of attraction, while the second biggest has grown. Therefore the removed connection SBF → Cln1,2 contributes to the ability of the main fixed point to funnel trajectories. The third example (Table IV) is a model with four removed arrows which has the same main fixed point with a slightly higher probability. Also, the second largest fixed point is same as in the wild-type model. This indicates that the effect of some mutations can be canceled by further mutations. While such cases exist, we found that networks with several removed links that preserving the unperturbed cell-cycle behavior are rare. IV. FIG. 2: Entropy, number of attractors and ∆ after class II perturbation (removed arrows). Colors represent different number of removed arrows: black for one removed arrows, red for 2, green for 3, turquoise for 4 and yellow for more than 4. A: same figure as for the class I perturbation, the range of possible H values is indicated by the dashed gray lines, the open blue circles represent the reference network. B: Distribution of ∆. C: ∆ vs. H plot, the dashed gray line represents the entropy of the reference network. efficiently. 2. If ∆ is small and the entropy increases, the probability that the dynamics ends in the reference attractor decreases demonstrating that the removed arrows contributed to the channeling properties of the system. 3. If ∆ is small and the entropy decreases, the main attractor of the perturbed network has a stronger attraction property. Some of these networks could be considered as alternative cell-cycle models. We illustrate these three regimes by examples: CONCLUSION We have proposed a systematic approach for studying the dynamical attractor landscape of biological networks, and their response to structural perturbations. In particular, we introduced a low dimensional representation of the system of attractors, the entropy, and a probabilistic measure in the perturbation size ∆. This enabled us to study the global characteristics of network perturbation in a compact and visually effective form. In a biological context, this can provide hints to elucidate the dynamical role of specific network links. Alternatively, the function Cln3 MBF SBF Cln1,2 Cdh1 Swi5 Cdc20,14 Clb5,6 Sic1 Clb1,2 Mcm1 % 0 0 0 0 0 1 1 0 1 1 1 0.880 0 0 0 0 1 0 0 0 1 0 0 0.054 0 0 1 1 0 0 0 0 0 0 0 0.027 0 1 0 0 1 0 0 0 1 0 0 0.015 0 0 0 0 1 1 1 0 1 1 1 0.010 0 0 0 0 0 0 0 0 1 0 0 0.004 0 0 0 0 0 0 0 0 0 0 0 0.003 0 1 0 0 0 0 0 0 1 0 0 0.003 0 0 0 0 1 0 0 0 0 0 0 0.000 TABLE II: Basins of attraction with their respective probabilities, when (Cdc20,Cdc14) → Clb1,2 and Sic1 → Clb1,2 are removed. Entropy = 0.549, Number of attractors = 9, ∆ = 0.41. 5 Cln3 MBF SBF Cln1,2 Cdh1 Swi5 Cdc20,14 Clb5,6 Sic1 Clb1,2 Mcm1 % 0 0 0 0 1 0 0 0 1 0 0 0.6669 0 1 1 0 1 0 0 0 1 0 0 0.1762 0 0 1 0 1 0 0 0 1 0 0 0.0654 0 1 0 0 1 0 0 0 1 0 0 0.0532 0 1 1 0 0 0 0 0 1 0 0 0.0180 0 0 1 0 0 0 0 0 1 0 0 0.0043 0 0 0 0 0 0 0 0 1 0 0 0.0043 0 0 0 0 0 0 0 0 0 0 0 0.0034 0 0 1 0 0 0 0 0 0 0 0 0.0034 0 1 0 0 0 0 0 0 1 0 0 0.0034 0 0 0 0 1 0 0 0 0 0 0 0.0004 0 0 1 0 1 0 0 0 0 0 0 0.0004 work was presented. The simulations were performed on an Itanium2 cluster from HP/Intel at the Vital-IT facilities. FN ad GS acknowledge funding from the NCCR Molecular Oncology program and NIH administrative supplement to parent grant GM54339. APPENDIX A: THE YEAST CELL-CYCLE NETWORK OF LI ET AL. TABLE III: Basins of attraction with their respective probabilities, when SBF → Cln1,2 is removed. Entropy H = 1.096, Number of attractors K = 12, ∆ = 0.05. Cln3 MBF SBF Cln1,2 Cdh1 Swi5 Cdc20,14 Clb5,6 Sic1 Clb1,2 Mcm1 % 0 0 0 0 1 0 0 0 1 0 0 0.8793 0 0 1 1 0 0 0 0 0 0 0 0.0507 0 1 0 0 1 0 0 0 1 0 0 0.0356 0 0 0 0 0 0 0 0 1 0 0 0.0268 0 1 0 0 0 0 0 0 1 0 0 0.0034 0 0 0 0 0 0 0 0 0 0 0 0.0034 0 0 0 0 1 0 0 0 0 0 0 0.0004 TABLE IV: Basins of attraction with their respective probabilities, when (Cdc20,Cdc14) → Clb1,2, Clb1,2 → Mcm1, Clb1,2 → Cdh1 and Clb1,2 → Swi5 are removed. Entropy H = 0.523, Number of attractors K = 7, ∆ = 0.025. of new and yet unobserved links can be predicted as in [26], and imperfect starting models can be improved. We applied this method to a model of the yeast cellcycle by Li et al. Using the measures introduced, we have generalized the dynamical characterization of the model using a broad range of perturbations. This has enabled us to emphasize the breadth of dynamical behavior (Figure 2) induced by only few mutated links. Interestingly, we observed (Figure 2C) that the structure of the system of attractors (H) behaves quite robustly compared to the modification in the final states (∆), especially when the number of removed links is small (< 3). We illustrated through examples the consequences of removing individual or groups of links. Interestingly it was possible to remove up to four links while not affecting the basin structure significantly. Tracking the dynamical changes in the activity levels of proteins in a network is a very high-dimensional problem. It therefore important to be have few informative variables which allow one to efficiently assess a large number of perturbed models at once. We believe that basin entropy and distance to a reference attractor are well suited for this purpose. The following two tables are recapitulated from [22]. + + + + + + + + + + + + 1→2 1→3 2→8 3→4 6→9 7→5 7→6 7→9 8→10 8→11 10→7 10→11 + + + − − − − − − − − − 11→6 11→7 11→10 4→9 4→5 5→10 7→8 7→10 8→5 8→9 9→8 9→10 − − − − − − − − − − 10→2 10→3 10→5 10→6 10→9 1→1 4→4 6→6 7→7 11→11 TABLE V: Adjacency matrix of the Yeast cell-cycle network. The numbers refer to the ordering of the nodes as used in Tables I-IV,VI. + (respectively −) represent activating (respectively repressing) links. t 1 2 3 4 5 6 7 8 9 10 11 12 13 Cln3 MBF SBF Cln1,2 Cdh1 Swi5 C20,14 Clb5,6 Sic1 Clb1,2 Mcm1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 1 1 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 Phase START G1 G1 G1 S G2 M M M M M G1 G1* TABLE VI: This table represents the discrete time evolution of the boolean states of the YCC network as it traverses the different cell-cycle phases. Cdc20.14 has been abbreviated C20,14; G1* indicates the stationary G1 phase. Acknowledgments We thank the organizers of the CompBioNets ’04 conference (Recife, Brazil) at which an initial version of this [1] Hartwell, L. H., Hopfield, J. J., Leibler, S., and Murray, A. W. (1999). From molecular to modular cell biology, Nature 402, C47-52. [2] Alm, E., and Arkin, A. P. (2003). Biological networks, Curr Opin Struct Biol 13, 193-202. [3] Oltvai, Z. N., and Barabasi, A. L. (2002). Systems biol- 6 ogy. Life’s complexity pyramid, Science 298, 763-4. [4] Alon, U. (2003). Biological networks: the tinkerer as an engineer, Science 301, 1866-7. [5] Barabasi, A. L., and Oltvai, Z. N. (2004). Network biology: understanding the cell’s functional organization, Nat Rev Genet 5, 101-13. [6] Barabasi, A. L., and Albert, R. (1999). Emergence of scaling in random networks, Science 286, 509-12. [7] Milo, R., Shen-Orr, S., Itzkovitz, S., Kashtan, N., Chklovskii, D., and Alon, U. (2002). Network motifs: simple building blocks of complex networks, Science 298, 824-7. [8] Ihmels, J., Friedlander, G., Bergmann, S., Sarig, O., Ziv, Y., and Barkai, N. (2002). Revealing modular organization in the yeast transcriptional network, Nat Genet 31, 370-7. [9] Arkin, A., Ross, J., and McAdams, H. H. (1998). Stochastic kinetic analysis of developmental pathway bifurcation in phage lambda-infected Escherichia coli cells, Genetics 149, 1633-48. [10] Novak, B., Csikasz-Nagy, A., Gyorffy, B., Chen, K., and Tyson, J. J. (1998). Mathematical model of the fission yeast cell cycle with checkpoint controls at the G1/S, G2/M and metaphase/anaphase transitions, Biophys Chem 72, 185-200. [11] Cross, F. R., Archambault, V., Miller, M., and Klovstad, M. (2002). Testing a mathematical model of the yeast cell cycle, Mol Biol Cell 13, 52-70. [12] Barkai, N., and Leibler, S. (2000). Circadian clocks limited by noise, Nature 403, 267-8. [13] Vilar, J. M., Kueh, H. Y., Barkai, N., and Leibler, S. (2002). Mechanisms of noise-resistance in genetic oscillators, Proc Natl Acad Sci U S A 99, 5988-92. Epub 2002 Apr 23. [14] Leloup, J. C., and Goldbeter, A. (2003). Toward a detailed computational model for the mammalian circadian clock, Proc Natl Acad Sci U S A 100, 7051-6. Epub 2003 May 29. [15] von Dassow, G., Meir, E., Munro, E. M., and Odell, G. M. (2000). The segment polarity network is a robust developmental module, Nature 406, 188-92. [16] Mello, B. A., and Tu, Y. (2003). Quantitative modeling of sensitivity in bacterial chemotaxis: the role of coupling among different chemoreceptor species, Proc Natl Acad Sci U S A 100, 8223-8. Epub 2003 Jun 25. [17] Mello, B. A., and Tu, Y. (2003). Perfect and near-perfect adaptation in a model of bacterial chemotaxis, Biophys J 84, 2943-56. [18] Hoffmann, A., Levchenko, A., Scott, M. L., and Baltimore, D. (2002). The IkappaB-NF-kappaB signaling module: temporal control and selective gene activation, Science 298, 1241-5. [19] Brown, K. S., and Sethna, J. P. (2003). Statistical mechanical approaches to models with many poorly known parameters, Phys Rev E Stat Nonlin Soft Matter Phys 68, 021904. Epub 2003 Aug 12. [20] Aldana, A., Coppersmith, S., Kadanoff, L.P. (2002), Boolean dynamics with random couplings (http://arxiv.org/abs/nlin/0204062). [21] Kauffman, S., Peterson, C., Samuelsson, B., and Troein, C. (2003). Random Boolean network models and the yeast transcriptional network, Proc Natl Acad Sci U S A 100, 14796-9. Epub 2003 Dec 1. [22] Li, F., Long, T., Lu, Y., Ouyang, Q., and Tang, C. [23] [24] [25] [26] [27] (2004). The yeast cell-cycle network is robustly designed, Proc Natl Acad Sci U S A 101, 4781-6. Epub 2004 Mar 22. Albert, R., and Othmer, H. G. (2003). The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in Drosophila melanogaster, J Theor Biol 223, 1-18. Chaves, M., Albert, R., and Sontag, E. D. (2005). Robustness and fragility of Boolean models for genetic regulatory networks, J Theor Biol 235, 431-449. Vazquez, A., Dobrin, R., Sergi, D., Eckmann, J.-P., Oltvai, Z. N., and Barabasi, A.-L. (2004). The topological relationship between the large-scale attributes and local interaction patterns of complex networks, Proc Natl Acad Sci U S A 101, 17940-5. Stoll G., Rougemont J. and Naef F. (2006). Few crucial links assure checkpoint efficiency in the yeast cell-cycle network, Bioinformatics 22(20):2539-46. Cover, T. M. and Thomas, J.A. (1991). Elements of information theory, Wiley.