Exploring the organisation of complex systems through the dynamical interactions among their relevant subsets

Complex systems often show forms of organisation where a clear-cut hierarchy of levels with a well-defined direction of information flow cannot be found. In this paper we propose an information-theoretic method aimed at identifying the dynamically relevant parts of a system along with their relationships, interpreting in such a way the system’s dynamical organisation. The analysis is quite general and can be applied to many dynamical systems. We show here its application to two relevant biological examples, the case of mammalian cell cycle network and of Mitogen Activated Protein Kinase (MAPK) cascade. The result of our analysis shows that the elements of the mammalian cell cycle network act as a single compact group, whereas the MAPK system can be decomposed into two dynamically distinct parts, with asymmetric information flows.


Introduction
Complex systems often show forms of organisation that cannot be understood by referring to simple hierarchical models (like e.g. a tree).In most interesting cases one is faced with complicated situations, sometimes referred to as "tangled hierarchies", where a clear-cut hierarchy of levels, with a unique well-defined direction of information flow cannot be found.
In previous work (Villani et al., 2013;Villani and Serra, 2013;Villani et al., 2015) we have introduced the Dynamical Cluster Index method (shortly DCI1 ), which makes use of a measure based on information theory that appears well suited to explore the organisation of such complex systems.The DCI makes it possible to identify Candidate Relevant Subsets of variables (CRS) that show an integrated behaviour, while interacting more weakly with the rest of the system.The DCI has been applied with interesting results to several systems: some of them were artificially designed in order to test the effectiveness of the technique, while others refer to interesting chemical or biological systems (Filisetti et al., 2011;Sarma and Ghosh, 2012;Chaos et al., 2006;Farmer and Kauffman, 1986).Moreover, ongoing and yet unpublished work is also showing promising results in the study of social systems such as the detection of dynamical communities in regional innovation projects studied within the "Emergence by Design"2 European project, the Tuscany innovation system and the UK consumer price indices.
In many cases, though, the identified CRS have an intricate nested structure, so that it might not be clear which groups of variables are really important.To overcome this problem, a filtering algorithm is needed, which should be able to select those CRS responsible for the dynamics of the system.To this aim, in this paper we introduce a sieving procedure able to extract, from the set of candidate CRS, a smaller number of disjoint or partially overlapping sets, that capture the basic features of the groups of integrated variables, thus reducing the need for user intervention in the interpretation of the outcomes of the method.
In this work we also introduce a further improvement by considering the flow of information among the identified subsets of variables (i.e.relevant subsets RS) that is estimated on the basis of transfer entropy (Schreiber, 2000).In this way it is possible to determine whether there is an information flow from set A to set B, and vice-versa.Moreover, also the direction of the net information flow is defined.We use this method to analyse some artificially designed cases and also two biological models.
The paper is structured as in the following.In the subsequent sections we summarise the theoretical background and describe the sieving procedure and the technique used to detect the information flow among RS.Then, several case studies are presented, while in the final section conclusions and further developments are discussed.

The Dynamical Cluster Index method
In this section we summarise the main features of the DCI, as presented in (Villani et al., 2013(Villani et al., , 2015)).The DCI is an extension of the Functional Cluster Index (CI) introduced by Edelman and Tononi in 1998 (Tononi et al., 1998) and aimed at detecting functional clusters in brain regions.In our work, we relaxed the stationary constraint and extended the CI to actual dynamical systems, so as to apply it to a wide range of system classes, from abstract models to biological and socio-economic models.
The DCI is an information theoretical measure based on the well-know Shannon Entropy (Cover and Thomas, 2006).The entropy H(X) of a random variable X is computed as: The joint entropy of a pair of variables H(X, Y ) as: where Y assumes values y.Note that eq. 2 can be naturally extended to sets of k elements.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 measures based on the aforementioned entropy: the integration (I) and the Mutual Information (MI).I(S k ) measures the statistical independence of the k elements in S k and is defined as: Then mutual information M I(S k ; U \S k ) measures the mutual dependence between subset S k and the rest of the system U \S k and it is defined as: The value of DCI is not defined if M I(S k ; U \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 Villani et al. (2015) 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: 3 3 Introduced in (Tononi et al., 1998).
where DCI h and σ(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 ν = M I h / I h is the normalisation 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.
A list of CRS can be in principle obtained by computing 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 heuristic algorithms, such as genetic algorithms.However, since in this work we are interested in uncovering some properties concerning the information exchanged among RS, only small systems suitable for an exhaustive assessment will be tackled, postponing to a future contribution the presentation of the application of the DCI to large-size systems.Finally, random perturbations are applied to the dynamical system in order to observe their dynamical feedbacks also in case of stable situations.

Sieving the candidate subsets
The list of CRS can be ranked according to the significance of their DCI.In previous works we directly used this ranking to identify by hand the CRS relevant 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 ⊂ 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 organisation of the system.The pseudo-code of the algorithm is presented in Algorithm 1.

Assessing the temporal correlation among the subsets: the D index
Although by the application of the DCI, CRS are detected, this measure does not provide indications neither on the quantity nor on the direction of the entropy, i.e. the information flowing among subsets.To this aim we applied the directionality index proposed by Lautier and Raynaud (2014).Let X and Y be two random variables-or, equivalently, two sets of variables.We can define the entropy rate of X as the average number of bits needed to encode a successive Algorithm 1 Sieving algorithm Input: The array C of all the CRS ranked by their t c (DCI) 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 time t, eq.7a, or just on the value of X at time t, eq.7b.4 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. 8.
Thus, T Y →X can be described in term of entropy as: and, since the T Y →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.10.
Once that T Y →X and T X→Y are known, the directionality D of the information flow between X and Y can be measured as: where D X→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→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.
As introduced above, these considerations are also valid for groups of variables; however, the number of available observations limits the number of joint variables that may be analyzed (e.g.see Kraskov et al. (2004)): this is because spurious correlations are more frequently "detected" as this limit is approached, increasing the error in measurement.In order to avoid as much as possible this effect, each perturbation of the following analyses is introduced after the system has recovered a stable dynamical condition.Finally, the small size of the relevant subsets presented in this work allows a direct use of transfer entropy; in case of bigger subsets it is nevertheless possible to use sampling techniques, as in Lizier et al. (2011).

Analysis on test cases
In order to test the presented methodology, in the following we analyse several cases whose dynamical mechanisms are precisely known.We consider here boolean networks, a modelling framework that despite its apparent simplicity has obtained remarkable results in simulating real gene regulatory networks (Serra et al., 2004;Shmulevich et al., 2005;Serra et al., 2007;Villani et al., 2011).In particular, we present a system composed of 12 boolean nodes updated on the basis to either a boolean function or a random boolean value generator.Nodes update their state in parallel and synchronously.We illustrate the results of five instances of this network, defined in Table 1. 5he five 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 sieving 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 behaviour of (a part of) group B. We observe that the whole group B (Figure 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. 6  In case 4 five heterogeneously linked nodes are influenced by a triplet identical to that of groupA.The combination of the dynamical rules of the nodes and their initial condition makes the dynamical behaviour 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 the method detects their influence by returning some masks 7 having high t c values, partially 6 The interesting study of this effect is out of the scope of this paper and is subject of future work; in any case, a brief exploration of a limited set of initial conditions supports this preliminary claim (data not shown).
7 CRSs can be represented by rows, where entries corresponding 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 1b).
Case 5 derives from case 1 by adding two nodes whose dynamical behaviour 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 1a).

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 Fauré et al. (2006)-see Table 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 (Widmann et al., 1999).In Fauré et al. (2006) 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 (Fauré et al., 2006) 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 perturbed 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 Figure2a.
Alessandro Filisetti, Marco Villani, Andrea Roli, Marco Fiorucci, Roberto Serra (2015) Exploring the organisation of complex systems through the dynamical interactions among their relevant subsets.Proceedings of the European Conference on Artificial Life 2015, pp.286-293

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 (Sarma and Ghosh, 2012).Sarma and Ghosh 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 3a.The three chemicals identified as the core of this three-layered system are the M AP KKK, M AP KK and M AP K kinases (respectively M 3K , M 2K and MK for short) (Widmann et al., 1999).M 3K is activated by means of a single phosphorylation whereas M 2K 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 3a shows the schema of the M AP K cascade where each layer of the cascade is dephosphorylated by a specific phosphatase).Note that superimposed on the threelayered structure of substrates-product reactions there is the properly called M AP K signalling cascade, a linear chain of catalysis (dashed lines in Figure 3a) that transmit the external signal from M 3K * to MK * * .8When the external signal and the concentrations of the phosphatases are kept constant, a top-down reaction scheme as the one described in Figure 3a leads to fixed-point solutions.the other hand, oscillations have been reported in the MAPK cascade (Shin et al., 2009) and, in order to account for them, Sarma and Ghosh adopt a models with feedback, one of which is described in Figure 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 simulated the Sarma and Ghosh models with the CellDesigner software (Funahashi et al., 2008(Funahashi et al., , 2003)), 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 perturbed the asymptotic states of these models: in particular, we focused our analysis on kinases perturbations. 9The stable situations that are reached can  and 3, and 1 and 3. 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 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).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.

Conclusions
In this paper we have presented an improvement of a method based on information-theoretical measures able to identify the most relevant groups of variables that impact the dynamics of a system.We first introduced a sieving algorithm which makes it possible to identify the groups of variable that are responsible for system dynamics.Moreover, we show that it is possible to recover the direction of information flow among these groups, thus characterising the dynamical organisation of the system.
With respect to this, we noted that better normalisation methods can be applied to improve the algorithms efficacy: these methods are subject of ongoing work.
Finally, the effectiveness of this combination of techniques has been validated on test cases and subsequently applied to two prominent biological models, i.e. the mammalian cell cycle network and of Mitogen Activated Protein Kinase (MAPK) cascade.Alessandro Filisetti, Marco Villani, Andrea Roli, Marco Fiorucci, Roberto Serra (2015) Exploring the organisation of complex systems through the dynamical interactions among their subsets.Proceedings of the European Conference on Artificial Life 2015, pp.286-293  Alessandro Filisetti, Marco Villani, Andrea Roli, Marco Fiorucci, Roberto Serra (2015) Exploring the organisation of complex systems through the dynamical interactions among their relevant subsets.Proceedings of the European Conference on Artificial Life 2015, pp.286-293 Figure 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.

Figure 4 :
Figure 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 3a; M2 group involves the second layer, whereas groups M2 3 and M1 3 involve respectively the layers 2and 3, and 1 and 3.

Figure 2 :
Figure 2: 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.

Figure 3 :
Figure 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*).
Proceedings of the European Conference on Artificial Life 2015, pp.286-293