Topological Information Data Analysis

TL;DR

This paper introduces topological information data analysis methods to quantify complex statistical interactions in data, revealing insights into independence, co-variation, and clustering, with applications in genetics and biology.

Contribution

It develops new topological and information-theoretic tools for analyzing multivariate data, especially in biological systems, addressing estimation biases and demonstrating practical applications.

Findings

Vanishing I_k indicates statistical independence.

Maximal I_k identifies co-varying variables.

Application to genetic data reveals biologically relevant subtypes.

Abstract

This paper presents methods that quantify the structure of statistical interactions within a given data set, and was first used in \cite{Tapia2018}. It establishes new results on the k-multivariate mutual-informations (I_k) inspired by the topological formulation of Information introduced in. In particular we show that the vanishing of all I_k for 2\leq k \leq n of n random variables is equivalent to their statistical independence. Pursuing the work of Hu Kuo Ting and Te Sun Han, we show that information functions provide co-ordinates for binary variables, and that they are analytically independent on the probability simplex for any set of finite variables. The maximal positive I_k identifies the variables that co-vary the most in the population, whereas the minimal negative I_k identifies synergistic clusters and the variables that differentiate-segregate the most the population. Finite data size effects and estimation biases severely constrain the effective computation of the information topology on data, and we provide simple statistical tests for the undersampling bias and the k-dependences following. We give an example of application of these methods to genetic expression and unsupervised cell-type classification. The methods unravel biologically relevant subtypes, with a sample size of 41 genes and with few errors. It establishes generic basic methods to quantify the epigenetic information storage and a unified epigenetic unsupervised learning formalism. We propose that higher-order statistical interactions and non identically distributed variables are constitutive characteristics of biological systems that should be estimated in order to unravel their significant statistical structure and diversity. The topological information data analysis presented here allows to precisely estimate this higher-order structure characteristic of biological systems.

Equations

$$H(P_X)=-\sum _{x\in E}{P}_X(x){\log }_2{P}_X(x).$$

$$H_X(Y;P)=H(Y_{\ast }(P)),$$

$$Y_{\ast }(P)(y)=\sum _{x|Y(x)=y}P(x).$$

$$\forall x|Y(x)=y,$$

$$\forall x|Y(x)\ne y,$$

$$Y.H(Z;P)=\sum _{y\in E_Y}H(Z;P|(Y=y))Y_{\ast }(P)(y).$$

$$(W,Y).H(Z;P)=W.(Y.H)(Z;P).$$

$$H((Y,Z);P)=H(Y;P)+Y.H(Z;P).$$

$$P_X\mapsto F_X(Y_1;...;Y_k;P_X)$$

$$\forall X^{\prime },X\le X^{\prime },\forall P_{X^{\prime }},F_X(Y_1;...;Y_k;X_{\ast }(P_{X^{\prime }}))=F_{X^{\prime }}(Y_1;...;Y_k;P_{X^{\prime }}).$$

$$Y.F_X(Y_1;...;Y_k;P)=\sum _{y\in E_Y}{F}_X(Y_1;...;Y_k;P|(Y=y))Y_{\ast }(P)(y).$$

$$H_k(X_1;...;X_k;P)=H((X_1,...,X_k);P).$$

$$H_j(X_1;...;X_k;P)=\sum _{I\subset [k];card(I)=j}H(X_I;P),$$

$$I_k(X_1;...;X_k;P)=\sum _{j=1}^{j=k}(-1{)}^{j-1}{H}_j(X_1;...;X_k;P),$$

$$H(X_1;X_2;...;X_n;\mathbb{P})=\sum _{k=1}^{k=n}(-1{)}^{k-1}\sum _{I\subset [n];card(I)=k}{I}_k(X_{i_1};X_{i_2};...;X_{i_k};\mathbb{P})$$

$$I_{k,l}(Y_1;...;Y_k;P_X|Z_1,...,Z_l)=(Z_1,...,Z_l).I_k(Y_1;...;Y_k;P_X).$$

$$I_{k,l}(X_{I_1};...;X_{I_k};(\mathbb{P}|X_J))=X_J.I_k(X_{I_1};...;X_{I_k};\mathbb{P}),$$

$$I_{k,l}(X_{I_1};...;X_{I_k};(\mathbb{P}|X_J))=\phi ({\Sigma }_{I_1}\cap ...\cap {\Sigma }_{I_k}\backslash {\Sigma }_J),$$

$$\begin{gathered} H_{k+1}(X_0;X_1;...;X_k;\mathbb{P})=H_k((X_0,X_1);X_2;...;X_k;\mathbb{P}) \\ =H_k(X_1;...;X_k;\mathbb{P})+X_0.H_k(X_1;...;X_k;\mathbb{P}). \end{gathered}$$

$$I_k((X_0,X_1);X_2;...;X_k;\mathbb{P})=I_k(X_0;X_2;...;X_k;\mathbb{P})+X_0.I_k(X_1;X_2;...;X_k;\mathbb{P})$$

$$I_{k+1}(X_0;X_1;...;X_k;\mathbb{P})=I_k(X_1;X_2;...;X_k;\mathbb{P})-X_0.I_k(X_1;X_2;...;X_k;\mathbb{P}).$$

$$I_2(X;(Y,Z))\ne I_2(X;Y)+I_2(X;Z),$$

$$I_2(X;(Y,Z))+I_2(Y;Z)=I_2((X,Y);Z)+I_2(X;Y),$$

$$(A\cap (B\cup C))\cup (B\cap C)=(A\cap B)\cup ((A\cup B)\cap C).$$

$${\eta }_I(P)=\eta (X_{i_1};...;X_{i_k};P)=X_{[n]\backslash I}.I_k(X_{i_1};...;X_{i_k};P),$$

$$H_k(X_{i_1};...;X_{i_k};P)=\sum _{J\subset [n]|\exists m,i_m\in J}{\eta }_J(P),$$

$$I_k(X_{i_1};...;X_{i_k};P)=\sum _{J\supset I}{\eta }_J(P).$$

$${\eta }_J(P)=\sum _{k\ge l}\sum _{I\supseteq J|dimI=k}(-1{)}^{k-l}{H}_k(X_I;P).$$

$$I_k(X_{i_1};X_{j_1},...;X_{j_{k-2}};X_{i_n})=I_2(X_{i_1};X_{i_n}).$$

$$G_k(X_1;...;X_k;P)=\sum _{i=1}^{k}H(X_i)-H(X_1;...;X_k).$$

$$H(X_1,...,X_n;P)=H(X_1;P)+...+H(X_n;P).$$

Full text

”When you use the word information, you should rather use the word form”

René Thom

1. 1 Introduction

2. 2 Results

3. 2.1 Entropy and Information functions

4. 2.2 The independence criterion

5. 2.3 Information coordinates

6. 2.4 Application to Gene expression data - detection of cell types and gene modules

7. 2.5 Mutual-information negativity, clusters and links

8. 2.6 Cell type detection - comparison with previous MaxEnt studies

9. 3 Discussion

10. 4 Material and Methods

11. 4.1 The dataset: quantified genetic expression in two cell types

12. 4.2 Probability estimation

13. 4.3 Computation of k-Entropy, k-Information landscapes and paths

14. 4.4 Estimation of the undersampling dimension: statistical result

15. 4.5 Estimation of the undersampling dimension: Computational result

16. 4.6 k-dependence test

17. 4.6.1 Sampling size and graining landscapes - stability of minimum energy complex estimation

18. A Appendix: Bayes free energy and Information quantities

19. A.1 Parametric modelling

20. A.2 Bethe approximation

1 Introduction

This note presents a method of statistical analysis of a set of collected characters in a population, describing a kind of topology of the distribution of information in the data. New theoretical results are developed to justify the method. The data that concern us are represented by certain (observed or computed) parameters $s_{1},...,s_{n}$ belonging to certain finite sets $E_{1},...,E_{n}$ of respective cardinalities $N_{1},...,N_{n}$, which depend on an element $z$ of a certain set $Z$, representing the tested population, of cardinality $N_{Z}$. In other terms we are looking at $N$ ”experimental” functions $X_{i}:Z\rightarrow E_{i},i=1,...,n$, then we will refer to the data by the letters $(Z,X)$, where $X$ is the product function of the $X_{i}$, going from $Z$ to the product $E$ of all the sets $E_{i},i=1,...,n$.

For instance, as in $[..]$, each $E_{i},i=1,...,n$ has cardinality $8$ and is identified with the subset of integers $[8]=\{1,...,8\}$, each $s_{i}=S_{i}(z),i=1,...,n$ measures the level of expression of a gene $g_{i}$ in a neuron $z$ belonging to a set $Z$ of classified dopaminergic neurons (DA); another system to be compared to this one is made by the analog measurements for neurons $z^{\prime}$ belonging to a set $Z^{\prime}$ of classified non-dopaminergic neurons (NDA). To be precise, in this example, $n=41$, $N_{Z}=111$, $N_{Z^{\prime}}=37$.

The most usual hypothesis for interpreting the data is the existence of an objective and immutable joint probability $\mathbb{P}_{X}$ on the set $E$, coming from a hypothetical set $\Omega$ enclosing $Z$ and governing the results of the sample $(Z,X)$. Then, without any other knowledge, assuming that different experimental $z$ are independently chosen, the strong law of large numbers tells that the better possible approximation of $\mathbb{P}_{X}(s)$ is given by the number of $z$ such that $s=X(z)$ divided by the cardinality of $Z$. However, the standard inequalities of probability theory (Bienaymé-Tchebicheff, Markov, Chernov Kolmogorov, …) show that the confidence that we can have in this approximation depends heavily on $\mathbb{P}_{X}$ itself. This generates a risk of circularity.

Several approaches can be followed to escape circularity; part of them, maintaining the frequentist point of view, uses the Fisher information metric (cf. [34]), part of them, using a Bayesian approach, puts probability laws on the set of probabilities themselves, and studies the evolution of these choices with the introduction of new data. Once the total joint probability is found, it is theoretically possible to verify its agreement with the marginal laws on $X_{i},i=1,...,n$, or any partial joint variable $X_{I}=(X_{i_{1}},...,X_{i_{k}})$ for $I=\{i_{1},...,i_{k}\}\subset[n]$, but practically it is not so easy, being one of the principal problems in Statistical Mechanics or in Bayesian Analysis.

We will follow here a different approach, which consists in describing the manner the variables $X_{i},i=1,...,n$ distribute the Information on $(Z,X)$. The experimented population $Z$ has its own characteristics that the data explore, and the frequency of every value $s_{I}$ of each one of the variables $X_{I},I\subset[n]$ is an information important by itself, without considering the hypothetical law on the whole set $E$. The information quantities, derived from the Shannon entropy, offer a natural way for describing all these frequencies. In fact they define the form of the distribution of information contained in the raw data. For instance, the individual entropies $H(X_{i}),i=1,...,n$ tell us the shape of the individual variables: if $H(X_{i})$ is small (with respect to its capacity $\log_{2}N_{i}$), then $X_{i}$ corresponds to a well defined characteristic of $Z$; to the contrary if $H(X_{j})$ is close to the capacity, that is the value of the entropy of the uniform distribution, the function $X_{j}$ corresponds to a non-trivial partition of $Z$, and does not correspond to a well defined invariant. At the second degree, we can consider the entropies $H(X_{i},X_{j})$ for every pair $(i,j)$, giving the same kind of structures as before, but for pairs of variables. To get a better description of this second degree with respect to the first one, we can look at the mutual information as defined by Shannon, $I(X_{i};X_{j})=H(X_{i})+H(X_{j})-H(X_{i},X_{j})$. If it is small, near zero, the variables are not far from being independent, if it is maximal, i.e. not far from the minimum of $H(X_{i})$ and $H(X_{j})$, this means that one of the variables is almost determined by the other. In fact $I(X_{i};X_{j})$ can be taken as a measure of dependence, due to its universality and its invariance. Consider the graph with vertices $X_{i},i=1,...,N$ and edges $(X_{i},X_{j}),i\neq j,i,j=1,...,N$; by labeling each vertex with the entropy and each edge with the mutual information, we get a sort of one-dimensional skeleton of the data $(Z,X)$. The information of higher degrees define in an analogous manner the higher dimensional skeletons of the data $(Z,X)$ (see figure 3 for example).

Works of Clausius, Boltzmann, Gibbs and Helmholtz underlined the importance of entropy and free energy in Statistical Physics. In particular, Gibbs gave the general definition of the entropy for the distribution of microstates, cf. [18]. Later Shannon recognized in this entropy the basis of Information theory in his celebrated work on the mathematical theory of communication [52] (equation 11), and then further developed their structure in the lattice of variables [51]. Note that this kind of lattice takes its roots in the work of Boole on Logic and Probability [8]. Defining the communication channel, information transmission and its capacity, Shannon also introduced to degree two (pairwise) mutual information functions [52].

The expression and study of multivariate higher degree mutual-informations (equation 12) was achieved in two seemingly independent works: 1) McGill (1954) [38] (see also Fano (1961) [16]) with a statistical approach, who called these functions ”interaction information”, and 2) Hu Kuo Ting (1962) [23] with an algebraic approach who also first proved the possible negativity of mutual-informations for degrees higher than $2$. The study of these functions was then pursued by Te Sun Han [21, 22].

Higher-order mutual-informations were then rediscovered in several different contexts, notably by Matsuda in 2001 in the context of spin glasses, who showed that negativity is the signature of frustrated states [36] and by Bell in the context of Neuroscience, Dependent Component Analysis and Generalised Belief Propagation on hypergraphs [6]. Brenner and colleagues have observed and quantified an equivalent definition of negativity of the $3$ variables mutual information, noted $I_{3}$, in the spiking activity of neurons and called it synergy [11]. Anastassiou and colleagues unraveled $I_{3}$ negativity within gene expression, corresponding in that case to cooperativity in gene regulation [66, 30].

Another important family of information functions, named ”total correlation”, which corresponds to the difference between the sum of the entropies and the entropy of the joint, was introduced by Watanabe in 1960 [65]. These functions were also rediscovered several times, notably by Tononi and Edelman who called them ”integrated information” [61] in the context of consciousness quantification, and by Studený and Vejnarova [56] who called them ”multi-information” in the context of graphs and conditional independences. Closely related to what we present here with respect to both the kind of data analyzed and the conclusions, Margolin and colleagues [35] used these functions of Watanabe to quantify higher order statistical dependences within genetic expression.

In our approach, for any data $(Z,S)$, the full picture can be represented by a collection of numerical functions on the faces of a simplex $\Delta([n])$ having vertices corresponding to the random variables $X_{1},...,X_{n}$. We decided to focus on two subfamilies of Information functions: the first is the collection of entropies of the joint variables, denoted $H_{k},k=1,...,n$, giving the numbers $H_{k}(X_{i_{1}};...;X_{i_{k}})$, and the degree $k$ information of the joint variables, denoted $I_{k},k=1,...,n$, and giving the numbers $I_{k}(X_{i_{1}};...;X_{i_{k}})$; see the following section for their definition and their elementary properties. In particular, the value on each face of a given dimension of these functions gives interesting curves (histograms, see section on statistics 2.2) for testing the departure from independence, and their means over all dimensions for testing the departure from uniformity of the variables. These functions are information co-chains of degree $k$ (in the sense of ref [4]) and have nice probabilistic interpretations. By varying in all possible manners the ordering of the variables, i.e. by applying all the permutations $\sigma$ of $[n]=\{1,...,n\}$, we obtain $n!$ paths $H_{k}(\sigma)$, $I_{k}(\sigma)$, $k=1,...,n$. They constitute respectively the $H_{k}$-landscape and the $I_{k}$-landscape of the data.

When the data correspond to uniform and independent variables, that is the uninteresting null hypothesis, each path is monotonic, the $H_{k}$ growing linearly and the $I_{k}$ being equal to zero for $k$ between $2$ and $n$. Any departure from this behavior (estimated for instance in Bayesian probability on the allowed parameters) gives a hint of the form of information in the particular data.

Especially interesting are the maximal paths, where $I_{k}(\sigma)$ decreases, being strictly positive, or strictly negative after $k=3$. Other kinds of paths could also be interesting, for instance the paths with the maximal total variation as they can be oscillatory. In the examples provided here and in [58], we proposed to stop the empirical exploration of the information paths to their first minima, a condition of vanishing of conditional mutual-informational (conditional independence).

As a preliminary illustration of the potential interest of such functions for general Topological Data Analysis, we quantify the information structures for the empirical measures of the expression of several genes in two pre-identified populations of cells presented in [58], and we consider here both cases where genes or cells are considered as variables for gene or cell unsupervised classification tasks respectively.

In practice, the cardinality $N_{Z}$ of $Z$ is rather small with respect to the number of free parameters of the possible probability laws on $E$, that is $N-1=N_{1}...N_{n}-1$, then the quantities $H_{k},I_{k}$ for $k$ larger than a certain $k_{u}$ have in general no meaning, a phenomenon commonly called undersampling or curse of dimensionality. In the example, $n$ is $20$, but $k_{u}$ is $11$. Moreover, the permutations $\sigma$ of the variables values can be applied to test the estimation of the dependences quantified by the $I_{k}$ against the null hypothesis of randomly generated statistical dependences. In this approach describing the raw data for themselves, undersampling is not a serious limitation. However, it is better to test the stability of the shape of the landscapes by studying random subsets of $Z$. Moreover, the analytic properties of $H_{k}$ and $I_{k}$ considered as functions of $P$ in a given face of the simplex of probabilities $\Delta([n])$ ensure that, if $P_{X}$ tends to $\mathbb{P}$ in this face, the shape is preserved.

The originality of our method is the systematic consideration of the entropy and the information landscapes and paths that can be associated to all possible permutations of the basic variables, and the extraction of exceptional paths from them, in order to define the overall form of the distribution of information among the set of variables. This new perspective has its origin in the local (topos) homological theory introduced in [4]. Moreover this method was successfully applied to a concrete problem of gene expression in [5, 58].

In the present article, we first remind the definitions and basic properties of the entropy and information chains and functions. We give equivalent formulations of the Hu Kuo Ting theorem [23], which allows to express every partial mutual conditioned higher information of collections of joint variables from elementary higher entropies $H_{k}(X_{I})$ or by elementary higher mutual information functions $I_{k}(X_{I})$, i.e. the functions that form the entropy landscape and information landscape, respectively.

Second we establish that these ”pure” functions are analytically independent as functions of the probability laws, in the interior of the large simplex $\Delta([n])$. This follows from the fact we also prove here, that these functions constitute coordinates (up to a finite ambiguity) on $\Delta([n])$ in the special case of binary variables $X_{i},i=1,...,n$. In addition, we demonstrate that, for every set of numbers $N_{i},i=1,...,n$, the cancellation of the functions $I_{k}(X_{I}),k\geq 2,I\subset[n]=\{1,...,n\}$ is a necessary and sufficient condition of the set of variables $X_{1},...,X_{n}$ to be statistically independent. We were not able to find these results in the literature. They generalize results of Te Sun Han [21, 22].

Then this article not only presents a method of analysis but it gives proofs of basic results on information quantities that, to our knowledge, were not available until now in the literature.

Third we study the statistical properties of the entropy and information landscapes and paths, and present the computational aspects. The mentioned examples of genetic expression are developed. Finally in an appendix, we show how these functions appear in the theory of Free energies, in Statistical Physics and in Bayesian Variational Analysis.

2 Results

2.1 Entropy and Information functions

Given a probability law $P_{X}$ on a finite set $E=E_{X}$, Shannon defined the information content of this law by the Boltzmann-Gibbs entropy [52]:

[TABLE]

Shannon himself gave an axiomatic justification of this choice, that was developed further by Khinchin, Kendall and other mathematicians, see [29].

The article [4] (Baudot and Bennequin) presented such a set of axioms inspired by algebraic topology, see also [63] (Juan-Pablo Vigneaux). In all these approaches, the fundamental ingredient is the decomposition of the entropy for the joint variable of two variables. To better formulate this decomposition, we have proposed to consider the entropy as a function of three variables: first a finite set $E_{X}$, second a probability law $P$ on $E_{X}$ and third, a random variable on $E_{X}$, i.e. a surjective map $Y:E_{X}\rightarrow E_{Y}$, considered only through the partition of $E_{X}$ that it gives, indexed by the elements $y$ of $E_{Y}$. In this case we say that $Y$ is less fine than $X$, and write $Y\leq X$, or $X\rightarrow Y$. Then we define the entropy of $Y$ for $P$ at $X$:

[TABLE]

where $Y_{*}(P)$ is the image law, also named the marginal of $P$ by $Y$:

[TABLE]

Remark**.**

Frequently, when the context is clear, we simply write $H_{X}(Y;P)=H(Y;P)$ or even $H(Y)$, as everybody does, however the ”homological nature” of $H$ can only be understood with the index $X$, because it is here that the topos theory appears, see [4, 63].

The second fundamental operation on probabilities (after marginalization) is the conditioning: given $y\in E_{Y}$, such that $Y_{*}(P)(y)\neq 0$, the conditional probability $P|(Y=y)$ on $E_{X}$ is defined by the following rules:

[TABLE]

This allows to define the conditional entropy, as Shannon has done, for any $Z$ and $Y$ both less fine that $X$,

[TABLE]

Note that if $P|(Y=y)$ is not well defined, we can simply forget it in the formula, because it appears multiplied by zero.

This operation is associative (see [4, 63]), i.e. for any triple $W,Y,Z$ of variables less fine than $X$,

[TABLE]

With these notations, the fundamental functional equation of Information Theory, or its first axiom, according to Shannon, is

[TABLE]

Remark**.**

In [4, 63] it is shown that this equation can be understood as a co-cycle equation of degree one of a module in a topos, in the sense of Grothendieck and Verdier [1], and why the entropy is generically the only universal generator of the first co-homology functor.

More generally, we consider a collection of sets $E_{X},X\in\mathcal{C}$, such that each time $Y,Z$ are less fine than $X$ and belong to $\mathcal{C}$, then $(Y,Z)$ also belongs to $\mathcal{C}$; in this case we name $\mathcal{C}$ an information category. An example is given by the joint variables $X=(X_{i_{1}},...,X_{i_{m}})$ of $n$ basic variables $X_{1},...,X_{n}$ with values in finite sets $E_{1},...,E_{n}$; the set $E_{X}$ being the product $E_{i_{1}}\times...\times E_{i_{m}}$.

Then for every natural integer $k\geq 1$, we can consider families indexed by $X$ of (measurable) functions of the probability $P_{X}$ that are indexed by several variables $Y_{1},...,Y_{k}$ less fine than $X$

[TABLE]

satisfying the compatibility equations;

[TABLE]

We call these functions the co-chains of degree $k$ of C for the probability laws. An equivalent axiom is that $F_{X}(Y_{1};...;Y_{k};P_{X})$ only depends on the image of $P_{X}$ by the joint variable $(Y_{1},...,Y_{k})$. We call this property locality of the family $F=(F_{X},X\in\mathcal{C})$.

The action by conditioning extends verbally to the co-chains of any degree:

if $Y$ is less fine than $X$,

[TABLE]

It satisfies again the associativity condition.

Higher mutual information quantities were defined by Hu Kuo Ting [23] and McGill [38], generalizing the Shannon mutual information [4, 58]:

in our terms, for $k$ random variables $X_{1},...,X_{k}$ less fine than $X$ and one probability law $P$ on the set $E_{X}$,

[TABLE]

And more generally, for $j\leq k$, we define

[TABLE]

where $X_{I}$ denotes the joint variable of the $X_{i}$ such that $i\in I$.

We name these functions of $P$ the joint entropies.

Then the higher information functions are defined by

[TABLE]

In particular we have $I_{1}=H$, the usual entropy.

Reciprocally the functions $I_{k}$ decompose the entropy of the finest joint partition:

[TABLE]

The following result is immediate from the definitions, and the fact that $H_{X},X\in\mathcal{C}$ is local:

Proposition 1**.**

*The joint entropies $H_{k}$ and the higher information quantities $I_{k}$ are information co-chains, i.e. they are local functions of $P$.

Remark**.**

From the computational point of view, locality is important, because it means that only the less fine marginal probability has to be taken into account.

The definition of $H_{j},j\leq k$ and $I_{k}$ makes evident that they are symmetric functions, i.e. they are invariant by every permutation of the letters $X_{1},...,X_{k}$.

The particular case $I_{2}(S;T)=H(S)+H(T)-H(S,T)$ is the usual mutual information defined by Shannon.

Using the concavity of the logarithm, it is easy to show that $I_{1}$ and $I_{2}$ have only positive values, but this ceases to be true for $I_{k}$ as soon as $k\geq 3$ [23, 36].

Hu kuo Ting defined in [23] other information quantities, by the following formulas:

[TABLE]

For instance, considering a family of basic variables $X_{i},i=1,...,n$,

[TABLE]

for the joint variables $X_{I_{1}},...,X_{I_{k}},X_{J}$, where $I_{1},...,I_{k},J\subset[n]$.

The following remarkable result is due to Hu Kuo Ting [23]:

Theorem 2.1**.**

Let $X_{1},...,X_{n}$ be any set of random variables and $\mathbb{P}$ a given probability on the product $E_{X}$ of the respective images $E_{1},...,E_{n}$, then there exist finite sets $\Sigma_{1},...,\Sigma_{n}$ and a numerical function $\varphi$ from the union $\Sigma$ of these sets to $\mathbb{R}$, such that for any collection of subsets $I_{m};m=1,...,k$ of $\{1,...,n\}$, and any subset $J$ of $\{1,...,n\}$ of cardinality $l$, the following identity holds true

[TABLE]

where we have denoted $X_{I}=(X_{i_{1}},...,X_{i_{l}})$ and $\Sigma_{I}=\Sigma_{i_{1}}\cup...\cup\Sigma_{i_{l}}$ for $I=\{i_{1},...,i_{l}\}$, and where $\Omega\backslash\Sigma_{J}$ denotes the set of points in $\Omega$ that do not belong to $\Sigma_{J}$, i.e. the set $\Omega\cap(\Sigma\backslash\Sigma_{J})$, named subtraction of $Y=\Sigma_{J}$ from $\Omega$.

The Hu Kuo Ting theorem says that for a given joint probability law $\mathbb{P}$, and from the point of view of the information quantities $I_{k,l}$, the joint operation of variables corresponds to the union of sets, the graduation $k$ corresponds to the intersection, and the conditioning by a variable corresponds to the difference of sets. This can be precisely formulated as follows:

Corollary 1.1**.**

Let $X_{1},...,X_{n}$ be any set of random variables on the product $E_{X}$ of the respective goals $E_{1},...,E_{n}$, then for any probability $\mathbb{P}$ on $E_{X}$, every universal identity between disjoint sums of subsets of a finite set that are obtained, starting with $n$ subsets $\Sigma_{1},...,\Sigma_{n}$, by 1) forming collections of reunions, 2) taking successive intersections of these unions, and 3) subtracting by one of them, gives an identity between sums of information quantities, by replacing the union by the joint variables $(.,.)$, the intersections by the juxtaposition $(.;.;.)$ and the subtraction by the conditioning.

Remark**.**

Conversely the corollary implies the Theorem.

This corollary is the source of many identities between the information quantities.

For instance, the fundamental equation (6) corresponds to the fact that the union of two sets $A,B$ is the disjoint union of one of them, say $A$ and of the difference of the other with this one, say $B\backslash A$.

The following formula follows from (6):

[TABLE]

The two following identities are also easy consequences of the Corollary $1$; they are important for the method of data analysis presented in this article:

Proposition 2**.**

Let $k$ be any integer

[TABLE]

Proposition 3**.**

Let $k$ be any integer

[TABLE]

Remark**.**

Be careful that some universal formulas between sets do not give identities between information functions; for instance $A\cap(B\cup C)=(A\cap B)\cup(A\cap C)$ but in general we have

[TABLE]

What is true is the following identity:

[TABLE]

This corresponds to the following universal formula between sets

[TABLE]

*The formula (21) follows directly from the definition of $I_{2}$, by developing the four terms of the equation. It expresses the fact that $I_{2}$ is a simplicial co-cycle, being the simplicial co-boundary of $H$ itself.

However, although this formula between sets is true, it is not of the form authorized by the corollary $1$.

Consequently, some identities of sets that are not contained in the Theorem 2.1 correspond to information identities, but, as we saw just before with the false formula (20), not all identities of sets correspond to information identities.*

As we already said, the set of joint variables $X_{I}$, for all the subsets $I$ of $[n]=\{1,...,n\}$, is an information category, the set $\mathcal{C}$ being the $n-1$-simplex $\Delta([n])$ of vertices $X_{1},...,X_{n}$. In what follows we do not consider more general information categories.

We can paraphrase the Theorem 2.1, by a combinatorial Theorem on the simplex $\Delta([n])$:

Definition 4**.**

Let $X_{1},...,X_{n}$ be a set of random variables with respective goals $E_{1},...,E_{n}$, and let $X_{I}=\{X_{i_{1}},...,X_{i_{k}}\}$ be a face of $\Delta([n])$, we define, for a probability $P$ on the product $E$ of all the $E_{i},i=1,...,n$,

[TABLE]

Remark**.**

With the exception $J=[n]$, the function $\eta_{J}$ is not an information co-chain of degree $k$. But it is useful in the demonstrations of some of the following results.

Embed $\Delta([n])$ in the hyperplane $x_{1}+...+x_{n}=1$ as the standard simplex in $\mathbb{R}^{n}$ (the intersection of the above hyperplane with the positive cone, where $\forall i=1,...,n,x_{i}\geq 0$), and consider the balls $\Sigma_{1},...,\Sigma_{n}$ of radius $R$ strictly larger than $\sqrt{(n-1)/n}$ that are centered on the vertices $X_{j};j=1,...,n$; they have all possible non-empty intersections convex. The subsets $\Sigma^{\prime}_{I}=\Sigma_{I}\setminus\Sigma_{[n]\backslash I}$ are the connected components of complementary set of the unions of the boundary spheres $\partial\Sigma_{1},...,\partial\Sigma_{n}$ in the total union $\Sigma$ of the balls $\Sigma_{1},...,\Sigma_{n}$.

Proposition 5**.**

*For every $k+1$ subsets $I_{1},..,I_{k},K$ of $[n]$, if $l$ denotes the cardinality of $K$, the information function $I_{k,l}(X_{I_{1}};...;X_{I_{k}};P|X_{K})$ is equal to the sum of the functions $\eta_{J}(P)$, where $J$ describes all the faces such that $\Sigma^{\prime}_{J}$ is one of the connected components of the set $(\Sigma_{I_{1}}\cap...\cap\Sigma_{I_{k}})\backslash\Sigma_{K}$.

Proof.

Every subset that is obtained from the $\Sigma_{J};J\subset[n]$ by union, intersection and difference, repeated indefinitely (i.e. every element of the Boolean algebra generated by the $\Sigma_{i};i=1,...,n$), is a disjoint union of some of the sets $\Sigma^{\prime}_{J}$. This is true in particular for the sets obtained by the succession of operations $1,2,3$ in the order prescribed by the corollary 1.1 above. Then the proposition follows from the corollary 1.1. ∎

We define the elementary (or pure) joint entropies $H_{k}(X_{I})$ and the elementary (or pure) higher information functions $I_{k}(X_{I})$ as $H_{k}(X_{i_{1}};...;X_{i_{k}};P)$ and $I_{k}(X_{i_{1}};...;X_{i_{k}};P)$ respectively, where $I=\{i_{1},...,i_{k}\}\subset[n]$ describes the subsets of $[n]$. In the following of the text, we will consider only these pure quantities. We will frequently denote them simply by $H_{k}$ (resp. $I_{k}$). The other information quantities use joint variables and conditioning, but the preceding result tells that they can be computed from the pure quantities.

For the pure functions, the decompositions in the basis $\eta_{I}$ are simple:

Proposition 6**.**

If $I=\{i_{1},...,i_{k}\}$, we have

[TABLE]

and

[TABLE]

In other terms, the function $H_{k}$ evaluated on a face $X_{I}$ of dimension $k$ is given by the sum of the functions $\eta_{J}$ over all the faces $X_{J}$ connected to $X_{I}$. And the function $I_{k}$ evaluated on $X_{I}$ is the sum of the functions $\eta_{J}$ over all the faces $X_{J}$ that contain $X_{I}$.

Proposition 7**.**

For any face $J$ of $\Delta([n])$, of dimension $l$, and any probability $P$ on $E_{X}$, we have

[TABLE]

Proof.

This follows from the Moebius inversion formula (Rota 1964 [47]). ∎

Corollary 7.1**.**

(Te Sun Han): Any Shannon information quantity is a linear combination of the pure functions $I_{k},k\geq 1$ (resp. $H_{k}k\geq 1$), with coefficients in $\mathbb{Z}$, the ring of relative integers.

Proof.

This follows from the proposition 5. ∎

Hu Kuo Ting [23] also proved a remarkable property of the information functions associated to a Markov process:

Proposition 8**.**

The variables $X_{1},...,X_{n}$ can be arranged in a Markov process $(X_{i_{1}},...,X_{i_{n}})$ if and only if, for every subset $J=\{j_{1},...,j_{k-2}\}$ of $\{i_{2},...,i_{n-1}\}$ of cardinality $k-2$, we have

[TABLE]

This implies that, for a Markov process between $(X_{i_{1}},...,X_{i_{n}})$, all the functions $I_{k}(X_{I})$ involving $i_{1}$ and $i_{n}$, are positive.

2.2 The independence criterion

The total correlations were defined by Watanabe as the difference of the sum of entropies and the joint entropy, noted $G_{k}$ [65] (see also [61, 56, 35]):

[TABLE]

Total correlations are Kullback-Leibler divergences, cf. appendix A on bayes free energy; and $I_{2}=G_{2}$. It is well known (cf. the above references or [13]) that for $n\geq 2$, the variables $X_{1},...,X_{n}$ are statistically independent for the probability $P$, if and only if $G_{n}(X_{1};...;X_{n})=0$, i.e.

[TABLE]

Remark**.**

The result is proved by induction using repetitively the case $n=2$, which comes from the strict concavity of the function $H(P)$ on the simplex $\Delta([n])$.

Theorem 2.2**.**

For every $n$ and every set $E_{1},...,E_{n}$ of respective cardinalities $N_{1},...,N_{n}$, the probability $P$ renders the $n$ variables $X_{i},i=1,...,n$ statistically independent if and only if the $2^{n}-n-1$ quantities $I_{k}$ for $k\geq 2$ are equal to zero.

Proof.

For $n=2$ this results immediately from the above criterion and the definition of $I_{2}$. Then we proceed by recurrence on $n$, and assuming that the result is true for $n-1$ we deduce it for $n$.

The definition of $I_{n}$ is

[TABLE]

By recurrence, the quantities $I_{k}$ for $2\leq k\leq n-1$ are all equal to zero if and only if, for every subset $I=\{i_{1},...,i_{k}\}\subset[n]$ of cardinality $k$ between $2$ and $n-1$, the variables $X_{i_{1}},...,X_{i_{k}}$ are independent. Suppose this is the case. In the above formula (30), we can replace all the intermediary higher entropies $H(X_{I};P)$ for $I$ between $2$ and $n-1$ by the corresponding sum of the individual entropies $H(X_{i_{1}})+...+H(X_{i_{k}})$. By symmetry each term $H(X_{i})$ appears the same number of time, with the same sign each time. The total sum of signs is obtained by replacing each $H(X_{i})$ by $1$; it is

[TABLE]

However, as a polynomial in $x$, we have

[TABLE]

thus

[TABLE]

therefore

[TABLE]

because $n\geq 2$.

Then we obtain

[TABLE]

Therefore, if the variables $X_{i};i=1,...,n$ are all independent, the quantity $I_{n}$ is equal to [math]. And conversely, if $I_{n}=0$, the variables $X_{i};i=1,...,n$ are all independent. ∎

Te Sun Han established that, for any subset $I_{0}$ of $[n]$ of cardinality $k_{0}\geq 2$, there exist probability laws such that all the $I_{k}(X_{I}),k\geq 2$ are zero with the exception of $I_{k_{0}}(X_{I_{0}})$ [21, 22]. Consequently, in the equations of the theorem 2.2, no one can be forgotten.

The unique equation (29) also characterizes the statistical independence, but its gradient with respect to $P$ is strongly degenerate along the variety of independent laws. As shown by Te Sun Han [21, 22], this is not the case for the $I_{k}$.

2.3 Information coordinates

The number of different functions $\eta_{I}$, resp. pure $I_{k}$, resp. pure $H_{k}$, is $2^{n}-1$ in the three cases. It is natural to ask if each of these families of functions of $P_{X}$ are analytically independent; we will prove here that this is true. The basis of the proof is the fact that each family gives finitely ambiguous coordinates in the case of binary variables, i.e. when all the numbers $N_{i},i=1,...,n$ are equal to $2$. Then we begin by considering $n$ binary variables with values [math] or $1$.

Let us look first at the cases $n=1$ and $n=2$. And consider only the family $H_{k}$, the other families being easily deduced by linear isomorphisms.

In the first case the only function to consider is the entropy

[TABLE]

if we put $x=p_{0}$. To get a probability $(p_{0},p_{1})$ we impose that $x$ belongs to $[0,1]$. As a function of $x$, $h$ is strictly concave, attaining all values between [math] and $1$, but it is not injective, due to the symmetry $x\mapsto 1-x$, which corresponds to the exchange of the values [math] and $1$.

For $n=2$, we have two variables $X_{1},X_{2}$ and three functions $H(X_{1};P)$, $H(X_{2};P)$, $H(X_{1},X_{2};P)$. These functions are all concave and real analytic in the interior of the simplex of dimension $3$.

Let us describe the probability law by four positive numbers $p_{00},p_{01},p_{10},p_{11}$ of sum $1$. The marginal laws for $X_{1}$ and $X_{2}$ are described respectively by the following coordinates:

[TABLE]

For the values of $H(X_{1};P)$ and $H(X_{2};P)$ we can take independently two arbitrary real numbers between [math] and $1$. Moreover, from the case $n=1$, if two laws $P$ and $P^{\prime}$ give the same values $H_{1}$ and $H_{2}$ of $H(X_{1};P)$ and $H(X_{2};P)$ respectively, we can reorder [math] and $1$ independently on each variable in such a manner that the images of $P$ and $P^{\prime}$ by $X_{1}$ and $X_{2}$ coincide, i.e. we can suppose that $p_{0}=p^{\prime}_{0}$ and $q_{0}=q^{\prime}_{0}$, which implies $p_{1}=p^{\prime}_{1}$ and $q_{1}=q^{\prime}_{1}$, due to the condition of sum $1$. It is easy to show that the third function $H(X_{1},X_{2};P)$ can take any value between the maximum of $H_{1}$, $H_{2}$ and the sum $H_{1}+H_{2}$.

Lemma 9**.**

There exist at most two probability laws that have the same marginal laws under $X_{1}$ and $X_{2}$ and the same value $H$ of $H(X_{1},X_{2})$; moreover, depending on the given values $H_{1},H_{2},H$ in the allowed range, both cases can happen in open sets of the simplex of dimension seven.

Proof.

When we fix the values of the marginals, all the coordinates $p_{ij}$ can be expressed linearly in one of them, for instance $x=p_{00}$:

[TABLE]

Note that $x$ belongs to the interval $I$ defined by the positivity of all the $p_{ij}$:

[TABLE]

The fundamental formula gives the two following equations:

[TABLE]

We define the functions $f_{1}(x)$ and $f_{2}(x)$ by the two above formulas respectively. As a function of $x$, each one is strictly concave, being a sum of strictly concave functions, thus it cannot take the same value for more than two values of $x$.

This proves the first sentence of the lemma; to prove the second one, it is sufficient to give examples for both situations.

Remark that the functions $f_{1},f_{2}$ have the same derivative:

[TABLE]

This results from the formula $h^{\prime}(u)=-\log_{2}(u/1-u)$ of the derivative of the entropy.

Then the maximum of $f_{1}$ or $f_{2}$ on $[0,1]$, is attained for $p_{01}p_{10}=p_{00}p_{11}$, that is when

[TABLE]

which we could have written without computation, because it corresponds to the independence of the variables $X_{1}$, $X_{2}$.

Then the possibility of two different laws $P,P^{\prime}$ in the lemma is equivalent to the condition that $p_{0}q_{0}$ belongs to the interior of $I$. This happens for instance for $1>p_{0}>q_{0}>q_{1}>p_{1}>0$, where $I=[q_{0}-p_{1},q_{0}]$, because in this case $p_{0}q_{0}<q_{0}$ and $p_{1}>p_{1}q_{0}$ i.e. $p_{0}q_{0}=q_{0}-p_{1}q_{0}>q_{0}-p_{1}$. In fact, to get $P\neq P^{\prime}$ with the same $H$, it is sufficient to take $x$ different from $p_{0}q_{0}$ but sufficiently close to it, and $H=f_{2}(x)+H_{2}$.

However, even in the above case, the values of $f_{1}$ (or $f_{2}$) at the extremities of $I$ do not coincide in general. Let us prove this fact. We have

[TABLE]

When $p_{1}=0$, the interval $I$ is reduced to the point $q_{0}$, and $F(0)=G(0)=0$. Now fix $q_{0},q_{1}$, and consider the derivatives of $F,G$ with respect to $p_{1}$ at every value $p_{1}>0$:

[TABLE]

Therefore $F^{\prime}(p_{1})<G^{\prime}(p_{1})$ if and only if $p_{0}-q_{0}<q_{0}-p_{1}$, i.e. $q_{0}>1/2$. Then, when $q_{0}>1/2$,for $p_{1}>0$ near [math], we have $F(p_{1})<G(p_{1})$.

Consequently, any value $f_{2}(x)$ that is a little larger than $F(p_{1})$ determines a unique value of $x$. It is in the vicinity of $q_{0}$. Which ends the proof of the lemma. ∎

From this lemma, we see that there exist open sets where $8$ or $4$ different laws give the same values of the three functions $H(X_{1}),H(X_{2}),H(X_{1},X_{2})$. In degenerate cases, we can have $4$, $2$ or $1$ laws giving the same three values.

Theorem 2.3**.**

*For $n$ binary variables $X_{1},...,X_{n}$, the functions $\eta_{I}$, resp. pure $I_{k}$, resp. pure $H_{k}$, characterize the probability law on $E_{X}$ up to a finite ambiguity.

Proof.

From the preceding section, it is sufficient to establish the theorem for the functions $H_{k}(X_{I})$, where $k$ goes from $1$ to $n$, and $I$ describes all the subsets of cardinality $k$ in $[n]$.

The proof is made by recurrence on $n$. We just have established the cases $n=1$ and $n=2$.

For $n>2$ we use the fundamental formula

[TABLE]

By the Marginal Theorem of H.G. Kellerer [28] (see also F. Matus [37]), knowing the $2^{n}-2$ non-trivial marginal laws of $P$, there is only one resting dimension, thus one of the coordinates $p_{i}$ only is free, that we denote $x$. Suppose that all the values of the $H_{k}$ are known, the hypothesis of recurrence tells that all the non-trivial marginal laws are known from the values of the entropy, up to a finite ambiguity. We fix a choice for these marginals. The above fundamental formula expresses $H(X_{1},...,X_{n})$ as a function $f(x)$ of $x$, which is a linear combination with positive coefficients of the entropy function $h$ applied to various affine expressions of $x$; therefore $f$ is a strictly concave function of one variable, then only two values at most are possible for $x$ when the value $f(x)$ is given. ∎

The group $\{\pm 1\}^{n}$ of order $2^{n}$ that exchanges in all possible manners the values of the binary variables $X_{i},i=1,...,X_{n}$ gives a part of the finite ambiguity. However, even for $n=2$, the ambiguity is not associated to the action of a finite group, contrarily to what was asserted in [4] section $1.4$. What replaces the elements of a group are partially defined operations of permutations that deserve to be better understood.

Theorem 2.4**.**

The functions $\eta_{I}$, resp. the pure $I_{k}(X_{I})$, resp. the pure $H_{k}(X_{I})$, have linearly independent gradients in an open dense set of the simplex $\Delta([n])$ of probabilities on $E_{X}$.

Proof.

Again, it is sufficient to treat the case of the higher pure entropies.

We write $N=N_{1}...N_{n}$. The elements of the simplex $\Delta(N)$ are described by vectors $(p_{1},...,p_{N})$ of real numbers that are positive or zero, with a sum equal to $1$. The expressions $H_{k}(X_{J})$ are real analytic functions in the interior of this simplex. The number of these functions is $2^{n}-1$. The dimension $N-1$ of the simplex is larger (and equal only for the fully binary case), then to establish the result, we have to find a minor of size $2^{n}-1$ of the Jacobian matrix of the partial derivatives of the entropy functions with respect to the variables $p_{i},i=1,...,N-1$ that is non identically zero. For any index $j$ between $1$ and $n$ choose two different values of the set $E_{j}$. Then apply the theorem $2$. ∎

Remark**.**

This proves the fact mentioned in $1.4$ of $[3]$.

Te Sun Han established that the quantities $I_{k}(X_{I})$ for $k\geq 2$ are functionally independent [21, 22].

Remark**.**

The formulas of $H_{k}(X_{I})$, then of $I_{k}(X_{i})$ and $\eta_{I}$, extend analytically to the open cone $\Gamma([n])$ of vectors with positive coordinates. On this cone we pose

[TABLE]

*This is the natural function to consider to account for the empty subset of $[n]$.

Be careful that the functions $K_{k}$ for $k>0$ are no more positive in the cone $\Gamma([n])$, because the function $-x\ln x$ becomes negative for $x>0$. In fact we have, for $\lambda\in]0,\infty[$, and $P=(p_{1},...,p_{n})\in\Gamma([n])$,*

[TABLE]

The above theorems extend to the prolonged functions to the cone, by taking into account $H_{0}$.

Notice further properties of information quantities:

For $I_{k}$, due to the constraints on $I_{2}$ and $I_{3}$, see Matsuda [36], we have for any pair of variables

[TABLE]

and any triple $X_{1},X_{2},X_{3}$:

[TABLE]

It could be that interesting inequalities also exist for $k\geq 4$, but it seems that they are unknown.

Contrarily to $H_{k}$, the behavior of the function $I_{k}$ is not the same for $k$ even and $k$ odd. In particular, as functions of the probability $P_{X}$, the odd functions $I_{2m+1}$, for instance $I_{1}=H_{1}=H$, or $I_{3}$ (ordinary synergy), have properties of the type of pseudo-concave functions (in the sense of [4]), and the even functions $I_{2m}$, like $I_{2}$ (usual mutual information) have properties of the type of convex functions (see [4] for a more precise statement). Note that this accords well with the fact that the total entropy $H(X)$, which is concave, is the alternate sum of the $I_{k}(X_{I})$ over the subsets $I$ of $[n]$, with the sign $(-1)^{k-1}$ (cf. appendix A).

Another difference is that each odd function $I_{2m+1}$ is an information co-cycle, in fact a co-boundary if $m\geq 1$ (in the information co-homology defined in [4]), but each odd function $I_{2m+1}$ is a simplicial co-boundary in the ordinary sense, and not an information co-cycle.

Remark**.**

From the quantitative point of view, we have also considered and quantified on data the pseudo-concave function $(-1)^{k-1}I_{k}$ (in the sense of [4]) as a measure of available information in the total system and considered total variation along paths. Although such functions are sounding and appealing, we have chosen to illustrate here only the results using the function $I_{k}$ as they respect and generalize the usual multivariate statistical correlation structures of the data and provide meaningful data interpretation of positivity and negativity, as it will become obvious in the following application to data. However, what really matters is the full landscape of information sequences, showing that information is not well described by a unique number, but rather by a collection of numbers indexed by collections of joint variables.

2.4 Application to Gene expression data - detection of cell types and gene modules

The developments and tests of the estimation of simplicial information topology on data is made on a genetic expression dataset of two cell types obtained as described in the section material and methods 4.1. The result of this quantification of gene expression is represented in ”Heat maps” and allows two kinds of analysis:

• •

The analysis with genes as variables: in this case the ”Heat maps” correspond to $(m,n)$ matrices $D$ (presented in the section 4.2) together with the labels (population A or population B) of the cells. The data analysis consists in the detection of gene modules.

• •

The analysis with cells (neurons) as variables: in this case the ”Heat maps” correspond to the transposed matrices $D^{T}$ (presented in the Figure 3) together with the labels (population A or population B) of the cells. The data analysis consists in the detection of cell types.

2.5 Mutual-information negativity, clusters and links

As information negativity posed problems of interpretation as recalled in introduction and conclusion, we now illustrate what the negative and positive information values quantify on data with theoretical (taking the case of the binary variable case previously exposed) and empirical N-ary case examples of gene expression. Let us consider three ordinary biased coins $X_{1},X_{2},X_{3}$, we will denote by [math] and $1$ their individual states and by $a,b,c,...$ the probabilities of their possible configurations three by three; more precisely:

[TABLE]

We have

[TABLE]

The following identity is easily deduced from the definition of $I_{3}$ (cf. 18):

[TABLE]

Of course the identities obtained by changing the indices are also true. This identity interprets the information shared by three variables as a measure of the lack of information in conditioning. We notice a kind of intrication of $I_{2}$: conditioning can increase the information, which interprets rightly the negativity of $I_{3}$. Another useful interpretation of $I_{3}$ is given by

[TABLE]

In this case negativity is interpreted as a synergy, i.e. the fact that two variables can give more information on a third variable than the sum of the two separate information.

Several inequalities are easy consequences of the above formulas and of the positivity of mutual information of two variables (conditional or not), as shown in [36].

[TABLE]

[TABLE]

, and the analogs that are obtained by permuting the indices.

Let us remark that this immediately implies the following assertions:

1. when two variables are independent the information of the three is negative or zero;

2. when two variables are conditionally independent with respect to the third, the information of the three is positive or zero.

By using the positivity of the entropy (conditional or not), we also have:

[TABLE]

[TABLE]

We deduce from here

[TABLE]

[TABLE]

In the particular case of three binary variables, this gives

[TABLE]

Proposition 10**.**

The absolute maximum of $I_{3}$, equal to $1$, is attained only in the four cases of three identical or opposite unbiased variables. That is $H(X_{1})=H(X_{2})=H(X_{3})=1$, and $X_{1}=X_{2}$ or $X_{1}=1-X_{2}$, and $X_{1}=X_{3}$ or $X_{1}=1-X_{3}$, that is $a=h=1/2$ or $b=g=1/2$ or $c=f=1/2$ or $d=e=1/2$ and in each case all the other variables are equal to [math] (cf. Figure 1 a-c).

Proof.

First it is evident that the example gives $I_{3}=1$. Second, consider three variables such that $I(X_{1};X_{2};X_{3})=1$. We must have $H(X_{1})=H(X_{2})=H(X_{3})=1$, and also $I(X_{i};X_{j})=1$ for any pair $(i,j)$, thus $H(X_{i},X_{j})=1$, $H(X_{i}|X_{j})=0$, and the variable $X_{i}$ is a deterministic function of the variable $X_{j}$, which gives $X_{i}=X_{j}$ or $X_{i}=1-X_{j}$. ∎

Proposition 11**.**

The absolute minimum of $I_{3}$, equal to $-1$, is attained only in the two cases of three two by two independent unbiased variables satisfying $a=1/4,b=0,c=1/4,d=0,e=1/4,f=0,g=1/4,h=0$, or $a=0,b=1/4,c=0,d=1/4,e=0,f=1/4,g=0,h=1/4$. These cases correspond to the two borromean links, the right one and the left one (cf. Figure 1 d-f).

Proof.

First it is easy to verify that the examples give $I_{3}=-1$. Second consider three variables such that $I(X_{1};X_{2};X_{3})=-1$. The inequality (62) implies $H(X_{1})=H(X_{2})=H(X_{3})=1$, and the inequality (60) shows that $H(X_{i}|X_{j})=1$ for every pair of different indices, so $H(X_{1},X_{2})=H(X_{2},X_{3})=H(X_{3},X_{1})=2$, and the three variables are two by two independent. Consequently the total entropy $H_{3}$ of $(X_{1},X_{2},X_{3})$, given by $I_{3}$ minus the sum of individual entropies plus the sum of two by two entropies is equal to $2$. Thus

[TABLE]

But we also have

[TABLE]

that is

[TABLE]

Now we subtract (66) from (64), we obtain

[TABLE]

However each of the four quantities $-4a\lg 4a-4b\lg 4b,-4c\lg 4c-4d\lg 4d,\\ -4e\lg 4e-4f\lg 4f,-4g\lg 4g-4h\lg 4h$ is $\geq 0$ because each of the four sums $4a+4b,4c+4d,4e+4f,4g+4h$ is equal to $1$, so each of these quantities is equal to zero, which happens only if $ab=cd=ef=gh=0$. But we can repeat the argument with any permutation of the three variables $X_{1},X_{2},X_{3}$. We obtain nothing new from the transposition of $X_{1}$ and $X_{3}$. From the transposition of $X_{1}$ and $X_{3}$ we obtain $ae=bf=cg=dh=0$. From the transposition of $X_{2}$ and $X_{3}$, we obtain $ac=bd=eg=fh=0$. So from the cyclic permutation $(1,2,3)$ (resp. $(1,3,2)$, we get $ae=bf=cg=dh=0$ (resp. $ac=bd=eg=fh=0$).) If $a=0$ this gives necessarily $b,e,c$ nonzero, thus $d=f=g=0$, and $h\neq 0$, and if $a\neq 0$ this gives $b=e=c=0$, thus $d,f,g$ nonzero and $h=0$. ∎

Figure 1 illustrates the probability configurations giving rise to the maxima and minima of $I_{3}$ for 3 binary variables.

In the much more complex case of gene expressions, the statistical analysis shown in [58] exhibited also a combination of positivity and negativity of the information quantities $I_{k};k\geq 3$. In this analysis, the minimal negative information configurations provide a clear example of purely emergent and collective interactions analog to Borromean links in topology, since it cannot be detected from any pairwise investigation or 2-dimensional observations. In these Borromean links the variables are pairwise independent but dependent at 3. In general $I_{k}$ negativity detects such effects of their projection on lower dimensions, this illustrates the main difficulty when going from dimension $2$ to $3$ in information theory. The example given in Figure 1 provides a simple example of this dimensional effect in the data space: the alternated clustering of the data corresponding to $I_{3}$ negativity cannot be detected by the projections onto whichever subspace of pair of variables, since the variables are pairwise independent. For N-ary variables the negativity becomes much more complicated, with more degeneracy of the minima and maxima of $I_{k}$.

In order to illustrate the theoretical examples of Figure 1 on real data, considering the data set of gene expression (matrix $D$), we plotted some quadruplets of genes sharing some of the highest (positive) and lowest (negative) $I_{4}$ values in the data space of the variables (Figure 2). Figure 2 shows that in the data space, $I_{k}$ negativity identifies the clustering of the data points, or in other words, the modules (k-tuples) for which the data points are segregated into condensate clusters. As expected theoretically, $I_{k}$ positivity identifies co-variations of the variables, even in cases of non-linear relations, as shown by Reshef and colleagues [46] in the pairwise case. It can be easily shown in the pairwise case that $I_{k}$ positivity generalizes the usual correlation coefficient to non-linear relations. As a result, the interpretation of the negativity of $I_{k}$ is that it provides a signature and quantification of the variables that segregate or differentiate the measured population.

2.6 Cell type detection - comparison with previous MaxEnt studies

Example of cell type recognition with a low sample size $m=41$, dimension $n=20$, and graining $N=9$.

As introduced in previous section 2.5, the k-tuples presenting the highest and lowest information ($I_{k}$) values are the most biologically relevant modules and identify the variables that are the most dependent or synergistic (respectively ”entangled”). We call information landscape the representation of the estimation of all $I_{k}$ values for the whole simplicial lattice of k-subfaces of the $n$-simplex of variables ranked by their $I_{k}$ values in ordinate. In general the null hypothesis against whom are tested the data is the maximal uniformity and independence of the variables $X_{i},i=1,...,n$. Below the undersampling dimension $k_{u}$ presented in methods 4.5, this predicts the following standard sequence for any permutation of the variables $X_{i_{1}},...,X_{i_{n}}$;

[TABLE]

that is linearity (with $N_{1}=...=N_{n}=r$).

What we observed in the case where independence is confirmed, for instance with the chosen genes of the population B (NDA neurons) in [58], is linearity up to the maximal significant $k$, then stationarity. But where independency is violated, for example with the chosen genes of the population A (DA neurons) in [58], some permutations of $X_{1},...,X_{n}$ give sequences showing strong departures from the linear prediction.

This departure and the rest of the structure can also be observed on the sequence $I_{k}$ as shown in Figure 3 and 4, which present the case where cells are considered as variables. In the trivial case, i.e. uniformity and independence, for any permutation, we have

[TABLE]

As detailed in materials and methods 4.3, we further compute the longest information paths (starting at 0 and that go from vertex to vertex following the edges of the simplicial lattice) with maximal or minimal slope (with minimal or maximal conditional mutual-information) that end at the first minimum, a conditional-independence criterion (a change of sign of conditional mutual-information). Such paths select the biologically relevant variables that progressively add more and more dependences. The paths $I_{k}(\sigma)$ that stay strictly positive for a long time are especially interesting, being interpreted as the succession of variables $X_{\sigma_{1}},...,X_{\sigma_{k}}$ that share the strongest dependence. However, the paths $I_{k}(\sigma)$ that become negative for $k\geq 3$ through $I_{2}\approx 0$ are also interesting, because they exhibit a kind of frustration in the sense of Matsuda [36] or synergy in the sense of Brenner [12].

The information landscape and path analysis corresponding to the analysis with cells as variables are illustrated in Figure 3. It comes to consider the cells as a realization of gene expression rather than the converse, cf. [14]. In this case, the data analysis task is to recover blindly the pre-established labels of cell types (population A and population B) from the topological data analysis, an unsupervised learning task. The heat-map transpose matrix of $n=20$ cells with $m=41$ genes is represented in Figure 3a. We took $n=20$ neurons among the $148$ within which $10$ were pre-identified as population A neurons (in green) and $10$ were pre-identified as population B neurons (in dark red), and ran the analysis on the $41$ gene expression with a graining of $N=9$ values (cf. section 4.1). The dimension above which the estimation of information becomes too biased due to the finite sample size is given by the undersampling dimension $k_{u}=11$ (p value 0.05, cf. section 4.5). The landscapes turn out to be very different from the extremal (totally disordered and totally ordered) homogeneous (identically distributed) theoretical cases. The $I_{k}$ landscape shown in Figure 3c exhibits two clearly separated components. The scaffold below represents the tuple corresponding to the maximum of $I_{10}$: it corresponds exactly to the $10$ neurons pre-identified as being population A neurons.

The maximum (in red) and minimum (in blue) $I_{k}$ information paths identified by the algorithm are represented in Figure 3d. The scaffold below represents the two tuples corresponding to the two longest maximum paths in each component: the longest (noted Max $IP_{11}$ in green) $IP_{11}$ contains the 10 neurons pre-identified as population A and 1 ”error” neuron pre-identified as population B. We restricted the longest maximum path to the undersampling dimension $k_{u}=11$, but this path reached $k=14$ with erroneous classifications. The second longest maximum path (noted Max $IP_{11}$ in red) $IP_{11}$ contains the 10 neurons pre-identified as population B and 1 neuron pre-identified as population A that is hence erroneously classified by the algorithm. Altogether the information landscape shows that population A neurons constitute a quite homogenous population, whereas the population B neurons correspond to a more heterogeneous population of cells, a fact that was already known and reported in the biological studies of these populations. The histograms of the distributions of $I_{k}$ for $k=1,..,12$, shown in Figure 3e are clearly bimodal and the insets provide a magnification on the population A component. As detailed in the section materials and methods 4.6, we developed a test based on the random shuffles of the data points that leave the marginal distributions unchanged, as proposed by [43]. It estimates if a given $I_{k}$ significantly differs from a randomly generated $I_{k}$, a test of the specificity of the k-dependence. The shuffled distributions and the significance value for $p=0.1$ are depicted by the black lines and the doted lines, as in Figure 8. As illustrated in the histograms of Figure 3e and in [57], these results show that higher dependences can be important but they do not mean that pairwise or marginal informations are not: the consideration of higher dependences can only improve the efficiency of the detection obtained from pairwise or marginal considerations.

As illustrated in Figure 4 and expected from relative entropy positivity, the total correlation $G_{k}$ (see appendix A on bayes free-energy) is monotonically increasing with the order $k$, and quite linearly in appearance ($G_{k}\approx 2k$ asymptotically). The panel d quantifies this departure from linearity. However the $G_{k}$ landscape fails to distinguih as clearly as $I_{k}$ landscape does, the population A.

3 Discussion

During the last decades, there have been important efforts in trying to evaluate the pairwise and higher order interactions in neuronal and biological measurements, notably to extract the undergoing collective dynamics. Applying the Maximum of Entropy principle on Ising spin models to neural data [48, 59], a first series of studies concluded that pairwise interactions are mostly sufficient to capture the global collective dynamics, leading to the ”pairwise sufficiency” paradigm (see Merchan and Nemenman for presentation [39]). However, as shown by the Ising model itself, near a second order phase transition, elementary pairwise interactions are compatible with non-trivial higher-order dependences, and very large correlations at long distances. From the mathematical and physical point of view, this fact is nicely encoded in the normalization factor of the Boltzmann probability, namely the Partition Function $Z(\beta)$. As shown by the Ising model, the probability law can be factorized (up to the normalization number $Z$) on the edges and vertices of a graph, but the statistical clusters can have unbounded volumes. Moreover, subsequent studies notably of Tkac̆ik et al. [60] (see also [25]) have shown that for sufficiently large populations of registered neurons, the pairwise models are insufficient to explain the data as proposed in [2, 3] for example. Thus the dimension of the interactions to be taken into account for the models must be larger than two.

Note that most interactions in Biology are nowadays described in terms of networks, such that the concepts of protein networks, genetic networks or neural networks became familiar. However from the physical as well as the biological point of view, none of these systems are really $1$-dimensional graphs, and it is now clear for most researchers in the domain that higher order structures are needed for describing collective dynamics, cf. for instance [69], and [45]. Our Figure 5 clearly shows this point.

William Bialek and his collaborators have well explained the interest of a systematic study of joint entropies and general multi-modal mutual information quantities as an efficient way for understanding neuronal activities, networks of neurons, and gene expression [12, 49, 54]). They also developed approximate computational methods for estimating the information quantities. Mutual information analysis was applied for linking adaptation to the preservation of the information flow [10, 33]. Closely related to the present study, Margolin, Wang, Califano and Nemenman have investigated multivariate dependences of higher order [35] with MaxEnt methods, by using the total-correlation $G_{k}$ (cf. equation 28) in function of the integer $k\geq 2$. The apparent advantage is the positivity of the $G_{k}$.

In this respect, the originality of our method relies first on the systematic consideration of the entropy paths and the information paths that can be associated to all possible permutations of the basic variables, in arbitrary dimension, and the extraction of exceptional paths from them, in order to define the overall form of the distribution of information among the set of variables. Secondly, we used and proved the relevance of peculiar functions, multivariate mutual informations, where the previously cited works focused on total correlations, which fail to uncover the data structure as exemplified in Figure 4 or only explored pairwise or $I_{3}$. We named these tools the information landscapes of the data. This new perspective and mathematical justification of these functions has its origin in the local (topos) homological theory introduced in [4] developped and extended in several ways by Vigneaux [64]. In the present article, we also proved new theoretical results along this line, about the concrete structure of higher-order information functions. Moreover the method was successfully applied to a concrete problem of gene expression in [58].

Since their introduction, the possible negativity of the $I_{k}$ functions for $k\geq 3$ has posed serious problems of interpretation, and it was the main argument for many theoretical studies to discard such a family of functions for measuring information dependences and statistical interactions. Notably, it motivated the proposition of non-negative decomposition by Williams and Beer [68] and of ”unique information” by Bertschinger and colleagues [42, 7], or Griffith and Koch [19]. These partial decompositions of information are the subject of several recent investigations notably with applications to the development of neural network [67] and neuromodulation [27]. However, Rauh and colleagues showed that no non-negative decomposition can be generalized to multivariate cases for degrees higher than 3 [44] (th.2).

The present paper and [58] show that, to the contrary, the possible negativity is an advantage. The interest of this negativity was already illustrated in [12, 49, 36, 30], but we have further developed this topic with the study of complete $I_{k}$-landscapes, providing some new insights with respect to their meaning in terms of data point clusters.

The precise contribution of higher-order is indeed directly quantified by the $I_{k}$ values in the landscapes and paths. Figure 5 further illustrates the gain and the importance of considering higher statistical interactions, using the previous example of cells pre-identified as 10 population A and 10 population B cells ($n=20$, $m=47$, $N=9$). The plots are the finite and discrete analogs of Gibbs’s original representation of entropy vs. energy [17]. Whereas pairwise interactions ($k=2$) cannot (or very hardly) distinguish the population A and population B cell types, the maximum of $I_{10}$ unambiguously identifies the population A.

As illustrated in Figure 3, the present analysis shows that in the expression of 41 genes of interest of population A neurons, the higher-order statistical interactions are non-negligible and have a simple functional meaning of a collective module, a cell type. We believe such conclusion to be generic in biology. More precisely, we believe that even if related to physics, biological structures have higher-order statistical interactions defined by higher-order information and that these interactions provide the signature of their memory engramming. In fact ”information is physical” as stated by Bennett following Landauer [32]), in the sense of memory capacity and necessity of forgetting. The quantification of the information storage applied here to genes can be considered as a generic epigenetic memory characterization, resulting of a developmental-learning process. The consideration of higher dimensional statistical dependences increases combinatorially the number of possible information modules engrammed by the system. It hence provides an appreciable capacity reservoir for information storage and for differentiation, for diversity.

The critical points of the Ising model in dimension 2 and 3 show the difficulty to relate factorization (up to $Z(\beta)$), which describes the manner energy interactions localize, with the dependences structure, or in other words the manner information distributes itself, i.e. the form of information. Only few theoretical results relate the two notions. However, on the basis of several recent studies that we mentioned, particularly the studies of adaptive functions, and comforted by the analysis presented in this article, we can suggest that for biological systems, during development or evolution, the distribution of the information flow, as described in particular by higher order information quantities, participates in the generators of the dynamics, on the side of energy quantities coming from Physics.

4 Material and Methods

4.1 The dataset: quantified genetic expression in two cell types

The quantification of genetic expression was performed using microfluidic qPCR technique on single dopaminergic (DA) and non-dopaminergic (NDA) neurons isolated from two midbrain structures, the Substantia Nigra pars compacta (SNc) and the neighboring Ventral Tegmental Area (VTA), extracted from adult TH-GFP mice (transgenic mice expressing the Green Fluorescent Protein under the control of the Tyrosine Hydroxylase promoter). The precise protocols of extraction, quantification, and identification are detailed in [57, 58]. This technique allowed us to quantify in a single cell the levels of expression of 41 genes chosen for their implication in neuronal activity and identity of dopaminergic (DA) neurons. The SNc DA neurons were identified based on GFP fluorescence (TH expression). This identification was further confirmed based on the expression levels of Th and Slc6a3 genes, which are established markers of DA metabolism. The quantification of the expression of the 41 genes ($n=41$) was achieved in 111 neurons ($m=111$) identified as DA and in 37 neurons ($m=37$) identified as nDA. In this article, for readability purpose, we replaced the names of the genes by gene numbers and the cell type DA by population A, and the cell type nDA by population B. The dataset is available in supplementary material [57, 58].

4.2 Probability estimation

The presentation of the probability estimation procedure is achieved on matrices $D$ (genes as variables), and it is the same in the case of the analysis of the matrices $D^{T}$ (cells as variables). It is illustrated in Figure 6 for the simple case of 2 random variables taken from the dataset of gene expression presented in section 4.1, namely the expression of two genes Gene5 and Gene21 in $m=111$ population A cells. Our probability estimation corresponds to a step of the integral estimation procedure of Riemann.

We write the heatmap as a $(m,n)$ matrix $D$ and its real coefficients $x_{ij}\in\mathbb{R},~{}i\in\{1..m\},~{}j\in\{1...n\}$: the columns of $D$ span the $m$ repetitions-trials (here the $m$ neurons) and the rows of $D$ spans the $n$ variables (here the $n$ genes). We also note, for each variable $X_{j}$, the minimum and maximum values measured as $\min x_{j}=\min_{1\leq i\leq m}x_{ij}$ and $\max x_{j}=\max_{1\leq i\leq m}x_{ij}$.

We consider the space in the intervals $[\min x_{j},\max x_{j}]$ for each variable $X_{j}$ and divide it into $N_{1}.N_{2}...N_{n}$ boxes, on which it is possible to estimate the atomic probabilities by elementary counting. We note each $n$-dimensional box by an $n$-tuple of integers $\{a_{1},...,a_{n}\}$ where $\forall i\in\{1,...,n\},~{}a_{i}\in\{1,...,N_{i}\}$, and writing the min and the max of a box on each variable $X_{j}$ (the jth co-ordinate of the vertex of the box) as $\operatorname{bmin}_{j}=\min x_{j}+\frac{(a_{j}-1)(\max x_{j}-\min x_{j})}{N_{j}}$ and $\operatorname{bmax}_{j}=\min x_{j}+\frac{(a_{j})(\max x_{j}-\min x_{j})}{N_{j}}$, then the atomic probabilities can be defined using Dirac function $\delta$ as:

[TABLE]

For two variables, using the definition of conditioning $P_{X}(Y)=\frac{P(X.Y)}{P(X)}$ and in the general case using the theorem of total probability [31] ($P(X)=\sum_{i=0}^{N}P(A_{i}.X)=\sum_{i=0}^{N}P(A_{i}).P_{A_{i}}(X)$), we can marginalize, or geometrically project on lower dimensions, to obtain all the probabilities corresponding to subsets of variables, as illustrated in Figure 6. For example, with short notation, the probability associated to the marginal variable $X_{i}$ being in the interval $[\operatorname{bmin}_{i},\operatorname{bmax}_{i}]$ is obtained by direct summation:

[TABLE]

In the example of Figure 6, the probability of the level of Th being in the 4th box is:

[TABLE]

In geometrical terms, the set of total probability laws is an $N=N_{1}.N_{2}...N_{n}-1$ dimensional simplex $\Delta_{N_{1}.N_{2}...N_{n}-1}$ (the $-1$ accounts for the normalization equation $\sum P_{i}=1$, which embeds the simplex in an affine space). In the example of Figure 6, we have an 80-dimensional probability simplex $\Delta_{80}$, the set of sub-simplicies over the k-faces of the simplex $\Delta_{n}$, for every $k$ between [math] and $n$, represents the boolean algebra of the joint-probabilities, which is equivalent in the finite case to their sigma-algebra. In our analysis, we have chosen $N_{1}=N_{2}=...=N_{n}=9$ and this choice is justified in section 4.6.1 using Reshef and colleagues criterion [46] and undersampling constraints.

In summary, our probability estimation and data analysis depend on $n$ (the number of random variables), on $m$ (the number of observations), and on $N_{1},...,N_{i}$ (the graining). The merit of this method is its simplicity (few assumptions, no priors on the distributions) and low computational cost. There exist different methods that can significantly improve this basic probability estimation, but we leave this for future investigation. The graining given by the numbers $N=N_{1}.N_{2}...N_{n}$ and the sample size $m$ are important parameters of the analysis explored in this section.

4.3 Computation of k-Entropy, k-Information landscapes and paths

The computational exploration of the simplicial sublattice has a complexity in $\mathcal{O}(2^{n})$ ($2^{n}=\sum_{k=1}^{n}\binom{n}{k}$). In this simplicial setting we can exhaustively estimate information functions on the simplicial information structure, that is joint-entropy $H_{k}$ and mutual-informations $I_{k}$ at all dimensions $k\leq n$ and for every $k$-tuple, with a standard commercial personal computer (a laptop with processor Intel Core i7-4910MQ CPU @ 2.90GHz $\times$ 8, even though the program currently uses only one CPU) up to $k=n=21$ in a reasonable time ($\approx 3$ hours). Using the expression of joint-entropy (equation 11) and the probability obtained using equation 70 and marginalization, it is possible to compute the joint-entropy and marginal entropy of all the variables. The alternated expression of n-mutual information given by equation 12 then allows a direct evaluation of all these quantities. The definitions, formulas and theorems are sufficient to obtain the algorithm. We moreover provide the Information Topology program INFOTOPO-V1.2 under opensource licence on github depository at https://github.com/pierrebaudot/INFOTOPO. Information Topology is a program written in Python (compatible with Python 3.4.x), with a graphic interface built using TKinter [53], plots drawn using Matplotlib [26], calculations made using NumPy [62], and scaffold representations drawn using NetworkX [20]. It computes all the results on information presented in the current study, including the information paths, statistical tests of $I_{k}$ values described in the next sections and the finite entropy rate $\frac{H_{k}}{k}$. The input is an excel table containing the data values, e.g. the matrix $D$ with the first row and column containing the labels. Here, we limited our analysis to $n=21$ genes of specific biological interest.

4.4 Estimation of the undersampling dimension: statistical result

The information data analysis presented here depends on the two parameters $N$ and $m$. The finite size of the sample $m$ is known to impose an important bias in the estimation of information quantities: in high-dimensional data analysis, it is quoted as the Hugues phenomenon [24] and in entropy estimation it has been called the sampling problem since the seminal work of Strong and colleagues [55, 41, 39]. For the method we suggested, it is important to notice that the size $m$ of the population $Z$ is in general much smaller than the dimension of the probabilty simplex $N=N_{1}...N_{n}-1$. For instance, in the mentioned study of genes as variables [58], we had $m=111$ for $DA$ neurons (resp. $m^{\prime}=37$ for $NDA$ neurons) as respective number of neurons, but $N=9^{21}-1$, because we could only achieve the computation for the 21 most relevant genes. In the example considering cells as variables presented here in Figure 3, the situation is even worse, with a sample size of $m=41$ genes and a dimension of $N=9^{20}-1$ as only 20 cells were considered. Thus the pure entropies $H_{k},k=1,...,n$ must satisfy the following inequality:

[TABLE]

where equality is an extreme signature of undersampling. However, suppose that all the numbers $N_{i},i=1,...,n$ are equal to $r\geq 2$, the maximum value of $H_{k}$ is equal to $k\log_{2}r$, for instance $2k.\log_{2}(3)$ in the example.

Lemma 12**.**

Take the uniform probability on the simplex $\Delta([n])$ with affine coordinates, and take $\epsilon$ such that $0<\epsilon\leq 1/e\approx 0,367$; then the probability that $H_{k}(X_{J})$ is greater than $\epsilon k\log_{2}r$ is larger than $1-\epsilon$.

Proof.

Concerning $H_{k}$, the simplex $\Delta([n])$ is replaced by $\Delta([k])$; then consider the set $\Delta_{\epsilon}$ of probabilities such that $p_{j}\geq\epsilon r^{-k}$ for any coordinate $j$ between $1$ and $r^{k}$, this set is the complement of the union of the sets $X_{j}(\varepsilon),i=1,...,r^{k}$ where $p_{j}<\epsilon r^{-k}$. From the properties of volumes in affine geometry, the measure of each set $X_{j}(\varepsilon)$ is less than $\epsilon r^{-k}$, thus the probability of $\Delta_{\epsilon}$ is larger than $1-\epsilon$. And for any index $j$ the monotony of $-x\ln x$ between [math] and $1/e$ implies

[TABLE]

then by summation over all the indices we obtain the result. ∎

By example, for $r=9$, and $\epsilon=1/e$, this gives that $H_{k}\geq 2k\log_{2}(3)/e$ is two times more probable than the opposite.

Consequently, in the above experiment, the quantities $H_{k}$, then $I_{k}$ are not significant, except if they appear to be significantly smaller than $\log_{2}m$.

In counterpart, as soon as the measured $H_{k}$ is inferior than the predicted one for $m$ values, this is significant. Note that the lemma 12, with $n$ replaced by $m$, gives estimations for the entropies of raw data. In the next section, we propose a computational method to estimate the dimension $k_{u}$ above which information estimation ceases to be significant.

4.5 Estimation of the undersampling dimension: Computational result

Following the original presentation of the sampling problem by Strong and colleagues [55], the extreme cases of sampling are given by:

• •

When $N_{1}=N_{2}=...=N_{n}=1$, there is a single box $\Omega$ and $P(\Omega)=m/m=1$ and we have $H_{k}=I_{k}=0,\forall k\in{0,...,n}$. The case where $m=1$ is identical. This fixes the lower bound of our analysis in order not to be trivial; we need $m\geq 2$ and $N_{1}=N_{2}=...=N_{n}\geq 2$.

• •

When $N_{1}.N_{2}...N_{n}$ are such that only one data point falls into a box, $m$ of the values of atomic probabilities are $1/m$ and $N_{1}.N_{2}...N_{n}-m$ are null as a consequence of equation 71, and hence we have $H_{n}=\log_{2}m$.

Whenever this happens for a given k-tuple, all the $HP_{k}$ paths passing by this k-tuple will stay on the same information values since conditional entropy is non-negative: we have $H_{k}=H_{k+1}$ or equivalently $(X_{1},...,X_{k})H(X_{k+1})=0$, and all $k+l$-tuples are deterministic (a function of) with respect to the k-tuple. This is typically the case illustrated in Figure 3: adding a new variable to an undersampled k-tuple is equivalent to adding the deterministic variable ”0” since the probability remains unchanged ($1/m$).

Considering the analysis of cells as variables (matrix $D^{T}$), the signature of this undersampling is the saturation at $H_{k}=\log_{2}41$ observed in the $H_{k}$ landscape in Figure 3b, starting at $k=5$ for some 5-tuples of neurons. Considering the analysis of genes as variables (matrix $D$ [58]), the mean entropy computed also shows this saturation at $H_{k}=\log_{2}111$ for population A neurons and $H_{k}=\log_{2}37$ for population B neurons. We propose to define a dimension $k_{u}$ as the dimension for which the probability $p_{u}$ of having the $H_{k}$ at the biased value of $H_{k}=\log_{2}m$ is above $5$ percent ($p_{u}=0.05$). As shown for the analysis of cells as variables in Figure 7, this basic estimation gives here $k_{u}=6$ for population A neurons and $k_{u}=4$ for population B neurons. The information structures identified by our methods beyond these values can be considered as unlikely to have a biological or physical meaning and shall not be interpreted. Since undersampling mainly affects the distribution of $I_{k}$ values close to 0 value, the maxima and minima of $I_{k}$ and the maximal and minimal information paths below $k_{u}$ are the least affected by the sampling problem and the low sample size. This will be illustrated in the next sections.

4.6 k-dependence test

Pethel and Hahs [43] have constructed an exact test of 2-dependence for any pair of variables, not necessarily binary or iid. Indeed, the iid condition usually assumed for the $\chi^{2}$ test does not seem relevant for biological observations and the examples given here and in [57, 58] with genetic expression support such a general statement. It allows to test the significance of the estimated $I_{2}$ values given a finite sample size $m$, the null hypothesis being that $I_{2}=0$ (2-independence according to Pethel and Hahs). We follow here their presentation of the problem, and provide an extension of their test to arbitrary $k$ (higher dimensions), with the null hypothesis being the k-independence $I_{k}=0$. Even in the lowest dimensions, and below the undersampling bound, the values of $I_{k}$ estimated from a finite sample size $m$ are considered as biased [43]. If one considers an infinite sample ($m\rightarrow\infty$) of n independent variables, we then have for all $k\geq 2$ $I_{k}=0$. However, if we randomly shuffle the values such that the marginal distributions for each variable $X_{i}$ are preserved, the estimated $I_{k}$ can be very different from 0, with distributions of $I_{k}$ values not centered on 0. Figure 8 illustrates an example of such bias with $m=111$ for the analysis with genes as variables.

Reproducing the method of Pethel and Hahs [43], we designed a shuffling procedure of the $n$ variables, which consists in randomly permuting the measured values (co-ordinates) of each variable one by one in the matrix $D$ or $D^{T}$ (geometrically, a ”random” permutation of the co-ordinates of each data point, point by point). Such a shuffle leaves marginal probabilities invariant. Figure 8 gives an example of the joint and marginal distributions before and after shuffle for two genes. Extending the 2-test of [43] to $k\geq 2$, the $I_{k}$ values obtained after shuffling provide the distribution of the null hypothesis, k-independence ($I_{k}=0$) according to [43]. The task is hence to compute many shuffles, 10.000 in [43], in order to obtain these ”null” distributions. The exact procedure of Pethel and Hahs [43] would require to obtain such ”null” distribution for all the $2^{n}$ tuples, which would require a number of shuffled trials impossible to obtain computationally. We hence propose a global test that consists in computing 17 different shuffles of the 21 genes, giving ”null” distribution of shuffled $I_{k}$ values composed of $21\times\binom{n}{k}$. For example, the test of $2$-dependence and $3$-dependence will be against a null distribution with $21*210=3750$ $I_{2}$ values and $21*1330=22610$ $I_{3}$ values respectively. We fix a p value above which we reject the null hypothesis (a significance level, fixed at $p=0.05$ in [43]), allowing to determine the statistical significance thresholds as information values for which the integral of the null distribution reaches the significance level $p=0.05$. This holds for $k=2$, as described in [43], but since for $k\geq 2$ $I_{k}$ can be negative, the test becomes symmetric on the distribution, and hence for $k\geq 2$ we choose a significance level of $p=0.1$ in order to stay consistent with the 2-dependence test. The ”null” distributions and the threshold given by the significance p-value of rejection are illustrated in Figure 8d. If the observed values of $I_{k}$ are above or below these threshold values, we reject the null hypothesis.

In practice, a random generator is used to generate the random permutations (here the NumPy generator [62]), and the present method is not exempt from the possibility that it generates statistical dependences in the higher degrees.

Interpretation of the dependence test. The original interpretation of the test by Pethel and Hahs was that the null hypothesis corresponded to independent distributions, motivated by the statement that ”permutation destroys any dependence that may have existed between the datasets but preserves symbol frequencies”. However, considering simple analytical examples could not allow us to confirm their statement. We propose that for a given finite $m$, random permutations express all the possible statistical dependences that preserve symbol frequencies (cf. the discussion of E.Borel in [9]). This statement basically corresponds to what we observe in Figure 8. Hence we propose that in finite context the null-hypothesis corresponds to a random k-dependence. The meaning of the presented test is hence a selectivity or specificity test: a test of an $I_{k}$ of given k-tuple against a null hypothesis of ”randomly” selected k-statistical dependences that preserve the marginals and $m$.

4.6.1 Sampling size and graining landscapes - stability of minimum energy complex estimation

Figure 9 gives a first simple study of how robust the paths of maximum length are with respect to the variations of $m$ and $N$, in the case of the analysis of genes as variables. The limit $N\rightarrow\infty$ recovers Riemann integration theory and gives the differential entropy with the correcting additive factor $N$ (theorem 8.3.1 [13]).

The information paths of maximal length identified by our algorithm are relatively stable in the range of $N=5,7,9,11$ and $m=34,56,89,111$ where the $m$ cells were taken among the 111 neurons of population A. If we consider that the paths that only differ by the ordering of the variables are equivalent, then the stability of the two first paths is further and largely improved. The undersampling dimension obtained in these conditions is $k_{u}(m=34)=5,~{}k_{u}(m=56)=6,~{}k_{u}(m=89)=6,~{}k_{u}(m=111)=6$ and $k_{u}(N=5)=8,~{}k_{u}(N=7)=7,~{}k_{u}(N=9)=6,~{}k_{u}(N=11)=5$. In general, information landscapes can be investigated with the additional dimensions of $N$ and $m$ together with $n$. It allows to define our landscapes as iso-graining landscapes and to study the appearance of critical points in a way similar to what is done in thermodynamics. In practice, to study more precisely the variations of information depending on $N$ and $m$ and to obtain a 2-dimensional representation, we plot the mean information as a function of $N$ and $m$ together with $n$, as presented in Figure 10a. We call the obtained landscapes the iso-graining $I_{k}$ landscapes. The choice of a specific graining $N$ can be done using this representation: a ”pertinent” graining should be at a critical point of the landscape (a first minimum of an information path), consistent with the proposition of the work of Reshef and colleagues [46], who used maximal information coefficient ($MI_{2}C$) depending on the graining (with a more elaborated graining procedure) to detect pairwise associations. We have chosen to illustrate the landscapes with $N=9$ according to this criterion and the undersampling criterion, because the $I_{2}$ values are close to their maximal values and the sampling size is not too limiting, with a $k_{u}=6$ (see Figure 10a). Moreover, this choice of graining size $N=9$ is sufficiently far from the critical point to ensure that we are in the condensed phase where interactions are expected. It is well below the analog of the critical temperature (the critical graining size), which according to the Figure 10a happens at $N_{c}=3$ (the $N$ for which the critical points cease to be trivial). In general, there is no reason why there should be only one ”pertinent” graining.

The graining algorithm could be improved by applying direct methods of probability density estimation [50], or more promisingly persistent homology [15]. Finer methods of estimation (graining) have been developed by Reshef and colleagues [46] in order to estimate pairwise mutual-information, with interesting results. Their algorithm presents a lower computational complexity than the estimation on the lattice of partitions, but a higher complexity than the simple one applied here.

What we call the iso-sampling size $I_{k}$ landscapes is presented in Figure 10b for mean $I_{k}$. Such investigation is also important since it monitors what is usually considered as the convergence (or divergence) in probability of the informations. For the estimations below the $k_{u}$ represented here, the information estimations are quite constant as a function of $m$, indicating the stability of the estimation with respect to the sample size.

Acknowledgement

Acknowledgments: This work was funded by the European Research Council (ERC consolidator grant 616827 CanaloHmics to J.M.G.; supporting M.T. and P.B.) and Median Technologies, developed at Median Technologies and UNIS Inserm 1072 - Université Aix-Marseille, and at Institut de Mathématiques de Jussieu - Paris Rive Gauche (IMJ-PRG), and thanks previously to supports and hostings since 2007 of Max Planck Institute for Mathematic in the Sciences (MPI-MIS) and Complex System Instititute Paris-Ile-de-France (ISC-PIF). This work addresses a deep and warm acknowledgement to the researchers who helped its realization: G.Marrelec and J.P. Vigneaux; or supported-encouraged it: H.Atlan, F.Barbaresco, H.Bénali, P.Bourgine, F.Chavane, J.Jost, A.Mohammad-Djafari, JP.Nadal, J.Petitot, A.Sarti, J.Touboul.

Supplementary material, contributions and previous versions

Previous Version: a partial version of this work has been deposited in the method section of Bioarxiv 168740 in July 2017 and preprints [5].

Author contributions: D.B. and P.B. wrote the paper, P.B. analysed the data ; M.T. performed the experiments ; M.T. and J.M.G. conceived and designed the experiments ; D.B., P.B., M.T. and J.M.G participated in the conception of the analysis.

Supplementary material: The software Infotopo is available at https://github.com/pierrebaudot/INFOTOPO

Abbreviations

The following abbreviations are used in this manuscript:

[TABLE]

Appendix A Appendix: Bayes free energy and Information quantities

A.1 Parametric modelling

As we mentioned in the introduction, the statistical analysis of data $X$ is confronted to a serious the risk of circularity, because the confidence in the model is dependent on the probability law it assumes and reconstructs in part. Several approaches were followed to escape from this circularity; all of them rely on the choice of a set $\Theta$ of probability laws where $\mathbb{P}_{X}$ is researched. For instance, maintaining the frequentist point of view, the Fisher information metric on $\Theta$ (cf. [34]) determines bounds on the confidence. Another popular approach is to choose an a priori probability $\mathbb{P}_{\Theta}$ on $\Theta$, and to revise this choice after all the experiments $X(z),z\in Z$, by computing the probability on $E\times\Theta$, which better explains the results (the new probability on $\Theta$ is its marginal, and for each $\theta$ in $\Theta$, the probability $P_{\theta}$ on $E$ is its conditional probability). Here a more precise principle is necessary, which expresses a trade-off between the maximization of the marginal probability of the results under the constraint to be not too far from the prior. A popular example is the minimization of the Bayes Free energy $F_{V}(P)$, which appears as the maximum of entropy of the new a posteriori probability under the constraint to predict in the mean the data and to depart the less possible from the a priori probability on the probabilities. This function is given by a Kullback-Leibler distance $D_{KL}$. In the finite setting, with a uniform a priori, this consists in maximizing the entropy among the laws that predict the observed distribution. Remark that the two methods, Bayes and Fisher, are related, because in most cases the chosen a priori probability laws (and the data estimation) used in the function $F_{V}$ are given by frequencies, and because the distance $D_{KL}(P,Q)$ is approximated by the Fisher metric at $P$ when $Q$ approaches $P$.

A.2 Bethe approximation

Let us remind that for two probability laws $P,Q$ on the same finite set $\Omega$, the Kullback-Leibler divergence from $P$ to $Q$ is defined by

[TABLE]

Contrarily to its name, it is not a true distance, because it is not symmetric, however it is always positive and it is equal to zero if and only if $P=Q$. Another drawback is that it can be $+\infty$: it is so when there exists $x$ such that $Q_{x}=0$ but $P_{x}>0$, i.e. when $P$ is not absolutely continuous with respect to $Q$.

The Kullback-Leibler divergence permits to define the Bayes free energy functional as follows:

The unknown is the probability law $P_{b}$ on $E\times\Theta$.

[TABLE]

where $P_{a}(\theta)$ is the a priori on the probability laws and where $P_{L}(s)$ represents the new partial data, collected by a collection of variables $X_{L}$, and expressed by a probability law.

[TABLE]

This function looks like a free energy in Statistical Physics, that is the sum of the negentropy and the mean of an energy function.

Here we assume that $\Omega=E_{S}$ for a family of variables $S_{i},i=1,...,N$, and the states are the possible values of the joint variable $S$.

Due the strict convexity of the negentropy, $F_{V}$ has a unique minimum, that defines the equilibrium state.

Practically, the full entropy is difficult to estimate, thus approximations were introduced, following Bethe and Kikuchi (cf. Mori [40]), generalizing the Mean Field Theory. These approximations are no more convex in the unknown $P_{b}$, they are obtained by replacing the full entropy $H$ by a convenient linear combination of entropies of more accessible variables (observable quantities). It is here that the information functions $H_{k}$ and $I_{k}$ appear in the bayesian variational calculus (cf. Mori [40]):

Consider a simplicial complex $K$ in the simplex $\Delta([N])$, i.e. a collection of faces that contains every faces inside each face it contains, and assume $K$ a combinatorial ($PL$) manifold of dimension $d$, with possibly a boundary that is a combinatorial ($PL$) manifold $\partial K$; then the Bethe function associated to $K$ is given by the two equivalent following formulas:

[TABLE]

where the sum is taken over the set $K^{*}$ of faces not contained in $\partial K$, and $|I|$ denotes the dimension of the face $I$;

[TABLE]

where the sum is taken over all the faces of $K$, including the boundary, and $I_{|J|}(S_{J};Q)$ is the higher mutual information considered everywhere above in the text.