Functional Dynamical Structures in Complex Systems: an Information-Theoretic Approach

TL;DR

This paper introduces an information-theoretic method to analyze complex systems' dynamics without prior knowledge of their network topology, identifying influential groups and functional structures.

Contribution

It develops the Dynamical Cluster Index and D-index, novel measures for detecting relevant and influential subsets in complex systems using only dynamical data.

Findings

Effective detection of relevant subsets in artificial systems

Identification of influential groups in natural systems

Method does not require network topology knowledge

Abstract

Understanding the dynamical behavior of complex systems is of exceptional relevance in everyday life, from biology to economy. In order to describe the dynamical organization of complex systems, existing methods require the knowledge of the network topology. By contrast, in this thesis we develop a new method based on Information Theory which does not require any topological knowledge. We introduce the Dynamical Cluster Index to detect those groups of system elements which have strong mutual interactions, named as Relevant Subsets. Among them, we identify those which exchange most information with the rest of the system, thus being the most influential for its dynamics. In order to detect such Functional Dynamical Structures, we introduce another information theoretic measure, called D-index. The experimental results make us confident that our method can be effectively used to study both artificial and natural complex systems.

Figures

Tables

Table 1. Table 4.1: The update rules of the boolean networks discussed on the text. Random( 0.5 0.5 0.5 ) denotes a Bernoulli distribution with probability 0.5 0.5 0.5 .

Node	Node Rule
	Case 1	Case 2	Case 3	Case 4	Case 5
N01	Random($0.5$)	Random($0.5$)	Random($0.5$)	Random($0.5$)	Random($0.5$)
N02	Random($0.5$)	Random($0.5$)	Random($0.5$)	Random($0.5$)	Random($0.5$)
N03	(N04 $\oplus$ N05)	(N04 $\oplus$ N05)	N10$\wedge$(N04 $\oplus$ N05)	(N04 $\oplus$ N05)	(N04 $\oplus$ N05)
N04	(N03 $\oplus$ N05)	(N03 $\oplus$ N05)	(N03 $\oplus$ N05)	(N03 $\oplus$ N05)	(N03 $\oplus$ N05)
N05	(N03 $\oplus$ N04)	(N03 $\oplus$ N04)	(N03 $\oplus$ N04)	(N03 $\oplus$ N04)	(N03 $\oplus$ N04)
N06	Random($0.5$)	Random($0.5$)	Random($0.5$)	(N05 $\oplus$ N08)	(N05 $\oplus$ N08)
N07	Random($0.5$)	Random($0.5$)	Random($0.5$)	(N07$+$N08$+$N09$+$N10) $\ge 2$	$\neg$(N05 $\oplus$ N08)
N08	(N09 $\oplus$ N10)	N05$\wedge$(N09 $\oplus$ N10)	N05$\wedge$(N09 $\oplus$ N10)	N03 $\oplus$ N05	N09 $\oplus$ N10
N09	(N08 $\oplus$ N10)	(N08 $\oplus$ N10)	(N08 $\oplus$ N10)	(N04$+$N05$+$N07$+$N08)$\le 2$	(N08 $\oplus$ N10)
N10	(N08 $\oplus$ N09)	(N08 $\oplus$ N09)	(N08 $\oplus$ N09)	N06$\wedge$(N05 $\oplus$ N09)	(N08 $\oplus$ N09)
N11	Random($0.5$)	Random($0.5$)	Random($0.5$)	Random($0.5$)	Random($0.5$)
N12	Random($0.5$)	Random($0.5$)	Random($0.5$)	Random($0.5$)	Random($0.5$)

Table 2. Table 4.2: Adapted from  [ 6 ] , the boolean regulatory network of mammalian cell cycle network and a short description of each node of the system.

Product	Logical Rule leading to an activity of product	Legend
$CycD$	$CycD$	Cyclin D
$Rb$	$\left(\overline{CycD}\wedge \overline{CycE}\wedge \overline{CycA}\wedge \overline{CycB}\right)\vee \left(p27\wedge \overline{CycD}\wedge \overline{CycB}\right)$	Retinoblastoma Protein
$E2F$	$\left(\overline{rB}\wedge \overline{CycA}\wedge \overline{CycB}\right)\vee \left(p27\wedge \overline{rB}\wedge \overline{CycB}\right)$	Transcription factors
$CycE$	$\left(E2F\wedge \overline{rB}\right)$	Cyclin E
$CycA$	$\left(E2F\wedge \overline{rB}\wedge \overline{Cdc20}\wedge \overline{Cdh1}\wedge \overline{Ubc}\right)\vee \left(CycA\wedge \overline{RB}\wedge \overline{Cdc20}\wedge \overline{Cdh1}\wedge \overline{Ubc}\right)$	Cyclin A
$p27$	$\left(\overline{CycD}\wedge \overline{CycE}\wedge \overline{CycA}\wedge \overline{CycB}\right)\vee \left(p27\wedge \overline{CycE}\wedge \overline{CycA}\wedge \overline{CycB}\wedge \overline{CycD}\right)$	$p27$ enzyme inhibitor
$Cdc20$	$CycB$	Activators of the Anaphase Promoting Complex
$Cdh1$	$\left(\overline{CycA}\wedge \overline{CycB}\right)\vee \left(Cdc20\right)\vee \left(p27\wedge \overline{CycB}\right)$
$Ubc$	$\left(\overline{Cdh1}\right)\vee \left(Cdh1\wedge Ubc\wedge \left(Cdc20\vee CycA\vee CycB\right)\right)$	E2 ubiquitin conjungating enzyme
$CycB$	$\left(\overline{Cdc20}\wedge \overline{Cdh1}\right)$	Cyclin B

Equations

$$I_H(n,N)=nlo{g}_{k}N$$

$$I(p)=log\frac{1}{p}$$

$$H(X)=-\sum _{x\in \mathcal{X}}p(x)logp(x)$$

$$H(\mathbf{X})=-\sum _{\mathbf{x}\in \mathcal{X}}p(\mathbf{x})logp(\mathbf{x})$$

$$H(Y|X)=-\sum _{x\in \mathcal{X}}\sum _{y\in \mathcal{Y}}p(x,y)logp(y|x)$$

$$H(X,Y)=H(X)+H(Y|X)=H(Y)+H(X|Y)$$

$$H(\mathbf{X})=\sum _{i=1}^{n}H(X_i|X_{i-1},X_{i-2},\dots ,X_1)$$

$$\sum _{x\in \mathcal{X}}p(x)ln\frac{p(x)}{q(x)}\ge 0$$

$$0\le H(X)\le log|\mathcal{X}|$$

$$H(X)=-\sum _{x\in \mathcal{X}}p(x)logp(x)\le -\sum _{x\in \mathcal{X}}p(x)logq(x)$$

$$-\sum _{x\in \mathcal{X}}p(x)logp(x)\le \sum _{x\in \mathcal{X}}p(x)log|\mathcal{X}|=log|\mathcal{X}|\sum _{x\in \mathcal{X}}p(x)$$

$$\sum _{x\in \mathcal{X}}p(x)=1$$

$$H(X)\le log|\mathcal{X}|$$

$$H(X)=\sum _{x\in \mathcal{X}}\frac{1}{|\mathcal{X}|}log|\mathcal{X}|=log|\mathcal{X}|\sum _{x\in \mathcal{X}}\frac{1}{|\mathcal{X}|}=log|\mathcal{X}|$$

$$I(X,Y)=\sum _{x\in \mathcal{X}}\sum _{y\in \mathcal{Y}}p(x,y)log\frac{p(x,y)}{p(x)p(y)}$$

$$f(\lambda x_1+(1-\lambda )x_2)\le \lambda f(x_1)+(1-\lambda )f(x_2)$$

$$f(x)\ge f(x_0)+f^{\prime }(x_0)(x-x_0)x,x_0\in [a,b]$$

$$\Delta \left\{\mathbf{x}\in {\mathbb{R}}^n:\sum _{i=1}^n{x}_i=1\text{ and }{x}_i\ge 0\forall i=1,\dots ,n\right\}$$

$$f\left(\sum _{i=1}^k{\lambda }_i{x}_i\right)\le \sum _{i=1}^k{\lambda }_{i}f(x_i)$$

$$I(X,Y)\ge 0$$

$$I(X,Y)=\sum _{x\in \mathcal{X}}\sum _{y\in \mathcal{Y}}p(x,y)log\frac{p(x,y)}{p(x)p(y)}$$

$$-I(X,Y)\le log\sum _{x\in \mathcal{X}}\sum _{y\in \mathcal{Y}}p(x,y)\frac{p(x)p(y)}{p(x,y)}$$

$$-I(X,Y)\le log\sum _{x\in \mathcal{X}}\sum _{y\in \mathcal{Y}}p(x)p(y)$$

$$\sum _{x\in \mathcal{X}}\sum _{y\in \mathcal{Y}}p(x)p(y)=1$$

$$-I(X,Y)\le log1=0⟺I(X,Y)\ge 0$$

$$I(X,Y)=H(X)-H(X|Y)$$

$$I(X,Y)=H(X)-H(X|Y)=H(Y)-H(Y|X)=I(Y,X)$$

$$I(X,Y)=H(X)+H(Y)-H(X,Y)$$

$$H(X|Y)\le H(X)$$

$$I(X,Y)\le \min (H(X),H(Y))$$

Full text

Ca’ Foscari University of Venice

Department of Environmental

Sciences, Informatics and Statistics

MSc Thesis in Computer Science

**Functional Dynamical Structures

in Complex Systems:

an

Information-Theoretic Approach

**

Marco Fiorucci

Supervisor:

Prof. Irene Poli

Academic Year 2014/2015

Dedicated to Corinne

Abstract

Understanding the dynamical behavior of complex systems is of exceptional relevance in everyday life, from biology to economy. In order to describe the dynamical organization of complex systems, existing methods require the knowledge of the network topology. By contrast, in this thesis we develop a new method based on Information Theory which does not require any topological knowledge. We introduce the Dynamical Cluster Index to detect those groups of system elements which have strong mutual interactions, named as Relevant Subsets. Among them, we identify those which exchange most information with the rest of the system, thus being the most influential for its dynamics. In order to detect such Functional Dynamical Structures, we introduce another information theoretic measure, called D-index. The experimental results make us confident that our method can be effectively used to study both artificial and natural complex systems.

Acknowledgements

This thesis is the result of the research I performed at the European Centre for Living Technology (ECLT) in Venice under the supervision of Prof. Irene Poli.

First of all, I wish to express my sincerest recognition to my supervisor, Prof. Irene Poli, for her continuous guidance and for introducing me in an international highly reputed research institute.

During this thesis, I had the opportunity to collaborate with Prof. Roberto Serra and Dr. Marco Villani from the university of Modena and with Dr. Andrea Roli from the university of Bologna. Working with them was particularly fruitful and inspiring.

I enjoyed a lot the pleasant environment in the ECLT with Dr. Matteo Borrotti, Dr. Alessandro Filisetti, Dr. Davide De March, Dr. Debora Slanzi. I would like to thank in particular my tutor, Dr. Alessandro Filisetti, who helped me break through the difficulties in my research and friendly supported me all over the thesis period.

I would like to extend my thanks to Prof. Marcello Pelillo and Prof. Nicoletta Cocco for their valuable discussions and productive collaboration.

More personally, I want to express my gratitude to my family. Since my childhood my parents have always trusted me and stimulated to learn. I’m also grateful to my brothers for all the experiences we have shared together.

Finally, a lovely thanks goes to Corinne for her unconditional support and for being a constant source of encouragement during my challenges.

Contents

1. 1 Introduction

2. 1.1 Functional Dynamical Structures

3. 1.2 Contributions

4. 1.3 Outline of the Thesis

5. 2 Information Theory and Functional Structures

6. 2.1 Information and Entropy

7. 2.2 A Few Properties of Entropy

8. 2.3 Mutual Information

9. 2.4 Integration

10. 2.5 Entropy Rate and Transfer Entropy

11. 3 Dynamical Cluster Index Methodology

12. 3.1 The Dynamical Cluster Index

13. 3.2 DCI: remarks and observations

14. 3.3 Integration vs. mutual information

15. 3.3.1 Experiments on a tuneable model

16. 3.4 The Relevant Subset Detection Algorithm

17. 3.5 Temporal Correlation: the D Index

18. 4 Experimental Results

19. 4.1 Boolean Networks

20. 4.2 Gene regulatory networks

21. 4.2.1 Mammalian cell cycle networks

22. 4.2.2 Matabolic pathway MAPK

23. 5 Conclusion

24. 5.1 Contributions and Novelty of the Work

25. 5.2 Future Work

List of Figures

1. 2.1 Standard simplex for $n=2$ (a) and $n=3$ (b).
2. 2.2 Integration scaling with respect to the subset size. The analyzed system is a catalytic reaction networks [40] composed of $26$ elements.
3. 2.3 Mutual Information scaling with respect to the subset size. The analyzed system is a catalytic reaction networks [40] composed of $26$ elements.
4. 3.1 Typical statistical behaviour of integration for different subsystem sizes. The system analyzed is a catalytic reaction system [39] composed of 26 molecular species.
5. 3.2 Typical statistical behavior of mutual information for different subsystem sizes. The system analysed is a catalytic reaction system [md-ecal2013] composed of 26 molecular species.
6. 3.3 Data generated according to a simple leader-followers model.
7. 3.4 Integration of $G_{1}$ as a function of $p_{noise}$ and $p_{copy}$.
8. 3.5 Mutual information between $G_{1}$ and $G_{2}$ as a function of $p_{noise}$ and $p_{copy}$.
9. 4.1 (a) The dynamical organization of the systems described in the text; on the edges, the Transfer Entropy, and on the nodes the $D$ index. Group C in case 5 has two values for the $D$ index, each value depicting the group C role in the relation with group A and group B, respectively. (b) The CRS identified by the DCI only in some systems described in the text. Each row denotes one CRS, composed of nodes whose entries are marked in black; on the right, the value of $t_{c}$. Note that in case 2 the third row is the correct identification of group B, and that in case 3 the third and fourth rows are group A and group B. For case 4 the output of the combination of the DCI analysis and sieving algorithm are presented; note that besides the correct group formed by nodes N03–N06 other CRS have high $t_{c}$ values, highlighting the presence of a larger but less evident dynamical group.
10. 4.2 (a) The masks identify by the combination of the DCI and the sieving algorithm for (a) the Mitogen Activated Protein Kinase (MAPK) cascade and (b) the two dynamical conditions of the mammalian cell cycle network. In each situation only the masks having significantly high $tC$ values are represented.
11. 4.3 (A) Basic model: the scheme of the three layers MAPK cascade reaction pathway is represented ( “” stands for the phosphorylation). The signal catalyzes the phosphorylation of M3K to M3K that in turn catalyzes the phosphorylation of M2K to M2K* and the successive phosphorylation of M2K* to M2K**. Finally M2K** performs the double phosphorylation of MK in MK** that is the final output of the MAPK cascade. P1, P2 and P3 dephosphorylate M3K, M2K and MK kinases respectively. V1-V10 stand for the involved reactions. Dashed lines with circle head represent catalysis; the figure highlights the presence of the three  “layers” described on the text. (B) Two distinct positive and negative feedbacks are added to the basic model: the negative feedback goes from MK** to the second layer (M2K, M2K* and M2K**) while the positive feedback goes from MK** to the first layer (M3, M3K*).
12. 4.4 The dynamical organization of MAPK system (basic model and S2 version). In italic the Transfer Entropy, in bold the $D$ index associated to the interested group within the relation. M1 group involves the first layer of 4.3a; M2 group involves the second layer, whereas groups M2_3 and M1_3 involve respectively the layers 2 and 3, and 1 and 3.

List of Tables

1. 4.1 The update rules of the boolean networks discussed on the text. Random($0.5$) denotes a Bernoulli distribution with probability $0.5$.
2. 4.2 Adapted from [6], the boolean regulatory network of mammalian cell cycle network and a short description of each node of the system.

1   Introduction

The whole is more than the sum of its parts.

Aristotele

In recent years complex systems attracted increasing interest among scientists of different fields like physics, computer science, chemistry, biology, ecology, economy and social science, just to name a few. This interest arises from the general idea that several natural systems (such as gene regulatory networks, social systems, etc.) and artificial systems (such as economical system, power grid, etc.) are composed of a large number of simple elements which show a collective emergence behavior that cannot be predicted from the description of each particular element.

There is no universal definition of complex system, but most researchers in the field would probably agree that it involves numerous components (agents), which may be simple both in terms of their internal characteristics and in the way they interact. These agents exhibit emergence activities, i.e. not derivable from the summations of the activity of the individual components, and a tangled hierarchical self-organization which cannot be understood by referring to simple hierarchical models.

Understanding the behavior of complex systems is of exceptional relevance, from both practical and theoretical reasons. This problem has been addressed by researchers which follow two basic approaches. The first one involves the development of mathematical models which abstract the most relevant features of the system to gain scientific insight into its dynamic. The concepts used in the development of such models belong to dynamical system theory, information theory, game theory, network theory, numerical methods and cellular automata. The second approach, instead, relies on the development of a more detailed and realistic models in order to simulate the emergent behavior which appears when the system is observed for a long time and length scales. The tools used in this approach are computational methods, like agent-based simulation and Monte Carlo simulation.

In this thesis we study the dynamical organization of complex systems using methods based on information theory. To do so, we model a system using a random field evolving in time. In particular, the dynamical behavior of each system component is described by means of a discrete stochastic process. We observe the system evolution in a given time interval and we estimate the probability distributions of the variables belonging to the random field. Then, exploiting measures based on information theory, which are known as information-theoretic measures, we look for relevant subsets (RSs). A RS is a group of variables which have a strong mutual interaction and that weakly interact with variables which are not in the relevant subset.

In order to find a measure for the interaction among variables, we exploit a fundamental property of the dynamical system from a complex system perspective which was introduced by Kauffman [14]. He stated that there is a continuous exchange of information among the system constituents and also between the environment and the system.

Historically, Tononi and Edelman [35] introduced an information-theoretic measure, the Cluster Index, to study biological neural networks which are close to a stationary state. More recently, such Cluster Index was profitably applied to the study of dynamical systems by Villani et al. [39] and was referred to as Dynamical Cluster Index (DCI).

Starting from such ideas, we study the two measures that contribute to the DCI, which are the Integration and the Mutual Information. We show that the analysis of the different parts of the index is extremely useful to better characterize the nature of the identified relevant subsets. In many cases the detected RSs have an intricate nested structure, so that it might not be clear which groups of variables are really important. To overcome this problem, we introduce a sieving procedure to extract only disjoint or partial overlapping RSs [8].

Furthermore, we extend the DCI introducing causality [43] in order to investigate the influence of found relevant subsets on the system dynamics. To achieve this goal we use information theoretic concepts like Entropy Rate and Transfer Entropy [28],[18]. In particular, on the basis of these measures, we introduce a directionality index [8] to find the direction of information flows which can take place between a RS and other parts of the system. The direction and the amount of information flows are needed to select the leaders, i.e. the most influenced RSs. We call these subsets functional dynamical structures (FDS) to stress their functional role in the system dynamics.

Several different application domains are investigated to test the effectiveness and the robustness of this measures. These domains include: random boolean networks, mammalian cell cycle networks and the Mitogen Activated Protein Kinase (MAPK) cascade signaling pathway in eukaryotes.

Finally, it’s worth noting that our method does not require any previous knowledge of underlying network topology of the complex system under study, but relies only on the values assumed by the random fields during the observation time interval.

1.1 Functional Dynamical Structures

The Dynamical Cluster Index is an extension of a measure introduced by Tononi and Edelman [35] to detect clusters in biological neural networks. In particular, Tononi et al. start from the idea that there are functional bounded regions in the brain despite the fact that all brain regions are widespread connected. These brain regions, called functional clusters, are composed of neurons which interact more strongly with neurons belonging to the same region than with neurons belonging to another brain part. The functional aspect of these clusters is suggested by the fact that there is a very fast signal transmission among intra-cluster neurons and they can thus influence many aspects of brain behavior. It’s worth noting that this method was applied to networks considering only fluctuations around stable asymptotic state. The time system evolution was indeed not studied. The first generalization of Tononi’s method was made by Villani et al. [39] in order to detect intermediate-level emergent structures in complex systems. They extended the cluster index to the study of truly dynamical systems and the DCI was introduced.

In this thesis, we further develop the method based on the Dynamical Cluster Index and we also show new interesting applications that highlight its effectiveness and help uncovering some of its features [7, 40]. There are of course several methods to identify clusters of nodes in a network based on its topology, but the DCI can be applied also without any prior knowledge on the network topology. This property is shared also by nonlinear correlation methods [13]; with respect to the latter, the advantage of the DCI is that it is not limited to binary relations but it can be applied to clusters of any size. Finally, we study the influence that detected relevant subsets have on the system by focusing on the causal interactions among variables through the different dynamical states. For this purpose, we introduce a directionality index which, as previously mentioned, is an information theoretic measure based on the transfer entropy. To the best of our knowledge, such approach in unprecedented in the study of the dynamical behavior of complex systems.

1.2 Contributions

The contribution of this thesis may be summarized as follows:

• •

we studied the two measures that contribute to the DCI, which are the Integration and the Mutual Information. We showed that the analysis of the different parts of the index is extremely useful to better characterize the nature of the identified relevant subsets. This led to a first publication [7];

• •

we extended the Dynamical Cluster Index methodology introducing a sieving algorithm. The results of this extension will appear in [40][8];

• •

we further extended the DCI methodology introducing causality in order to investigate the influence of found relevant subsets on the system dynamics. The results will be presented at European Conference of Artificial Life [8].

The results obtained in this thesis will be the subject of my oral presentation at the Student Conference on Complexity Science (SCCS2015) Conference in Granada in September 2015.

1.3 Outline of the Thesis

Chapter 2 reviews the fundamental information theoretic notions like entropy, joint entropy, conditional entropy, relative entropy, mutual information, integration and transfer entropy, which provide a basis for the development of our indexes. Chapter 3 introduces our contribution to the description of complex system behavior, namely the Dynamical Cluster Index, the D-Index and the methodology to detect functional dynamical structures [7, 40, 8]. Chapter 4 shows the results of the method application to various models. In particular, they include what follows: the random boolean network framework (in order to verify that the method is able to identify subsets that make sense) and leaders-followers model (to test the method robustness with respect to increasing noise levels). Furthermore, two applications to real systems are presented, namely mammalian cell cycle networks and the Mitogen Activated Protein Kinase (MAPK) cascade signaling pathway in eukaryotes. In these two latter applications, a peculiar discrete coding which describes the signs of their first-order time derivatives is used. In chapter 5, the novelty of the work and the future perspectives are presented.

2   Information Theory and Functional Structures

This word ”information” in communications theory relates not so much to what you do say, as to what you could say.

Warren Weaver

A complex system is composed of a huge amount of interacting parts. The dynamical behavior of the whole system cannot be understood if the analysis just focuses on the properties of the single parts. As a consequence, statistical methods play a key role to understand how the whole system works. In this framework, information theory provides a collection of measures which can be useful to determine the degree of interaction among interconnected elements. In this chapter, we provide a brief background on information-theoretic measures. In particular, we introduce the integration and the transfer entropy which are exploited by our method. In particular, the first one is used to measure the degree of interaction among the system elements while the latter is used to find the direction of the information flow which takes place between two subsets of the system.

2.1 Information and Entropy

Information theoretic measures give us a formal definition of the dynamical structures notion. Even if the concept of information may be referred to different meanings, we consider the formal presentation proposed by Claude Shannon [31], who developed a theory based on a statistic description of a communication system. Shannon generalized the information measure which was presented by Ralph Hartley [12], who stated that the information content of a message of length $n$ composed of symbols chosen from an alphabet of cardinality $N$ is

[TABLE]

where $k$ is an arbitrary base. Shannon generalized this concept associating a probability distribution to the symbols source, i.e. a symbol from the alphabet is chosen with a certain probability $p$. He stated that the information associated to the observation of a symbol with occurrence probability $p$ is

[TABLE]

Hence, the information that is received from a symbol observation depends on its occurrence probability, i.e less likely is the symbol appearance and more is the information that is gained from its observation. As a consequence, a source is modeled by a discrete random variable $X$ which can assume value on a non empty finite set $\mathcal{X}$ according to a probability distribution $p(x)=Pr\{X=x\},\,x\in\mathcal{X}$. If the source follows a uniform distribution, the Shannon information coincides with the Hartley information. In order to measure the average amount of transmitted information, Shannon entropy [31] was introduced in 1949.

Definition 2.1**.**

The entropy $H(X)$ of a discrete random variable $X$ with alphabet $\mathcal{X}$ and a probability mass function $p(x)=Pr\{X=x\},\,x\in\mathcal{X}$ is defined as follows

[TABLE]

As an alternative notation, we use $H(p)$. We consider the log in base two (bits as units of entropy) and we define $0\,log0=0$(since $x\,logx\rightarrow 0$ as $x\rightarrow 0$).

We naturally extend the definition of entropy to a random vector, i.e. a vector composed of two or more random variables.

Definition 2.2**.**

The entropy $H(\bm{X})$ of a discrete random vector

$\bm{X}=(X_{1},X_{2},\ldots,X_{n})$ with alphabet $\bm{\mathcal{X}}=\mathcal{X}_{1}\,\times\,\mathcal{X}_{2}\times\,\cdots\times\,\mathcal{X}_{n}$ and a joint probability distribution $p(\boldsymbol{x})=p(x_{1},x_{2},\ldots,x_{n})$ is defined as follows

[TABLE]

By analogy with probability theory, we define the conditional entropy using the notion on conditional probability which is a measure of the probability of an event given the knowledge of the occurrence of another event.

Definition 2.3**.**

The conditional entropy of two discrete random variables $X$ and $Y$ with alphabets $\mathcal{X}$ and $\mathcal{Y}$ and a joint probability distribution $p(x,y)$ is defined as follows

[TABLE]

The following theorem describes the relationship between joint entropy and conditional entropy.

Theorem 2.1** (Chain Rule).**

*Given two discrete random variables $X$ and $Y$ with alphabets $\mathcal{X}$ and $\mathcal{Y}$ and a joint probability distribution $p(x,y)$, then

[TABLE]

We extend the previous result to a generic random vector.

Theorem 2.2**.**

*Given a discrete random vector $\mathbf{X}=(X_{1},X_{2},\ldots,X_{n})$ with alphabet $\bm{\mathcal{X}}=\mathcal{X}_{1}\,\times\,\mathcal{X}_{2}\times\,\cdots\times\,\mathcal{X}_{n}$ and a joint probability distribution $p(\boldsymbol{x})=p(x_{1},x_{2},\ldots,x_{n})$, then

[TABLE]

2.2 A Few Properties of Entropy

We present few properties of entropy which are exploited to construct our method. In order to prove these properties, we introduce one of the most important inequality in information theory.

Lemma 2.3** (Gibbs’ Inequality).**

*Given a random variables $X$ with alphabet $\mathcal{X}$ and two probability mass functions $p(x)\,,q(x)\,,x\in\mathcal{X}$, then

[TABLE]

with equality iff $p(x)=q(x)$, for all $x\in\mathcal{X}$.

Gibbs’ inequality allows us to introduce an upper bound for Shannon entropy.

Theorem 2.4**.**

*Given a discrete random variable $X$ with alphabet $\mathcal{X}$, which contains a number of symbols equals to $|\mathcal{X}|$, and a probability mass function $p(x),\,x\in\mathcal{X}$, then

[TABLE]

*with equality iff $p(x)$ is a uniform distribution. *

Proof.

$0\leq p(x)\leq 1$ implies that $log(\frac{1}{p(x)}\geq 0)$. It follows that $H(X)\geq 0$.

The equality holds iff $\exists x\in\mathcal{X}\,:\,p(x)=1$ because $\forall x\in\mathcal{X}$

$p(x)\,logp(x)=0\,\longleftrightarrow\,p(x)=0\,\vee,\ p(x)=1\,\longleftrightarrow,\ \exists x\in\mathcal{X}\,:\,p(x)=1.$

From Gibbs’ inequality we obtain

[TABLE]

Let $q(x)$ be uniformly distributed. Hence

[TABLE]

Since

[TABLE]

We obtain

[TABLE]

Finally we prove that if $X$ is uniformly distributed then $H(X)=log|\mathcal{X}|$.

[TABLE]

∎

For our purpose, we interpret the entropy as a measure of a lack of knowledge about the system. In particular, we compute the probability of each possible state of a dynamical system after observing its evolution through a certain time interval. The gained information depends on the observed trajectory. Therefore, we can see the entropy as the average information which is gained by the system observation. This approach is borrowed from statistical mechanics, where a macrostate of a system is described associating a probability distribution over the possible microstates which are unknown. In this contest, the entropy is seen as a measure of the disorder of the systems: a connection between the information theory and the thermodynamics can be thus grasped. [20].

2.3 Mutual Information

In the study of complex systems, we have to deal with nonlinear systems which do not satisfy the superposition principle. We are thus interested in extensive measures which are useful to characterize a complex systems at least in some particular cases. Intuitively, the information emitted from two independent sources should be equal to the sum of the two individual sources. The log-dependence in the definition of the Shannon entropy makes the entropy of two independent random variables an additive quantity. Hence, the system entropy is maximum when every pair of system elements is independent. Otherwise, the correlation between two or more elements of the systems causes an entropy reduction.

We introduce mutual information which can be used to measure the degree of dependence of two random variables.

Definition 2.4**.**

The mutual information $I(X,Y)$ of two discrete random variables $X$ and $Y$ with alphabets $\mathcal{X}$ and $\mathcal{Y}$ and a joint probability distribution $p(x,y)$ is defined as follows

[TABLE]

Hence, the mutual information is symmetric under the exchange of its arguments.

We present some relationships between entropy, conditional entropy and mutual information. In order to prove these properties, we introduce another extremely important information-theoretic inequality: the Jensen’s inequality. We recall some mathematical notions which are used in the inequality definition.

Definition 2.5** (Convex function).**

Given a function $f:\,[a,b]\rightarrow\mathbb{R}$. $f$ is said to be convex if for all $x_{1},x_{2}\in[a,b],\,\lambda\in[0,1]$

[TABLE]

Theorem 2.5**.**

*Given $f:\,[a,b]\rightarrow\mathbb{R}$. Let be $f\in\mathcal{C}^{1}$. $f$ is convex iff

[TABLE]

Definition 2.6** (Standard Simplex).**

[TABLE]

Theorem 2.6** (Jensen’s Inequality).**

*Given $f:\,[a,b]\rightarrow\mathbb{R}$, $x_{1},\ldots,x_{k}\in[a,b]$, $\boldsymbol{\lambda}=(\lambda_{1},\ldots,\lambda_{k})\in\Delta_{k}$. Let be $f$ convex, then

[TABLE]

whith equality iff $f$ is strictly convex.

Jensen’s inequality allows us to introduce a lower bound for mutual information.

Theorem 2.7**.**

*Let $X$ and $Y$ be two discrete random variables with alphabets $\mathcal{X}$ and $\mathcal{Y}$ and a joint probability distribution $p(x,y)$, then

[TABLE]

Proof.

By definition

[TABLE]

Since f(z) = log(z) is a concave function, we can apply Jensen’s inequality changing the inequality direction

[TABLE]

[TABLE]

From the normalization axiom of probability theory, it follows

[TABLE]

Hence

[TABLE]

∎

There are the following relationships between entropy and mutual information which can be proved using (2.19), (2.6).

[TABLE]

From (2.6) it follows that

[TABLE]

[TABLE]

The following theorem states that, on average, the knowledge about a random variable can reduce the uncertainty in another random variable.

Theorem 2.8**.**

*Let $X$ and $Y$ be two discrete random variables with alphabets $\mathcal{X}$ and $\mathcal{Y}$ and a joint probability distribution $p(x,y)$, then

[TABLE]

with equality iff X\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y

It’s worth nothing that $H(X|Y=y)$ may be greater than $H(X)$. The equation 2.29 is true only for average quantity, i.e $H(X|Y)$ [5].

The following theorem introduces an upper bound for mutual information.

Theorem 2.9**.**

*Let $X$ and $Y$ be two discrete random variables with alphabets $\mathcal{X}$ and $\mathcal{Y}$ and a joint probability distribution $p(x,y)$, then

[TABLE]

Proof.

Entropy is non negative quantity so, in order to find the maximum of the mutual information, we can find the maximum of the sum H(X)+H(Y) and the minimum of the joint entropy H(X,Y).

From (2.6) we obtain

[TABLE]

Hence, for every joint distribution p(x,y) holds

[TABLE]

Hence

[TABLE]

Finally we obtain

[TABLE]

∎

Hence, from the equations 2.20 and 2.30, we can state that I(X,Y) takes values over the interval $[0,\min(H(X),H(Y))]$.

2.4 Integration

As we previously hinted at, mutual information can be used to measure the degree of statistical dependence between two elements or subsets of elements when they are associated to two random variables or two random vectors respectively.

We can generalize this concept in order to measure the degree of interaction among $k$ elements of the system. To do so, let us consider a system $U$ composed of $N$ elements. The system dynamics is described using a random field, i.e. $\boldsymbol{X}=(X_{1},X_{2},\ldots,X_{N})$, where $X_{i}$ is a discrete random variable associated with the $i$-th element. $X_{i}$ assumes value in a finite alphabet $\mathcal{X}_{i}$. Let $S$ be a set of $k$ elements of the system, such that $k\leq N$, i.e. $S$ is a subset of $U$. Without loss of generality, let $X_{1},X_{2},\ldots,X_{k}$ be the random variables associated with the $k$ elements.

If all the random variables are mutually independent, the entropy of the whole subset $S$ can be computed as follows

[TABLE]

If there are interactions among the elements of $S$, then a degree of statistical dependence arises among the $k$ random variables. It follows an entropy reduction that can be measured by means of integration which is an information-theoretic measure defined as follows [37]

Definition 2.7**.**

[TABLE]

where $\mathbf{S}$ is the random vector composed of all the variables which are in $S$.

The following theorem introduces a lower and an upper bounds for integration.

Theorem 2.10**.**

*Given a discrete random vector $\mathbf{X}=(X_{1},X_{2},\ldots,X_{k})$ with alphabet $\bm{\mathcal{X}}=\mathcal{X}_{1}\,\times\,\mathcal{X}_{2}\times\,\cdots\times\,\mathcal{X}_{k}$ and a joint probability distribution $p(\boldsymbol{x})=p(x_{1},x_{2},\ldots,x_{k})$. Let be $\mathcal{X}_{1}=\mathcal{X}_{2}=\ldots=\mathcal{X}_{k}=\mathcal{X}$, then

[TABLE]

Proof.

From 2.7 we obtain

[TABLE]

From 2.29 we can state that for $i=1,\ldots,k$

[TABLE]

Hence

[TABLE]

To prove that $I(\mathbf{X})\leq(k-1)\,log|\mathcal{X}|$ we use 2.9

[TABLE]

If $X_{i}$ is uniformly distributed for $i=1,\ldots,k$ then $H(X_{i})$ assumes its maximum value (i.e. $log|\mathcal{X}_{i}|)$. In this case, the minimum value of $H(\mathbf{X})$ is achieved when $\mathbf{X}$ assumes, with probability greater than zero, only $|\mathcal{X}|$ out of the $|\mathcal{X}|^{k}$ possible values. Under this assumption it follows that $H(\mathbf{X})\leq log|\mathcal{X}|$ and hence

[TABLE]

[TABLE]

∎

For all of these features we can use integration and mutual information to detect relevant subsets (RS) in complex dynamical systems. As we previously hinted at, a RS is a subset of state variables which heavily influence the system dynamics. It is composed of variables which have a strong mutual interaction and that weakly interact with variables which do not belong to the subset. Integration can be used as internal criterion: it is a measure of the internal interaction among the elements of a RS. At the same time, mutual information can be used as an external criterion: it is a measure of the interaction between the RS and the rest of the system (or between two relevant subsets). A candidate relevant subset is characterized by a high integration value and a low mutual information value.

It’s worth noting that both integration and mutual information scale with the size of a relevant subset. We can see this behavior in 2.2, 2.3, where integration and mutual information values are computed analyzing the dynamics of a catalytic reaction network composed of 26 molecules (see chapter $5$ for more details). In order to compare RSs of different sizes, a normalization procedure has been developed. It will be presented in the next chapter.

2.5 Entropy Rate and Transfer Entropy

As previously explained, mutual information can be used to measure the interaction between two relevant subsets. Since mutual information is a symmetric measure, it does not provide indications on the information flows directions. Furthermore, it does not take into account dynamical information which are needed to detect the subsets which most influence the system dynamics. We say that a subset $A$ affects the behavior of a subset $B$ if the net information flow goes from $A$ to $B$. Hence, the direction and the amount of information flows is needed to select the leaders, i.e. the most influential RSs. We call these subsets Functional Dynamical Structures (FDS) to stress their functional role in the system dynamics. A FDS sends a high amount of information to the other subsets conditioning their behavior. In order to measure the dynamical information flow between two relevant subsets, we used a new measure based on the transfer entropy [28]. We recall some notions from stochastic processes theory [2] and we define the entropy rate which is used in the transfer entropy definition.

Definition 2.8** (Stochastic Process).**

A stochastic process $\{X_{i}\}$ is a sequence of discrete random variables with alphabet $\bm{\mathcal{X}}=\mathcal{X}^{k}$ and a joint probability distribution $p(\boldsymbol{x})=p(x_{1},x_{2},\ldots,x_{n})$ for $n=1,2,\ldots$.

Definition 2.9** (Markov Process).**

A stochastic process $\{X_{i}\}$ is called a Markov Process or Markov Chain if

[TABLE]

for $n=1,2,\ldots$ and for all $x_{1},\ldots,x_{n+1}$.

Definition 2.10** (Stationary Markov Process).**

A Markov Process $\{X_{i}\}$ is stationary if

[TABLE]

i.e. $p(x_{n+1}|x_{n})$ is independent with respect to $n$.

In order to introduce the notion of entropy rate, let us consider a system $U$ composed of $N$ elements. The system dynamics is described using a random field evolving in time, i.e. $\boldsymbol{X}(t)=(X_{1,t},X_{2,t},\ldots,X_{N,t})$, where $X_{i,t}$ is a discrete random variable associated with the $i$-th element at time $t$. $X_{i}$ assumes value in a finite alphabet $\mathcal{X}_{i}$. Hence, if we fix $t$, $\boldsymbol{X}$ is a random field while, if we fix $i$, $X_{i}(t)$ is a stochastic process. Let assume that $X_{i}(t)$ is a stationary Markov process of order k. Under this assumption, the probability to find $X_{i}$ in state $x_{i}$ at time $t+1$ satisfies the following relation

[TABLE]

We adopt the following shorthand notation proposed by Schreiber [28]

[TABLE]

for words of length $k$.

Let us now define the entropy rate as follows

Definition 2.11** (entropy rate).**

Given a stationary Markov process $\{X_{i}\}$ of order $k$, the entropy rate $h_{X_{i}}$ is

[TABLE]

for $n=1,2,\ldots$ and for all $x_{1},\ldots,x_{n+1}$.

Since $p(x_{i,t+1}|x_{i,t}^{k})=p(x_{i,t+1}^{k+1}/p(x_{i,t}^{k})$, we can express $h_{X_{i}}$ as follows

[TABLE]

where $H_{X_{i}^{k+1}}$ and $H_{X_{i}^{k}}$ are the entropies of the process ${X_{i}}$ considering two delay vectors of dimension $k+1$ and $k$ respectively.

In order to find the direction of the information flow which takes place between two elements $X_{i}$ and $X_{j}$, we have to generalize the entropy rate. The idea is to measure the deviation from the following generalized Markov property

[TABLE]

If there is an information flow from $X_{j}$ to $X_{i}$, then the transition probabilities of $X_{i}$ depends on the states of $X_{j,t}$. In this case the entropy rate is

[TABLE]

Otherwise, if there is not an information flow from $X_{j}$ to $X_{i}$, then the transition probabilities of $X_{i}$ are independent from the state of $X_{j}$. In this case the entropy rate becomes

[TABLE]

In order to measure the degree of statistical dependence of the transitional probabilities of $X_{i}$ on that of $X_{j}$, we can use the transfer entropy which is defined as follows

[TABLE]

At the same time, we can define the transfer entropy from $X_{i}$ to $X_{j}$ as follows

[TABLE]

We can express transfer entropy in terms of conditional entropies as follows

[TABLE]

[TABLE]

It’s worth noting that the transfer entropy is explicitly non-symmetric: $T_{X_{j}\rightarrow X_{i}}$ is in general different from $T_{X_{i}\rightarrow X_{j}}$. In the next chapter, we present an index to detect the directionality of information flow between two relevant subsets.

3   Dynamical Cluster Index

Methodology

An idea which can be used once is a trick. If it can be used more than once it becomes a method.

G. Polya and S. Szegö

3.1 The Dynamical Cluster Index

The DCI is an extension of the Functional Cluster Index (CI) introduced by Edelman and Tononi in 1998 [35] and aimed at detecting functional clusters in brain regions. In our work [40], we relax the stationary constraint and extend the CI to actual dynamical systems, so as to apply it to a wide range of system classes, from abstract models to biological models.

As in our work we rely on observational data, probabilities are estimated by the relative frequencies of the values observed for each variable.

Let us now consider a system $U$ composed of $K$ variables (e.g. agents, chemicals, genes, artificial entities) and suppose that $S_{k}$ is a subset composed of $k$ elements, with $k<K$. The value $DCI(S_{k})$ is defined as the ratio between two aforementioned measures: the integration ($I$) and the Mutual Information ($MI$).

As previously stated, $I(S_{k})$ measures the statistical independence of the $k$ elements in $S_{k}$ and is defined as:

[TABLE]

As previously stated, mutual information $MI(S_{k};U\backslash S_{k})$ measures the mutual dependence between subset $S_{k}$ and the rest of the system $U\backslash S_{k}$ and it is defined as:

[TABLE]

hence, the DCI,

[TABLE]

The value of DCI is not defined if $MI(S_{k};U\backslash S_{k})=0$. However, a vanishing mutual information is a sign of separation of the investigated subset from the rest of the system, and therefore the subset must be studied separately. It is worth noting that the DCI scales with the subset size. In [40] we show a procedure to normalize it, nevertheless a better approach is that of assessing the statistical significance of the DCI of $S_{k}$ by means of a statistical significance index [35]:

[TABLE]

where $\left\langle DCI_{h}\right\rangle$ and $\sigma(DCI_{h})$ are respectively the average and the standard deviation of the $DCI$ of a sample of subsets of size $k$ extracted from a reference system $U_{h}$ randomly generated according to the frequency of each single state in $U$ and $\nu=\left\langle MI_{h}\right\rangle/\left\langle I_{h}\right\rangle$ is the normalization constant. It is worth noting that the aim of the reference system is that of quantify the finite size effects affecting the information theoretical measures on a random instance of a system with finite dimensions.

3.2 DCI: remarks and observations

As previously stated, the two factors used to compute the $DCI$ of a candidate RS are the integration of its elements and the mutual information between the subset and the rest of the system. Therefore, the $DCI$ of a subset $S$ of variables in a system—i.e., a possible candidate RS—is then estimated by collecting a set of system states and then by computing the entropy values of the possible combinations of variables in $S$. Therefore, all we need is just a collection of observations of the system and, in principle, no information on the topology or on the internal mechanics of the system are required. The rationale behind the ratio between $I$ and $MI$ is that a candidate RS should express higher interactions among its components than with the rest of the system. In spite of the simplicity of the concept, one must carefully consider the meaning of the two quantities and the way they are combined. In fact, the implicit assumption here is that it holds $MI>0$ and that the orders of magnitude of $I$ and $MI$ are reasonably close. However, there might be cases in which $MI\approx 0\;\wedge\;I>0$, which denotes integrated subsets dynamically independent from the rest of the system. Furthermore, we may be interested in finding the most integrated subsystems among the ones that exchange information with the rest of the system. Moreover, it is important to remark that insightful information can be also provided by analyzing the two factors composing the $DCI$, before computing their ratio. For this reason, an assessment of the information brought by the individual values of $I$ and $MI$ can shed light on the potential of the application of these information-theoretic measures to detect candidate RSs.

3.3 Integration vs. mutual information

To assess the contribution of $I$ and $MI$ we first study their statistical behavior with respect to subsystem size. We analyzed several dynamical systems and the results obtained show that the general trend respects the theoretical findings [35], as shown in Figure 3.1 and Figure 3.2.

As we can observe, $I$ scales approximately linearly, whilst $MI$ is non-monotonic. We have previously proved that the maximal value of $I$ for a subset of size $n$ equals $n-1$. Hence, we can define a rescaled version of $I$, $rI=I/n-1$, useful to compare the integration of subsets of different dimensions. Moreover, the addition of a variable $x$ to the considered subset $S$ increases $I$ by 1 if $x$ deterministically depends on any variable in $S$, while it leaves $I$ unchanged if it assumes random values. These properties may be useful to reckon the relative importance of integration values computed for subsystems of different size.

3.3.1 Experiments on a tuneable model

We are interested in understanding the individual informative contribution of $I$ and $MI$. To this aim, we study these measures on a simple model in which the integration among variables in a subsystem under observation and its mutual information can be tuned by acting on two parameters. The model abstracts from specific functional relations among elements of the system and could resemble a basic leader-followers model (see Figure 3.3). The system is composed of a vector of $n$ binary variables $x_{1},x_{2},\ldots,x_{n}$, e.g., representing the opinion in favor or against a given proposal. The model generates independent observations of the system state, i.e. each observation is a binary $n$-vector generated independently of the others, on the basis of the following rules:

• •

Variables are divided into two groups, $G_{1}=[x_{1},\ldots,x_{k}]$ and $G_{2}=[x_{k+1},\ldots,x_{n}]$;

• •

$x_{1}$ is called the leader and it is assigned a random value in $\{0,1\}$;

• •

the value of the followers $x_{2},\ldots,x_{k}$ is set as a copy of $x_{1}$ with probability $1-p_{noise}$ and randomly with probability $p_{noise}$;

• •

the values of elements of $G_{2}$ are assigned as a copy of a random element in $G_{1}$ with probability $p_{copy}$, or a random value with probability $1-p_{copy}$.

It is possible to tune the integration among elements in $G_{1}$ and the mutual information between $G_{1}$ and $G_{2}$ by changing $p_{noise}$ and $p_{copy}$. Note that, given significant level of integration, we have two notable cases:

• (a)

$MI\approx 0\rightarrow$ isolated (possibly integrated) RS;

• (b)

$MI\gg 0\rightarrow$ integrated and segregated cluster.

The possible scenarios which can be obtained by tuning $p_{noise}$ and $p_{copy}$ can be conveniently illustrated by a 3-dimensional plot. Figure 3.4 shows the behavior of integration of $G_{1}$ as a function of $p_{noise}$ and $p_{copy}$. We can observe that it is a decreasing function of $p_{noise}$, while it is independent of $p_{copy}$ (by definition, indeed). The behavior of the mutual information between $G_{1}$ and $G_{2}$ is depicted in Figure 3.5. As we can observe, $MI$ increases fast with $p_{copy}$, as this parameters increases the correlation between variables in $G_{2}$ and $G_{1}$. Moreover, it also increases with $p_{noise}$, but the reason is that the correlation among variables in $G_{1}$ increases the randomness of variables in $G_{1}$, which behaves similarly to the variables in $G_{2}$.

The case of $MI\approx 0$ corresponds to the situation in which $G_{1}$ is almost completely independent of $G_{2}$ and can be easily detected by observing only $MI$. Conversely, if we are interested in discovering $G_{1}$ as significant RS, then we would consider cases in which both $I$ and $MI$ are significantly high. In our experiments, these cases correspond to $p_{noise}<0.3$ and $0<p_{copy}<0.5$. In these scenarios, we find that it is possible to detect $G_{1}$ by using only $I$, divided by $n-1$. It is important to mention that, when $p_{noise}$ is slightly higher, it is necessary to resort to the computation of the DCI to detect group $G_{1}$.

3.4 The Relevant Subset Detection Algorithm

As we previously hinted at, the dynamical cluster index is used to detect the Candidate Relevant Subsets (CRSs) of a system. In order to develop a methodology to detect CRSs, let us consider a system $U$ composed of $N$ elements. The system dynamics is described using a random field, i.e. $\boldsymbol{X}=(X_{1},X_{2},\ldots,X_{N})$, where $X_{i}$ is a discrete random variable associated with the $i$-th element. $X_{i}$ assumes value in a finite alphabet $\mathcal{X}_{i}$. Let $S$ be a set of $k$ elements of the system, such that $k\leq N$, i.e. $S$ is a subset of $U$. In order to get a list of Candidate Relevant Subsets (CRSs), we compute the $DCI$ of every possible subset of variables in $U$ and ranking the subsets by $DCI$ values. The subsets occupying the first positions are most likely to play a relevant role in system dynamics. For large-size systems, exhaustive enumeration is computationally impractical as it requires to enumerate the power set of $U$. In this case, we resort to a genetic algorithms.

The list of CRSs can be ranked according to the significance of their $DCI$. We can directly use this ranking to identify by hand the relevant CRSs for the dynamics of the system. Nevertheless, in many cases this analysis might return a huge list of entangled CRS, so that a direct inspection is required for explaining their relevance. To this aim, we present a DCI analysis post-processing sieving algorithm to reduce the overall number of CRS to manually tackle. The main assumption of the algorithm is that if a CRS $A$ is a proper subset of CRS $B$, i.e. $A\subset B$, then only the subset with the higher DCI is maintained between the two. Thus, only disjoint or partially overlapping CRSs are retained: the used assumption implies that the remaining CRSs are not further decomposable, forming in such a way the “building blocks” of the dynamical organization of the system. The pseudo-code of the algorithm is presented in Algorithm 1.

3.5 Temporal Correlation: the D Index

Although by the application of the DCI, CRSs are detected, this measure does not provide indications neither on the quantity nor on the direction of the information which flow among subsets. To this aim we apply the directionality index proposed in [18].

Let $X$ and $Y$ be two random variables—or, equivalently, two sets of variables. As mentioned in the previous chapter, we can define the entropy rate of $X$ as the average number of bits needed to encode a successive state of $X$ if all the previous states are known, considering that the value of $X$ at $t+1$ depends either on $X$ and $Y$ at the time $t$, eq. 3.5a, or just on the value of $X$ at the time $t$, eq. 3.5b.111Note that the temporal dependency is not necessarily of unitary lag, i.e. $t-1\rightarrow t$. For a complete assessment of the statistical dependency of $X$ on $Y$ one should sum over $t-1,t-2,\ldots,t-T$, where $T$ is the observation time limit. Nevertheless, note that (i) in this work we are analyzing Markovian systems, whose behavior depends only from the immediately previous step and (ii) although TE is not a direct measure of causal effect, the use of short history length alters the character of the measure towards inferring causal effect [22].

[TABLE]

then, the transfer entropy $T$ is defined as the difference between the aforementioned entropy rates. $T$ describes how the uncertainty of $X$ is reduced by knowing the previous states of $Y$ and $X$ itself, eq. 3.6.

[TABLE]

Thus, $T_{Y\rightarrow X}$ can be described in term of entropy as:

[TABLE]

and, since the $T_{Y\rightarrow X}$ describes the information moving from $Y$ to $X$, and the transfer entropy is not symmetric, the information from $X$ to $Y$ is computed as well, eq. 3.8.

[TABLE]

Once that $T_{Y\rightarrow X}$ and $T_{X\rightarrow Y}$ are known, the directionality $D$ of the information flow between $X$ and $Y$ can be measured as:

[TABLE]

where $D_{X\rightarrow Y}=1$ indicates that all the information moves from $X$ to $Y$, i.e. absence of information flow from $Y$ to $X$ and, conversely, $D_{X\rightarrow Y}=-1$ indicates that $X$, with respect to $Y$, is just an information receiver and not an information sender.

It is worthwhile to notice that $D$ does not provide any indication about the amount of information exchanged between the variables, but it only provides suitable indications on the direction of the information flow.

4   Experimental Results

It doesn’t matter how beautiful your theory is, it doesn’t matter how smart you are. If it doesn’t agree with experiment, it’s wrong.

Richard P. Feynman

4.1 Boolean Networks

We consider a system composed of $12$ boolean nodes updated on the basis of either a boolean function or a random boolean value generator. Nodes update their state in parallel and synchronously. We illustrate the results of $5$ instances of this network, defined in Table 4.1 111Note that the size of these systems allows for an exhaustive enumeration of all the possible groups, allowing their complete assessment. It is worth remarking that each perturbation is introduced after the system has recovered a stable dynamical situation.. The $5$ instances share a common structure but differ in specific dynamical organizations of some nodes.

• •

In case 1, we consider two integrated groups of three nodes (namely, group A and group B), by assigning at each node the $XOR$ function of the other two nodes in the group. In this case the seaving algorithm filtered the $94\%$ of the evaluated CRS, making it possible to easily identify the subsets responsible for the dynamical organization of the system. Only the two correct CRSs have high $t_{c}$ values, whereas all the other ones have $t_{c}$ values lower of more than 2 orders of magnitude. Moreover, no information exchange takes place between group A and group B, as they are structurally independent.

• •

Case 2 derives from case $1$ by introducing in the first node of group B a further dependence from the last node of group A, hence introducing information transfer from group A to group B. The combination of the DCI analysis and the sieving algorithm correctly recalls the dynamical organization of the system—i.e. group A influences the behavior of (a part of) group B. We observe that the whole group B (Figure 4.1b) was anyway ranked high w.r.t. its DCI significance, but it was discarded by the sieving algorithm because the dynamics of its first node is influenced also by group A and so the assessment of the whole group B is weakened. In general, the amount of this difference depends on the strength of the forces that influence the interface nodes and the elements interfacing different CRS can be detected by a simple comparison between the DCI analysis and the sieving algorithm outputs.

• •

Case 3 derives from case 2 by introducing a further dependence of the first node of group A from the last node of group B: again, the combination of the DCI and the sieving algorithm detects the interface nodes and the underlining dynamical system organization. Note that the asymmetry in transfer entropies (and therefore in $D$ index) are due to differences in the initial conditions of the boolean network trajectories: a shift in initial conditions can change the direction of this asymmetry.

• •

In case 4 five heterogeneously linked nodes are influenced by a triplet that is identical to that of group A. The combination of the dynamical rules of the nodes and their initial condition makes the dynamical behavior of the sixth node always in phase with the triplet, so our analysis individuates this quartet as the most relevant CRS. The other dynamical relations are not sufficiently strong to coordinate all the $8$ nodes, nevertheless their influence creates some masks222CRSs can be represented by rows, where entries corresponding to variables belonging to CRS are marked in black ( “masks”in the following) having high $t_{c}$ values, partially overlapped with the leading quartet. The overlap of these masks indirectly indicates the presence of a greater group with respect to the winning quartet, but having a less evident dynamical presence (see Figure 4.1b).

• •

Case 5 derives from case 1 by adding two nodes whose dynamical behavior directly depends on nodes of both group A and group B: these 8 nodes form therefore a group clearly different from the remaining 4 random nodes, as they are interdependent and ruled by deterministic functions. This group is identified by the plain DCI method, but the combination of DCI and sieving algorithm strikingly enlightens the interpretation of two triplets directly influencing a couple of nodes (Figure 4.1a).

4.2 Gene regulatory networks

In this section we show the application of our method to two models of regulatory networks: a model of mammalian cell cycle network (MCC in the following), as “booleanized” in [6]—see Table 4.2 for the chosen boolean model—, and a model of one of the major cellular signal transduction pathways, known as the Mitogen Activated Protein Kinase (MAPK) cascade [42].

4.2.1 Mammalian cell cycle networks

In [6] the authors provides a boolean dynamical model of the mammalian cell cycle, able to reproduce the main characteristics of the succession of molecular events leading to the reproduction of the genome of a cell and its division into two daughter cells.

Mammalian cell division must be coordinated with the overall growth of the organism; this coordination is achieved through extra-cellular signals whose balance decides whether a cell will divide or remain in a resting state. The positive signals or growth factors ultimately elicit the activation of Cyclin D (CycD) in the cell. In the proposed model CycD thus represents the input and its activity is considered constant. By pointing the interested reader to [6] for the details, for now it is enough to say that in absence of CycD the system presents a unique stable attractor where only Rb, p27 and Cdh1 are active, whereas in its presence E2F, CycE, CycA, Cdc20, Cdh1, UbcH10 and CycB cycle with a period of length 7. We perturb both asymptotic states, obtaining in each case only one group (composed of E2F, CycE, CycA, Cdc20, Cdh1, UbcH10 and CycB in the first case, and of Rb, E2F, p27, Cdc20, UbcH10 and CycB in the second case). The leading groups of the two situations are different, but in each case the other CRS identified overlap with these groups and their sum cover the whole system, indicating the presence of a single coordinated pattern. So, the analysis indicates that the elements of the mammalian cell cycle network act as a single compact engine, see Figure4.2b.

4.2.2 Matabolic pathway MAPK

The MAPK pathway (evolutionarily conserved from yeasts to humans) responds to a wide range of external stimuli, triggering growth, cell division and proliferation [27]. [27] also introduce the models considered in our analysis. The basic model is composed of reactions that are quite well-established from an experimental viewpoint, and it has the hierarchical structure shown in Figure 4.3a. The three chemicals identified as the core of this three-layered system are the $MAPKKK$, $MAPKK$ and $MAPK$ kinases (respectively $M3K$ , $M2K$ and $MK$ for short) [42]. $M3K$ is activated by means of a single phosphorylation whereas $M2K$ and $MK$ are both activated by double phosphorylation. Parallel to the phosphorylation by kinases, phosphatases present in the cellular volume can dephosphorylate the phosphorylated kinases (Figure 4.3a shows the schema of the $MAPK$ cascade where each layer of the cascade is dephosphorylated by a specific phosphatase). Note that superimposed on the three-layered structure of substrates-product reactions there is the properly called $MAPK$ signalling cascade, a linear chain of catalysis (dashed lines in Figure 4.3a) that transmit the external signal from $M3K^{*}$ to $MK^{**}$.333The symbol “$*$” indicates the phophorylated version of the molecule.

When the external signal and the concentrations of the phosphatases are kept constant, a top-down reaction scheme as the one described in Figure 4.3a leads to fixed-point solutions. On the other hand, oscillations have been reported in the MAPK cascade [32] and, in order to account for them, [27] adopt a models with feedback, one of which is described in Figure 4.3b. This variant (called S2 in the following) is characterized by a configuration of the activating and inhibiting interactions among layers that alters the  “layered” structure of the basic model, which is no longer strictly hierarchical. This alternative model is grounded on experimental data; we will not enter here a discussion about the merits and limits of the model, referring the interested reader to the original paper, but we will take it “for granted” and we will apply our method to test whether it can discover significant dynamical organization, without any prior knowledge of the interactions, but on the sole basis of the dynamics of concentrations. We simulate the [27] models with the CellDesigner software [9, 10], keeping the P1, P2 and P3 phosphatases as constant (as they do) obtaining the same asymptotic states shown by the authors. In order to apply our method we perturb the asymptotic states of these models: in particular, we focus our analysis on kinases perturbations. In particular, we perform $10$ perturbation cycles, each cycle involving the perturbation of each single kinase and the successive relaxing to a stable situation before the subsequent perturbation [40]. The stable situations that are reached can show both oscillating (S2 system) or constant concentrations (basic system). Concentration changes are more significant than their absolute values (it is important to monitor the variables that are changing in coordinate way); therefore the continuous concentration values are represented according to a three levels code related to the sign of the time derivatives at time $t$ ( “decreasing concentration”,  “no significant change”,  “increasing concentration”). The combination of DCI and sieving algorithm applied to the basic MAPK model detects two dynamical groups: the first including the first layer of Figure 4.3a and the second including the other two layers). The two groups exchange information, the second transmitting more information to the first one (see Figure 4.4). The introduction of the feedbacks changes system dynamics: there are still two dynamically relevant groups, now composed of the second layer and by the other two layers, respectively. The analysis therefore suggests that the MAPK system may be decomposed in two dynamically distinct parts.

5   Conclusion

This is not the end. It is not even the beginning of the end. But it is, perhaps, the end of the beginning.

Winston Churchill

In this thesis we introduced a methodology based on information theory to identify functional dynamical structures in complex systems. To do so we modeled a system using a random field evolving in time. In particular, the dynamical behavior of each system component is described by means of a discrete stochastic process. We observed the system evolution in a given time interval and we estimated the probability distributions of the variables belonging to the random field. Then, exploiting information-theoretic measures (Dynamical Cluster Index and D-index), we looked for functional dynamical structures.

The method does not require any previous knowledge of underlying network topology of the complex system, but relies only on the values assumed by the random fields during the observation time interval.

The effectiveness of the methodology has been validated on test cases and subsequently applied to two prominent biological models, i.e. the mammalian cell cycle network and Mitogen Activated Protein Kinase (MAPK) cascade.

1. (a)

Random boolean networks. 2. (b)

Mammalian cell cycle networks. 3. (c)

Mitogen Activated Protein Kinase (MAPK) cascade.

5.1 Contributions and Novelty of the Work

The main novelty of the thesis, in comparison to previous application of the cluster index and of similar measures [35] is that we use it to consider truly dynamical systems, and not only fluctuations around stable asymptotic states. In principle, different kinds of data can be considered. In the case of a deterministic dynamical system, attractors are the main candidates to provide the required time series and we have shown in case (a) that they can work effectively. Note however that the method is ineffective in a situation that is sometimes encountered, i.e. if there is just a single fixed point. Therefore we conclude that the use of attractor states is only effective if the attractor landscape is rich enough to show the main features of the system organization. This has usually to be evaluated a posteriori, with the exception of trivial cases like that of a single fixed point.

It is worth noting that the three-level coding used in case (c) regards the similarities of the derivatives rather than those of the values of the variables and that the models that were used in this work to generate the time series are all based on first-order ordinary differential or difference equations; it should be verified whether the approach is valid also when higher-order dynamical systems are considered. Of course, high-order ODE systems can be transformed in first-order systems by adding variables, but the auxiliary variables that are required might turn out to be unobservable.

When a system is subject to continuous external disturbances the time series directly provide the required data and the experimental results show that our treatment can reveal its organizational features even when a high noise level is present.

Actually, the range of applicability of this method is quite broad and it does not necessarily need to be limited to dynamical system. Indeed the method just needs a set of frequencies of co-occurences of the values of the system variables. So the method can be applied also to many other systems, since all that is required is a series of ”cases” associated to vectors of numerical variables that are not necessarily ordered in time (think for example of different patients, each one described by a vector of values of various symptoms).

A final comment is that the method is not a brute-force one: when there is a clear organization in the system, e.g. in case (a), the organization can be read out directly from the order list of candidate relevant sets, but this does not always happen. The study of case (b) and (c) shows that even in entangled real systems the method provides useful clues to uncover the system organization. In order to describe the dynamics of these entangled systems, we introduce a sieving algorithm which selects disjoint or partially overlapping CRSs with the most high DCI values and we used the D-index in order to recover the direction of information flow among these CRSs.

5.2 Future Work

As far as future directions of this thesis are concerned, let us only mention some major ones. Several aspects have been mentioned above and will be subject to further analysis, including the effective analysis of high-order dynamical systems with a simplified discrete coding and the comparison among different kinds of time series data, i.e. attractors, transients from arbitrary initial conditions and perturbations.

Our method needs discrete variables, so when dealing with continuous data we have introduced a three-level coding that necessarily misses some information. Different coding schemes could be compared, moreover one might also consider the continuous generalization of the proposed method.

Other research will be devoted to improvements of the method: it is apparent that it faces a huge computational problem for large systems, since the number of possible subsets increases extremely fast with the number of variables. We used a genetic algorithm to efficiently search the optimal candidate relevant subsets in terms of Dynamical Cluster Index. However, the evaluation of the fitness function can be very expensive when high dimensional systems are considered. In order to deal with high dimensionality, a more efficient fitness function needs to be developed and other heuristic algorithms need to be tested.

Several other specific cases need to be studied to confirm the validity of the approach, mainly from real-world data. Let us mention that, among others, we are considering applications to the study of innovation processes, where data come from the real world and not from models. In this respect, it is important to realize that the DCI methodology may be integrated with other approaches: when dealing with real-world problems and not from data generated by a perfectly known model, a typical situation that is often encountered is that one knows some relationships between variables, but not them all. In this case, the most promising approach would be based on a combination of the a priori knowledge with the DCI, in ways that still need to be tested.

Moreover, it should be recalled that the Dynamical Cluster Index is just one out of several information-theoretic measures that might be applied to analyze dynamical systems. It is worth noting that the integration and the mutual Information can be useful even if used in isolation and not combined together in the DCI. But other measures, e.g. those that refer to joint distributions at different times, might prove to be particularly useful for the study of dynamical systems.

Finally, it is worth noting that the information exchange assessment is just preliminary; we are now working on a method to normalize and then to evaluate the significance of the information transmitted among CRS.