Presentation slides base on my paper about Complex System and Information Theory

Successful network inference from time-series data using mutual information rate E. Bianco-Martinez, N. Rubido, Ch. G. Antonopoulos, and M. S. Baptista Citation: Chaos 26, 043102 (2016); doi: 10.1063/1.4945420 View online: http://dx.doi.org/10.1063/1.4945420 View Table of Contents: http://scitation.aip.org/content/aip/journal/chaos/26/4?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Possible signatures of dissipation from time-series analysis techniques using a turbulent laboratory magnetohydrodynamic plasma Phys. Plasmas 23, 055709 (2016); 10.1063/1.4948275 Measurements of DNA barcode label separations in nanochannels from time-series data Biomicrofluidics 9, 064119 (2015); 10.1063/1.4938732 3-dimensional (orthogonal) structural complexity of time-series data using low-order moment analysis AIP Conf. Proc. 1479, 670 (2012); 10.1063/1.4756223 Detecting essential nodes in complex networks from measured noisy time series AIP Conf. Proc. 922, 603 (2007); 10.1063/1.2759751 Synchronization and control of spatiotemporal chaos using time-series data from local regions Chaos 8, 300 (1998); 10.1063/1.166310

Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 78.147.89.183 On: Sun, 08 May

2016 16:26:09

Successful network inference from time-series data using mutual information rate

E. Bianco-Martinez,1 N. Rubido,1,2 Ch. G. Antonopoulos,3 and M. S. Baptista1,a) 1Institute for Complex Systems and Mathematical Biology, SUPA, University of Aberdeen, King’s College, AB24 3UE Aberdeen, United Kingdom 2Instituto de F�ısica, Facultad de Ciencias, Universidad de la Rep�ublica, Igu�a 4225, Montevideo 11200, Uruguay 3Department of Mathematical Sciences, University of Essex, Wivenhoe Park, CO4 3SQ Colchester, United Kingdom

(Received 9 January 2016; accepted 23 March 2016; published online 12 April 2016)

This work uses an information-based methodology to infer the connectivity of complex systems

from observed time-series data. We first derive analytically an expression for the Mutual

Information Rate (MIR), namely, the amount of information exchanged per unit of time, that can

be used to estimate the MIR between two finite-length low-resolution noisy time-series, and then

apply it after a proper normalization for the identification of the connectivity structure of small net-

works of interacting dynamical systems. In particular, we show that our methodology successfully

infers the connectivity for heterogeneous networks, different time-series lengths or coupling

strengths, and even in the presence of additive noise. Finally, we show that our methodology based

on MIR successfully infers the connectivity of networks composed of nodes with different time-

scale dynamics, where inference based on Mutual Information fails. VC 2016 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4945420]

The Mutual Information Rate (MIR) measures the time

rate of information exchanged between two non-random

and correlated variables. Since variables in complex sys-

tems are not purely random, the MIR is an appropriate

quantity to access the amount of information exchanged

in complex systems. However, its calculation requires

infinitely long measurements with arbitrary resolution.

Having in mind that it is impossible to perform infinitely

long measurements with perfect accuracy, this work

shows how to estimate the MIR taking into consideration

this fundamental limitation and how to use it for the

characterization and understanding of dynamical and

complex systems. Moreover, we introduce a novel nor-

malized form of MIR that successfully infers the struc-

ture of small networks of interacting dynamical systems.

The proposed inference methodology is robust in the

presence of additive noise, different time-series lengths,

and heterogeneous node dynamics and coupling

strengths. Moreover, it also outperforms inference meth-

ods based on Mutual Information when analysing net-

works formed by nodes possessing different time-scales.

I. INTRODUCTION

We understand a complex system as a system with a

large number of interacting components whose aggregated

behaviour is non-linear and undetermined from the behav-

iour of the individual components. 1

If we now consider these

components as nodes of a network, and the underlying physi-

cal interaction between any two nodes as links, one way to

understand these complex systems is by studying its

topological structure, namely, the network connectivity. In

natural complex systems, the connectivity of the components

is often unknown or is difficult to detect by physical methods

due to large system-sizes. Hence, it is of interest to infer the

network structure that represents the physical interaction

between time-series collected from the dynamics of the

nodes.

Although network inference in non-linear systems has

been extensively studied in recent years using Cross-

Correlation or Mutual Information (MI), 2–4

recurrences, 5–7

functional dynamics, 8–11

and Granger Causality, 12–14

to

name a few, it still presents open challenges. The fundamen-

tal reason is that non-linearities, even in the absence of noise,

produce behaviour that hinders the correct identification of

existing or non-existing underlying direct physical depend-

ence between any pair of nodes.

In this paper, we introduce an information-based meth-

odology to infer the structure of complex systems from time-

series data. Our methodology is based on a normalized form

of an estimated MIR, the rate by which information is

exchanged per unit of time between any two components.

The MIR is an appropriate measure to quantify the exchange

of information in systems with correlation. 15–17

In particular,

the authors in Ref. 15 show how to calculate the MIR in the

case a Markov partition is attainable, which is generally

extremely difficult to find or unknown. Here, we first show

how the MIR can be approximately calculated for time-

series data of finite length and low-resolution. Then, we pro-

pose a normalization of the estimated MIR that allows for a

successful inference about the dependence structure of small

networks of interacting dynamical systems, when Markov

partitions are unknown. Our findings show that the estimated

normalized MIR allows for a successful inference of thea)Electronic mail: [email protected]

1054-1500/2016/26(4)/043102/9/$30.00 VC 2016 AIP Publishing LLC26, 043102-1

CHAOS 26, 043102 (2016)

Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 78.147.89.183 On: Sun, 08 May

2016 16:26:09

structure of small networks even in the presence of additive

noise, parameter heterogeneities, and different coupling

strengths. Moreover, our normalized estimated MIR outper-

forms the use of MI based inference when different time-

scale dynamics are present in the networks.

The paper is organized as follows. In Sec. II, we intro-

duce two information-based measures, the MI and the MIR.

We discuss the theoretical aspects of their definitions and

show how they are related to each other. In Sec. III, we intro-

duce the models used to create the complex system dynamics

studied in this work. In Sec. IV, we explain our methodology

to calculate an approximation value of MIR and introduce

the normalized MIR. Section V shows how we apply our

methodology to different coupled maps and to a neural net-

work in which the dynamics of the nodes is described by the

Hindmarsh-Rose (HR) neuron model. 18,19

Finally, in Sec.

VI, we discuss our work and discuss our findings.

II. BACKGROUND

Information can be produced in a system, and it can be

transferred between its different components. 15–17,20–23

If

transferred, at least two components that are physically inter-

acting by direct or indirect links should be involved. In gen-

eral, these components can be time-series, modes, or related

functions of them, defined on subspaces or projections of the

state space of the system. In this work, we study the amount

of information transferred per unit of time, i.e., the MIR,

between any two components of a system, to determine if a

link between them exists. The existence of a link between

two units means there is a bidirectional connection between

them due to their interaction.

A. Mutual information

The MI 24

between two random variables, X and Y, of a system is the amount of uncertainty one has about X (Y) after observing Y (X). Specifically, the MI is given by24–26

IXYðNÞ¼ HX þ HY � HXY; (1)

where HX ¼� PN

i¼1 PXðiÞ logðPXðiÞÞ and HY ¼� PN

j¼1 PYðjÞlogðPYðjÞÞ are the marginal entropies of X and Y (Shannon entropies), respectively, and HXY ¼�

PN2 i;j¼1

PXYði; jÞ logðPXYði; jÞÞ is the joint entropy between X and Y. PXðiÞ is the probability of a random event i to happen in X, PYðjÞ is the probability of a random event j to happen in Y, and PX;Yði; jÞ is the joint probability of events i and j to occur simultaneously in variables X and Y. N is the number of ran- dom events in both variables X and Y.

In particular, Eq. (1) can be written equivalently as

IXY Nð Þ¼ XN

i

XN j

PXY i; jð Þlog PXY i; jð Þ

PX ið ÞPY jð Þ

! : (2)

This equation can be interpreted as the strength of the de-

pendence between two random variables X and Y.25 When IXY¼0, the dependence strength between X and Y is null, consequently, X and Y are independent.

The computation of IXYðNÞ from time-series is a subtle task. First, it requires the calculation of probabilities

computed on an appropriate probabilistic space on which

a partition can be defined. Second, IXYðNÞ is a measure suitable for the comparison between pairs of components

of the same system but not between different systems. The

reason is that different systems can have different correla-

tion decay times, 27–29

hence, different characteristic time-

scales.

There are three main approaches to compute MI, and the

variation resides in the different ways to compute the proba-

bilities involved in Eq. (2). The first one is the bin or histo-

gram method, which finds a suitable partition of the 2D

space on equal or adaptive-size cells. 30,31

The second one

employs density kernels, where a kernel estimation of the

probability density function is used. 32,33

The last one com-

putes MI by estimating probabilities from the distances

between closest neighbours. 34

In this work, we adopt the first

method and compute probabilities in a partition of equally

sized cells in the probabilistic space generated by two varia-

bles X and Y. It is well known that this approach, proposed in Ref. 4 and studied in Ref. 35, overestimates the value of

IXYðNÞ for random systems or non-Markovian partitions.35,37 In particular, the authors explain two basic reasons for the

overestimation of MI: The finite resolution of a non-

Markovian partition and the finite length of the recorded

time-series. According to Refs. 35 and 37, these errors are

systematic and are always present in the computation of MI

for arbitrary non-Markovian partitions. Here, we avoid these

systematic errors by creating a novel normalization when

dealing with the MIR.

For the numerical computation of IXYðNÞ [Eq. (2)], we use the approach reported in Refs. 4 and 15. We define a

probabilistic space X, where X is formed by the time-series data observed from a pair of nodes, X and Y, of a complex system. Then, we partition X into a grid of N�N fixed-sized cells. The length-side of each cell, �, is then set to � ¼ 1=N. Consequently, the probability of having an event i for vari- able X, PXðiÞ, is the fraction of points found in row i of the partition X. Similarly, PYðjÞ is the fraction of points that are found in column j of X, and PXYði; jÞ is the joint probability computed from the fraction of points that are found in cell

(i, j) of the same partition, where i; j ¼ 1; …; N. We empha- size here that IXYðNÞ depends on the partition considered for its calculation as PX, PY, and PXY attain different values for different cell-sizes �.

B. Mutual information rate

Due to the issues arising from the definition of MI in

terms of its partition dependence, the authors in Ref. 15 have

demonstrated how to calculate the MIR for two time-series

of finite length irrespective of the partitions, instead of using

the MI. This quantity is invariant with respect to the resolu-

tion of the partition. 15

In particular, and for infinitely long

time-series, the MIR is theoretically defined as the mutual in-

formation exchanged per unit of time between X and Y.24,26,36 Specifically

043102-2 Bianco-Martinez et al. Chaos 26, 043102 (2016)

Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 78.147.89.183 On: Sun, 08 May

2016 16:26:09

MIRXY ¼ lim N!1

lim L!1

XL�1 i¼1

IXY i þ 1; Nð Þ� IXY i; Nð Þ L

¼ lim N!1

lim L!1

IXY L; Nð Þ� IXY 1; Nð Þ L

¼ lim N!1

lim L!1

IXY L; Nð Þ L

; (3)

where IXYðL; NÞ represents the MI of Eq. (1) between ran- dom variables X and Y, considering trajectories of length L that follow an itinerary over boxes in a grid with an infinite

number of cells N. Since IXY is a symmetric function with respect to X and Y, MIRXY ¼ MIRYX. We also note that the term

IXYð1;NÞ L

tends to zero in the limit of infinitely long trajec-

tories, L !1. The authors in Ref. 15 show that if a partition with N

cells is a Markov partition of order T, then the MIR can be estimated from finite-length and low-resolution time-series

(since the limits in Eq. (3) are not necessary) by using

MIRXY ¼ IXY Nð Þ T Nð Þ

; (4)

where both T(N) and N are finite quantities. Notice that an order T partition can only generate statistically significantly probabilities if there is in each cell a sufficiently large

amount of points (see Eq. (17)). Besides, points in a cell

must spread over the probabilistic space X after T iterations. So, the length of the time-series must be reasonably larger

than T. In Sec. IV, we make a novel demonstration of Eq. (4),

from which it becomes clear why MIR can be estimated

from finite-length and low-resolution time-series. In this

equation, IXYðNÞ is the MI between X and Y, considering probabilities that are calculated in a Markov partition, and

T(N) represents the shortest time for the correlation between X and Y to be lost for that particular Markov partition. T(N) also represents the time after which the evolution of a cha-

otic system is unpredictable. Moreover, this time is of the

order of the shortest Poincar�e return-time29 and is related to the order O Markov partition, where O indicates that the future state of a random variable X is independent on its ðO � 1Þ s previous states and is independent on the states of X for an order O�T .

III. MODELS FOR OUR COMPLEX SYSTEMS

We adopt various topologies for the networks and vari-

ous dynamics for the components of the complex systems

considered. Hence, the network inference, which represents

the detection of the topological structure of the component’s

interactions, is done from the time-series that are recorded

for each component. In particular, we divide the analysis on

discrete and on continues time-series components.

A. Networks with discrete-time units

The dynamics of the class of discrete complex systems

that are of interest here are described by the following

equation: 42

xinþ1 ¼ f x i n; r

� � 1 � að Þþ

a ki

XM j¼1

Aijf x j n; r

� � ; (5)

where xin is the n-th iterate of map i, where i ¼ 1; …; M and M is the number of maps (nodes) of the system, a 2 ½0; 1� is the coupling strength, Aij is the binary adjacency matrix

(with entries 1 or 0, depending on whether there is a connec-

tion between nodes i and j or not, respectively) that defines the structural connectivity in the network, r is the dynamical parameter of each map, ki ¼

PM j¼1 Aij is the node-degree,

and fðxn; rÞ is the considered map. Particularly, we use

fðxn; rÞ¼ rxnð1 � xnÞ; and (6)

f xn; rð Þ¼ xn þ r � K

2p sin 2pxnð Þ mod 1: (7)

For the logistic map 38–40

of Eq. (6), we use r¼4 (if it is not explicitly mentioned), that corresponds to fully developed

chaos, whereas for the circle map 41

of Eq. (7) we use

r¼0.35 and K � 6:9115, following Ref. 2, for the same reason.

Figure 1(a) shows the network topology described by

the adjacency matrix Aij used to create a network where the

dynamics of each node is described either by logistic or

circle maps. We will use these networks to study the robust-

ness of our methodology for different coupling strengths,

observational noise, and data-length. We also use small-size

networks with discrete dynamics, with different decay of

correlation times for the nodes to test our methodology (see

Fig. 1(b)). In those networks, the dynamics of the nodes is

given by logistic maps. In particular, we construct a network

formed by two clusters of 3 nodes each. The clusters are con-

nected by a small-coupling strength link. Specifically, the

dynamics of Fig. 1(b) for the cluster formed by the nodes 1,

2, and 3 is constructed by using r¼4, and the dynamics of the cluster formed by the nodes 4, 5, and 6 is given by a

third-order composition of the logistic map, i.e.,

fðxiÞ� f � f � fðxiÞ, with r¼3.9. Consequently, both clus- ters are constructed by time-series with different correlation

decay times, creating a good example to understand how a

clustered network with different time-scales can affect the

inference capabilities of MI- or MIR-based methodologies.

B. Networks with continuous-time units

We consider continuous dynamics for the nodes of a net-

work described by the HR neuron model. 18

The particular

network we choose is shown in Fig. 1(c). The HR model is

given by

_p ¼ q � ap3 þ bp2 � n þ Iext; _q ¼ c � dp2 � q; _n ¼ h½sðp � p0Þ� n�;

(8)

where p is the membrane potential, q is associated with the fast currents (Naþ or Kþ), and n with the slow current, for example, Ca2þ. The rest of the parameters are defined as a¼1, b¼3, c¼1, d¼5, s¼4, p0 ¼�1:6, and Iext ¼ 3:25,

043102-3 Bianco-Martinez et al. Chaos 26, 043102 (2016)

Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 78.147.89.183 On: Sun, 08 May

2016 16:26:09

for which the system exhibits a multi-scale chaotic behaviour

with spike bursting. The parameter h¼0.005 modulates the slow dynamics of the system. The neural networks of M neu- rons connected by electrical (linear coupling) synapses is

described in Refs. 19 and 23 and corresponds to having

_pi ¼ qi � ap3i þ bp 2 i � ni þ Iext � gl

XM j¼1

CijHðpjÞ;

_qi ¼ c � dp2i � qi; _ni ¼ h½sðpi � p0Þ� ni�; i ¼ 1; …; M;

(9)

where M is the number of neurons, and HðpiÞ¼ pi. In Eq. (9), gl is the strength of the electrical synapses. We use as initial conditions for each neuron i: pi ¼�1:30784489 þ gri ; qi ¼�7:32183132 þ gri ; ni ¼ 3:35299859 þ gri , and /i ¼ 0,

where gri is a uniformly distributed random number in ½0; 0:5� for all i ¼ 1; …; Nn, following Refs. 19 and 23. Cij is a Laplacian matrix and accounts for the way neurons are electri-

cally (diffusively) coupled. Particularly, Cij ¼ Kij � Aij, where A is the binary adjacency matrix of the electrical connections, and K is the node degrees’ diagonal matrix based on A. If Aði; jÞ¼ 1, then neuron j perturbs neuron i with an intensity given by gl.

IV. METHODS

A. Calculation of the correlation decay time using the diameter of an itinerary network

To infer the topology of a network using MIR [Eq. (4)],

we need to compute the correlation decay time T(N). T(N) is difficult to calculate in practical situations, since it depends on

quantities such as Lyapunov exponents and expansion rates,

which demand a high computational cost. 15

Here, we estimate

it by the number of iterations that takes to points in cells of X to expand and completely cover X. This is a necessary condi- tion to determine the shortest time for the correlation to decay

to zero. In particular, we are introducing a novel way to calcu-

late T(N) from the diameter of a network G, which is based on the dynamics of points mapped from one cell of X to another, namely, a network with the connectivity given by the transi-

tions of points from cell to cell of X or an itinerary network. We construct G as follows. We assume that each equally

sized cell in X, occupied by at least one point, represents a node in G. Then, following the dynamics of points moving from one cell to another, we create the connections between

nodes, i.e., the links in G. Specifically, a link between nodes i and j exists if points in X travel from cell i to cell j. If the link exists the weight is equal to 1, if it is absent, then it is

equal to 0, therefore, G is defined as a binary matrix with ele- ments Gij 2f0; 1g. In this framework, a uniformly random time-series with no correlation results in a complete network,

namely, an all-to-all network.

We define T(N) as the diameter of G. The reason is that T(N) is the minimum time that takes for points inside any cell of X to spread to the whole extent of X. By definition, the diameter of a network is the maximum length for all

shortest-paths, i.e., the minimum distance required to cross

the entire network. Hence, our approach transforms the cal-

culation of T(N) into the calculation of the diameter of G. In particular, for the estimation of the network diameter, we use

the Johnson’s algorithm. 43–48

B. Calculation of MIR

To estimate the MIR from finite-length low-resolution

time-series data, we truncate the summation in Eq. (3) up to

a finite size N, depending on the resolution of data, and con- sider small trajectory pieces of the time-series with a length

L, which depends on the total length of the time-series and on Eq. (17), such that

MIRXY ffi 1

L

XL i¼1

IXY i þ 1; Nð Þ� IXY i; Nð Þ � �

: (10)

FIG. 1. Network topologies used to construct the complex systems. Panel

(a) shows a network with 16 nodes and with similar characteristics with a

scale-free network, where the dynamics of each node is either a logistic or a

circle map. Panel (b) shows a network composed of 2 clusters of 3 nodes

each, which is composed of nodes with different time-scales and a logistic

map dynamics. Panel (c) shows a network of 12 nodes, where the dynamics

of each node is described by the Hindmarsh-Rose dynamics [Eq. (8)].

043102-4 Bianco-Martinez et al. Chaos 26, 043102 (2016)

Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 78.147.89.183 On: Sun, 08 May

2016 16:26:09

In Eq. (10), left-hand and right-hand sides would be equal if

the partition, where probabilities are being calculated, is

Markov. The length L represents also the largest order T that a partition that generates statistically significant probabilities

can be constructed from these many trajectory pieces.

Assuming that the order of the partition constructed is T¼L (which also represents the time for the correlation in the par-

tition to decay to zero, if the partition would be Markov),

then Eq. (10) becomes

MIRXY ffi 1

T

XT i¼1

IXY i þ 1; Nð Þ� IXY i; Nð Þ � �

: (11)

Now, taking two partitions, K1 and K2, with different corre- lation decay times, T1 and T2, respectively, and different number of cells, N1 � N1 and N2 � N2, respectively, with N2 > N1, we have T2 ¼ T1 þ 1. Moreover, K1 generates K2 in the sense that F�1ðK1Þ¼ K2, where F is the evolution op- erator, and F�1ðK1Þ means the pre-iteration of partition K1. Then

IXYðT2; K1Þ¼ IXYðT1; K2Þ: (12)

Hence, we can write Eq. (11) as

MIRXY ffi 1

T1

XT1 i¼1

IXY i þ 1; K1ð Þ� IXY i; K1ð Þ½ �

ffi 1

T1

XT1 i¼1

IXY i; K2ð Þ� IXY i; K1ð Þ½ �: (13)

When the partition is a Markov generating partition, its

properties 15

fulfil

IXYði; KkÞ¼ IXYð1; Kkþi�1Þ: (14)

Then, if our partition is close to a Markov partition, Eq. (11)

results in

MIRXY ’ 1

T1 IXY 1; KT1þ1ð Þ� IXY 1; K1ð Þ � �

(15)

� 1

T1 IXY 1; KT1ð Þ; (16)

which is our demonstration for the validity of Eq. (4).

Therefore, in order to use Eq. (15), we must have parti-

tions for which Eq. (14) is approximately valid. This condi-

tion can be reached for partitions constructed with a

sufficiently large number of equally sized cells of length

� ¼ 1=N, exactly the type of partition considered here. Notice, however, that the partitions will typically not be

Markov nor generating, causing systematic errors in the esti-

mation of MIR. To correct these errors, we propose the nor-

malizations in Eqs. (18) and (19).

It is important to notice that MIRXY is always a

partition-independent quantity, if and only if, the partitions

are Markov. In order to calculate IXYð1; KT1Þ, we use Eq. (1), which requires the calculation of probabilities in X. Fulfilling the inequality

hN0ðNÞi Noc; (17)

where hN0ðNminÞi is the mean number of points inside all occupied cells of the partition of X, Eq. (17) guarantees that the probabilities are unbiased.

C. Network inference using MIR

For our analysis, using a non-Markovian partition allows

us to simplify the calculations of MIRXY , however, taking

this kind of partitions into consideration would make the

MIR values to oscillate around an expected value. Moreover,

the MIR for different non-Markovian partitions not only has

a non-trivial dependence with the number of cells in the par-

tition but also presents a systematic error. 35

Therefore, since

MIRXY for a non-Markovian partition of N�N equally sized cells [estimated by Eq. (4)] is expected to be partition-

dependent, we propose here a way to obtain a measure, com-

puted from MIRXYðNÞ, that is partition independent and that is suitable for network inference.

To infer the structure of a network, we calculate the MIR

for the MðM � 1Þ=2 different pairs of nodes in the network, which is all we need due to the symmetric property of MIR.

We also discard the MIR values for the same variable, i.e.,

MIRXX, because we are interested in the exchange of informa-

tion between different variables. We compute the MIRXY exchanged between any two nodes in a network by taking the

expected value over different partition sizes Ni, i.e., MIRXY ¼ EiðMIRXYðNiÞÞ, where E(X) is the expected value of X. In order to remove the systematic error35 in this calculation, we perform instead a weighted average, where the finer parti-

tions (larger N) contribute more to the MIRXY value than the coarser ones (smaller N). The reason is that a smaller N is likely to create a partition that is further away from a Markovian one

than a partition of larger N. Consequently, we resolve the sys- tematic error by weighing differently the different partitions.

Therefore, we propose a novel normalization for the

MIR as follows. First, we use an equally sized grid of size N, we subtract from MIRXYðNÞ, calculated for all pairs of nodes, its minimum value and denote the new quantity as

minðMIRXYðNÞÞ. Theoretically, a pair that is disconnected should have a MIR value close to zero; however, in practice,

the situation is different because of the systematic errors

coming from the use of a non-Markovian partition, as well

as from the information flow passing through all the nodes in

the network. For example, the effects of a perturbation in

one single node will arrive to any other node in a finite

amount of time. This subtraction is proposed to reduce these

two undesired overestimations of MIR. After this step, we

remain with MIR as a function of N. Normalizing then by maxðMIRXYðNÞÞ� minðMIRXYðNÞÞ, where again the maxi- mum and minimum are taken over all different pairs, we

construct a relative magnitude M̂IRXYðNÞ, namely,

M̂IRXY Nð Þ¼ MIRXY Nð Þ� minfMIRXY Nð Þg

maxfMIRXY Nð Þg� minfMIRXY Nð Þg ; (18)

where MIRXYðNÞ is the MIR between nodes X and Y, minfMIRXYðNÞg is the minimum with respect to the MðM � 1Þ=2 pairs, and maxfMIRXYðNÞg is the maximum

043102-5 Bianco-Martinez et al. Chaos 26, 043102 (2016)

Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 78.147.89.183 On: Sun, 08 May

2016 16:26:09

with respect to all MðM � 1Þ=2 pairs. This magnitude is still a function of N; however, we can now perform an average over different values of N without the systematic error.

Next, we apply Eq. (18) for different grids sizes Ni; i ¼ 1; …; m to obtain MIRXYðNiÞ, where Nm is the maximum number of cells per axis, resulting in a grid of Nm � Nm cells, and fulfilling at the same time Eq. (17). Then, similarly to

the idea used for Eq. (18), we make a second normalization

over M̂IRXYðNiÞ to obtain

MIRXY ¼ P

iM̂IRXY Nið Þ maxf

P iM̂IRXY Nið Þg

; (19)

where the maximum is being taken now over the Nm grids. Finally, applying Eq. (19) to each pair XY, we obtain its

average value, MIRXY . The higher the value of MIRXY , the

higher the amount of information exchanged between X and Y per unit of time. This allows us to identify pairs of nodes that exchange larger rates of information than others.

In order to perform the network inference from the MIR,

we fix a threshold in ½0; 1� and create a binary adjacency ma- trix Ac, where the entry AcX;Y is 1 if MIRXY is higher than the threshold, and 0 otherwise. Ac is then compared with the adjacency matrix A used to construct the dynamics of the nodes in Sec. III. Recording the threshold used to create Ac, and varying it in ½0; 1�, we obtain different inferred networks. Our results show that there is an interval of thresholds within

½0; 1� that fulfil Ac ¼ A, i.e., a band that represents a 100% successful network inference.

In general, the effectiveness of our network inference

methodology is measured by the absolute difference between

the real topology and the one inferred for different threshold

values. We find that whenever there is a band of threshold

values, there is successful inference without errors. In practi-

cal situations, where the underlying network is unknown and

the absolute difference is impossible to compute, the ordered

values of the MIR or other similarity measures 2,3

show a pla-

teau which corresponds to the band of thresholds aforemen-

tioned. In particular, if the plateau is small, the authors in

Ref. 49 propose a method to increase the size of the plateau

by “silencing” the indirect connections, hence, allowing for

a more robust reconstruction of the underlying network.

V. RESULTS FOR NETWORK INFERENCE

We now present our results for network inference using

the three models introduced in Sec. III.

A. Discrete-time systems

1. Different coupling strengths

Here, we study the performance of Eq. (19) for network

inference in the case where the dynamics of each node is

described by a circle or a logistic map. The network structure

that comprises our small-network of interacting discrete-

time systems is given in Fig. 1(a). Here, we analyze the

effectiveness of the inference as the coupling strength, a, between connected nodes is varied. In Ref. 2, the authors

have shown that, for the logistic [Eq. (6)] and circle maps

[Eq. (7)], assuming the same topology, the dynamics is

quasi-periodic for a > 0:15 and chaotic for 0 a 0:15. We, therefore, choose the coupling strength a in Eq. (5) to be equal to 0.03 and 0.12, both values corresponding to chaotic

dynamics.

Figure 2 shows the network inference results using

MIRXY . The wideness of the red band represents all possible

values a threshold can take to perform a 100% success net-

work inference, i.e., the correct identification of all physical

and non-physical links. The wider the band, the bigger the

probability to perform a complete reconstruction, therefore

the reconstruction is more robust. When we deal with experi-

mental data, and the correct topology is unknown, the opti-

mal threshold can be determined by the range of consecutive

thresholds for which the inferred topology is invariant, see

Ref. 2.

An error in the percentage of reconstruction comes from

links that were not inferred (false negatives) or inferred erro-

neously (false positives). In our current study, we avoid the

distinction between them, and we categorize both as recon-

struction errors. Then, the reconstruction percentage can

decrease by inferring non-existent links (non physical links)

or by missing them. Each time this happens, we decrease the

percentage by an amount e% ¼ 100 1 Nl

, where Nl is the num- ber of real links in the original network.

2. Different time-series lengths and noise strengths

We start by analysing the effectiveness of MIRXY for

different time-series lengths, using the dynamics of the logis-

tic map for each node and a coupling strength a 2 ½0; 0:17�. In Fig. 3(a), we observe that for a closer to 0.15, a relatively short length (of about 3000 points) is enough to infer cor-

rectly the original network, which is generated by the adja-

cency matrix A of Sec. III. However, when a is close to 0.03, a larger time-series (of about 30 000 points) is needed to

FIG. 2. Network inference for different coupling strengths and coupled maps. Panels (a) and (b) represent the MIRXY values between different pair of nodes in

a network composed of coupled logistic maps with coupling strengths � ¼ 0:03 and 0.12, respectively. Panels (c) and (d) are similar to panels (a) and (b), but for circle maps. The red band indicates the range of thresholds from which the original network is correctly inferred, namely, achieving 100% successful

inference.

043102-6 Bianco-Martinez et al. Chaos 26, 043102 (2016)

Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 78.147.89.183 On: Sun, 08 May

2016 16:26:09

achieve 100% successful reconstruction. The values of a¼0 and a 2 ½0:15; 0:17� are considered to test the effectiveness of the values MIRXY in the case of nodes being totally inde-

pendent and in the case of nodes having periodic dynamics.

In these regions, MIRXY is expected to be zero, a situation

evidenced in both panels of the figure. Our results so far sug-

gest that the successful reconstruction for short-length time-

series depends on the intensity of the coupling strength.

However, it is surprising to see that exact inference can

always be achieved for this dynamical regime if a suffi-

ciently large time-series is available.

Next, we apply our methodology for network inference

considering noisy time-series data. In particular, we intro-

duce additive normally distributed noise to the logistic map,

i.e.,

fnsðxn; rÞ¼ fðxn; rÞþ c � r; (20)

where fnsðxn; rÞ is the noisy dynamics, r is a random number drawn from the normal distribution with 0 mean and stand-

ard deviation of 1, i.e., @ð0; 1Þ, and c 2 ½0; 1� is the noise strength. Since @ðc; 1Þ¼ c �@ð0; 1Þ, the noise strength is the

standard deviation in the normal distribution. Fig. 3(b) shows

the parameter space for different coupling strengths versus c. We observe perfect inference for noise strengths c < 0:3, i.e., for @ð0; 1Þ. Moreover, the best reconstruction using MIRXY is for coupling strengths in ½0:6; 0:11�, a dynamical regime where chaotic behaviour is prevalent.

B. Neural networks

We also apply our methodology for the study of network

inference in the case of continuous dynamics given by the HR

system. We use two electrical couplings, gl ¼ 0:05 and 0.1, both considered for time-series of length 2 � 105. Figure 4 shows the band for 100% successful network inference, where

panel (a) corresponds to gl ¼ 0:05 and panel (b) to gl ¼ 0:1. This figure shows that MIRXY is able to infer the correct net-

work structure, in this case, for small networks of continuous-

time interacting components.

C. Comparison between mutual information and mutual information rate

Finally, we compare MI and MIRXY to assess the effec-

tiveness of our proposed methodology for network inference.

We apply the same normalization process used for MIR, Eq.

(19), to MI to have an appropriate comparison. In particular,

we infer the network structure of the system described in

Sec. III with the network shown in Fig. 1(b). As we have

explained in Sec. III, this system has two clusters of nodes

with different dynamics. The dynamics in the left cluster is

given by the 3rd-order composition of the logistic map,

whereas the dynamics of the right cluster is given by ordi-

nary logistic map dynamics. The different dynamics of the

two groups produces different correlation decay times, T(N), for nodes X and Y, in particular, when the pair of nodes comes from different clusters. The different correlation

decay times produce a non-trivial dynamical behaviour that

challenges the MI performance for network inference.

Figure 5 shows the results obtained for the normalized

MI, �I XY , and our normalized MIR, MIRXY , for each of the possible pairs of nodes. The purple bars correspond to the

pairs of nodes 1, 2, and 3 of the first cluster, the orange bars

correspond to the pairs of nodes 4, 5, and 6 of the second

cluster (3rd order composed dynamics), and the black bar

corresponds to the link between clusters (notice that due to

the small coupling strength between the two clusters this link

is not detected using any of the two methods). Nevertheless,

FIG. 3. Network inference based on logistic maps, for different coupling

and noise strengths. Panel (a) shows the parameter space of the percentage

of reconstruction (0%—blue and 100%—dark red) for different coupling

strengths versus data lengths. Panel (b) is similar to panel (a) for different

coupling strengths versus the standard deviation of the normal distribution

of the noise added.

FIG. 4. Network inference for a network of nodes with a HR neural dynam-

ics for different electrical couplings. Panels (a) and (b) show the bar plots of

the percentage of inference for gl ¼ 0:05 and gl ¼ 0:1, respectively. The red bands show the range of thresholds for which the original network is inferred

with a 100% success.

043102-7 Bianco-Martinez et al. Chaos 26, 043102 (2016)

Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 78.147.89.183 On: Sun, 08 May

2016 16:26:09

the MIR identifies correctly all intra links of the network

where MI fails to do so. We conclude that the normalized

MIR is preferable over the normalized MI when it comes to

the detection of links in a complex system with different cor-

relation decay times. The reason is that the normalized MIR

takes into consideration the correlation decay time associated

to each pair of nodes, contrary to the MI.

VI. CONCLUSIONS

In this paper, we have introduced a new information

based approach to infer the network structure of complex

systems. The MIR is an information measure that computes

the information transferred per unit of time between pairs of

components in a complex system. MIRXY , our novel normal-

ization for the MIR that is introduced in Eq. (18), is a mea-

sure based on MIR and developed for network inference. We

find that MIRXY is a robust measure to perform network in-

ference in the presence of additive noise, short time-series,

and also for systems with different coupling strengths. Since

MIR and MIRXY depend on the correlation decay time T, they are suitable for inferring the correct topology of net-

works with different time-scales.

In particular, we have explored the effectiveness of MIR

versus MI in terms of how successful they are in inferring

exactly the network of our small complex systems. In gen-

eral, we find that the MIR outperforms the MI when different

time-scales are present in the system. Our results also show

that both measures are sufficiently robust and reliable to infer

the networks analyzed whenever a single time-scale is pres-

ent. In other words, small variations in the dynamical param-

eters, time-series length, noise intensity, or topology

structure maintain a successful inference for both methods. It

remains to be seen the types of errors that are found in these

measures when perfect inference is missing or impossible to

be done.

ACKNOWLEDGMENTS

E.B.M., M.S.B., and C.G.A. acknowledge the financial

support provided by the EPSRC “EP/I032606/1” grant.

C.G.A. contributed to this work while working at the

University of Aberdeen and then, while working at the

University of Essex, United Kingdom. N.R. acknowledges

the support of PEDECIBA, Uruguay.

1 A. Barrat, M. Barthelemy, and A. Vespignani, Dynamical Processes on Complex Networks (Cambridge University Press, 2008).

2 N. Rubido, A. C. Mart�ı, E. Bianco-Martinez, C. Grebogi, M. S. Baptista, and C. Masoller, “Exact detection of direct links in networks of interacting

dynamical units,” New J. Phys. 16, 093010 (2014). 3 G. Tirabassi, R. Sevilla-Escoboza, J. M. Buld�u, and C. Masoller, “Inferring the connectivity of coupled oscillators from time-series statisti-

cal similarity analysis,” Sci. Rep. 5, 10829 (2015). 4 A. Butte and I. S. Kohane, “Mutual information relevance networks:

Functional genomic clustering using pairwise entropy measurements,”

Pac. Symp. Biocomput. 5, 415–426 (2000). 5 Y. Zou, M. C. Romano, M. Thiel, N. Marwan, and J. Kurths, “Inferring indi-

rect coupling by means of recurrences,” Chaos 21(4), 1099–1111 (2011). 6 J. Nawrath, M. C. Romano, M. Thiel, I. Z. Kiss, M. Wickramasinghe, J.

Timmer, J. Kurths, and B. Schelter, “Distinguishing direct from indirect

interactions in oscillatory networks with multiple time scales,” Phys. Rev.

Lett. 104, 038701 (2010). 7 S. Hempel, A. Koseska, J. Kurths, and Z. Nikoloski, “Inner composition

alignment for inferring directed networks from short time series,” Phys.

Rev. Lett. 107, 054101 (2011). 8 M. Timme, “Revealing network connectivity from response dynamics,”

Phys. Rev. Lett. 98, 224101 (2007). 9 S. G. Shandilya and M. Timme, “Inferring network topology from com-

plex dynamics,” New J. Phys. 13, 013004 (2011). 10

M. Timme and J. Casadiego, “Revealing networks from dynamics: An

introduction,” J. Phys. A: Math. Theor. 47, 343001 (2014). 11

X. H. Ta, Ch. N. Yoon, L. Holm, and S. K. Han, “Inferring the physical

connectivity of complex networks from their functional dynamics,” BMC

Syst. Biol. 4(70), 1–12 (2010). 12

C. W. J. Granger, “Investigating causal relations by econometric models

and cross-spectral methods,” Econometrica 37, 424–438 (1969). 13

S. Guo, Ch. Ladroue, and J. Feng, “Granger causality: Theory and

applications,” in Frontiers in Computational and Systems Biology (Springer, London, 2010), Vol. 15, Chap. 5, pp. 83–111.

14 S. L. Bressler and A. K. Seth, “Wiener–Granger causality: A well estab-

lished methodology,” Neuroimage 58(2), 323–329 (2010). 15

M. S. Baptista, R. M. Rubinger, E. R. Viana, J. C. Sartorelli, U. Parlitz,

and C. Grebogi, “Mutual information rate and bounds for it,” PLoS One

7(10), e46745 (2012). 16

M. Palus, “Coarse-grained entropy rates for characterization of complex

time series,” Physica D 93, 64–77 (1996). 17

T. Schreiber, “Measuring information transfer,” Phys. Rev. Lett. 85,

461–464 (2000). 18

J. L. Hindmarsh and R. M. Rose, “A model of neuronal bursting using

three coupled first order differential equations,” Proc. R. Soc. London B

221, 87–102 (1984). 19

M. S. Baptista, F. M. Kakmeni, and C. Grebogi, “Combined effect of chemi-

cal and electrical synapses in Hindmarsh-Rose neural networks on synchro-

nization and the rate of information.,” Phys. Rev. E 82(3), 036203 (2010). 20

M. A. Machiori, R. Fariello, and M. A. M. de Aguiar, “Energy transfer dy-

namics and thermalization of two oscillators interacting via chaos,” Phys.

Rev. E 85, 041119 (2012).

FIG. 5. Network inference for the composed dynamical system introduced

in Sec. III A. Panel (a) plots IXY of Eq. (1) for all links. This is a case where a complete network inference cannot be achieved (indicated by the absence

of any red band). Panel (b) is the same as before but for MIRXY . The color

code corresponds to the same color code identifying different nodes in Fig.

1(b). The darkest color is the link connecting the two clusters.

043102-8 Bianco-Martinez et al. Chaos 26, 043102 (2016)

Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 78.147.89.183 On: Sun, 08 May

2016 16:26:09

21 D. Mandal, H. T. Quan, and C. Jarzynski, “Maxwell’s refrigerator: An

exactly solvable model,” Phys. Rev. Lett. 111(3), 030602 (2013). 22

Ch. G. Antonopoulos, E. Bianco-Martinez, and M. S. Baptista,

“Production and transfer of energy and information in Hamiltonian sys-

tems,” PLoS One 9(2), e89585 (2014). 23

Ch. G. Antonopoulos, S. Srivastava, S. E. de, S. Pinto, and M. S. Baptista,

“Do brain networks evolve by maximizing their information flow capaci-

ty?,” PLoS Comput. Biol. 11(8), e1004372 (2015). 24

C. E. Shannon, “A mathematical theory of communication,” Bell Syst.

Tech. J. 27, 379–423, 623–656 (1948). 25

S. Kullback, Information Theory and Statistics (Wiley, New York, 1959). 26

R. L. Dobrushin, “General formulation of Shannon’s main theorem of in-

formation theory,” Usp. Mat. Nauk. 14, 3–104 (1959) [Am. Math. Soc.

Transl. 33(2), 323–438 (1959)]. 27

J. P. Eckmann and D. Ruelle, “Ergodic theory of chaos and strange

attractors,” Rev. Mod. Phys. 57, 617–656 (1985). 28

M. S. Baptista, J. X. De Carvalho, and M. S. Hussein, “Finding quasi-

optimal network topologies for information transmission in active

networks,” PLoS One 3(10), e3479 (2008). 29

M. S. Baptista, N. Eulalie, P. R. F. Pinto, M. Brito, and J. Kurth, “Density

of first Poincare returns, periodic orbits, and Kolmogorov-Sinai entropy,”

Phys. Lett. A 374, 1135–1140 (2010). 30

G. A. Darbellay and I. Vajda, “Estimation of the information by an adapt-

ive partitioning of the observation space,” IEEE Trans. Inf. Theory 45(4),

1315–1321 (1999). 31

A. M. Fraser and H. L. Swinney, “Independent coordinates for strange

attractor from mutual information,” Phys. Rev. A 33(2), 1134 (1986). 32

Y. Moon, B. Rajagopalan, and U. Lall, “Estimation of mutual information

using kernel density estimators,” Phys. Rev. E 52(3), 2318 (1995). 33

A. Gretton, R. Herbrich, A. Smola, O. Bousquet, and B. Sch€olkopf, “Kernel methods for measuring independence,” J. Mach. Learn. Res. 6,

2075–2129 (2005). 34

A. Krakov, H. St€ogbauer, and P. Grassberger, “Estimating mutual information,” Phys. Rev. E 69, 066138 (2004).

35 R. Steuer, J. Kurths, C. O. Daub, J. Weise, and J. Selbig, “The mutual in-

formation: Detecting and evaluating dependencies between variables,”

Bioinformatics 18(2), 231–240 (2002).

36 R. M. Gray and J. C. Kieffer, “Asymptotically mean stationary measures,”

IEEE Trans. Inf. Theory 26, 412–421 (1980). 37

H. Herzel, A. O. Schmitt, and W. Ebeling, “Finite sample effects in

sequence analysis,” Chaos, Solitons Fractals 4(1), 97–113 (1994). 38

R. M. May, “Simple mathematical models with very complicated dynami-

cs,” Nature 261(5560), 459–467 (1976). 39

G. Benettin, L. Galgani, A. Giorgilli, and J. M. Strelcyn, “Lyapunov characteris-

tic exponents for smooth dynamical systems and for Hamiltonian systems; a

method for computing all of them. Part 1: Theory,” Meccanica 15, 9–20 (1980). 40

G. Benettin, L. Galgani, A. Giorgilli, and J. M. Strelcyn, “Lyapunov char-

acteristic exponents for smooth dynamical systems and for Hamiltonian

systems; a method for computing all of them. Part 2: Numerical

application,” Meccanica 15, 21–30 (1980). 41

P. Collett and J. P. Eckmann, Iterated Maps on the Interval as Dynamical Systems (Springer Science & Business Media, 2009).

42 K. Kaneko, “Clustering, coding, switching, hierarchical ordering, and con-

trol in a network of chaotic elements,” Physica D 41(2), 137–172 (1990). 43

J. Garc�ıa-Ojalvo, M. B. Elowitz, and S. H. Strogatz, “Modeling a synthetic multicellular clock: Repressilators coupled by quorum sensing,” Proc.

Natl. Acad. Soc. 101(30), 10955–10960 (2004). 44

Ch. Skokos, “The Lyapunov characteristic exponents and their

computation,” in Dynamics of Small Solar System Bodies and Exoplanets (Springer, Berlin, Heidelberg, 2010), pp. 63–135.

45 K. Wendel, N. G. Narra, M. Hannula, J. Hyttinen, and J. Malmivuo, “The

influence of electrode size on EEG lead field sensitivity distributions,” Int.

J. Bioelectromagn. 9(2), 116–117 (2007). 46

J. M. Stinnett-Donnelly, N. Thompson, N. Habel, V. Petrov-Kondratov, D.

D. Correa de Sa, J. H. Bates, and P. S. Spector, “Effects of electrode size

and spacing on the resolution of intra-cardiac electrograms,” Coron.

Artery Dis. 23(2), 126–132 (2012). 47

N. K. Chen, C. C. Dickey, S. S. Yoo, C. R. Guttmann, and L. P. Panych,

“Selection of voxel size and slice orientation for fMRI in the presence of

susceptibility field gradients: Application to imaging of the amygdala,”

Neuroimage 19(3), 817–825 (2003). 48

D. B. Johnson, “Efficient algorithms for shortest paths in sparse networks,”

J. ACM 24(1), 1–13 (1997). 49

B. Barzel and A.-L. Barab�asi, Nat. Biotechnol. 31(8), 720–725 (2013).

043102-9 Bianco-Martinez et al. Chaos 26, 043102 (2016)

Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 78.147.89.183 On: Sun, 08 May

2016 16:26:09