# S$\Omega$I: Score-based o-information Estimation

Mustapha Bounoua  Giulio Franzese  Pietro Michiardi

###### Abstract

The analysis of scientific data and complex multivariate systems requires
information quantities that capture relationships among multiple random
variables. Recently, new information-theoretic measures have been developed to
overcome the shortcomings of classical ones, such as mutual information, that
are restricted to considering pairwise interactions. Among them, the concept
of information synergy and redundancy is crucial for understanding the high-
order dependencies between variables. One of the most prominent and versatile
measures based on this concept is o-information, which provides a clear and
scalable way to quantify the synergy-redundancy balance in multivariate
systems. However, its practical application is limited to simplified cases. In
this work, we introduce s$\Omega$i, which allows — for the first time — to
compute o-information without restrictive assumptions about the system. Our
experiments validate our approach on synthetic data, and demonstrate the
effectiveness of s$\Omega$i in the context of a real-world use case.

Machine Learning, Mutual information, Score based models, Diffusion models

##  1 Introduction

Mutual Information (mi) is a fundamental measure which allows investigation of
the non-linear dependence between random variables (Shannon, 1948; MacKay,
2003). Despite its success in various domains, classical mi suffers from
limitations when analyzing systems composed by more than two variables. This
constitutes an important limitation, considering that many scientific
endeavors aim at an accurate statistical characterization of systems which are
composed of many random variables. Examples includes neuroscience (Latham &
Nirenberg, 2005; Ganmor et al., 2011; Gat & Tishby, 1998), climate models
(Runge et al., 2019), econometrics (Dosi & Roventini, 2019), and machine
learning (Tax et al., 2017), to name a few.

A recent attempt to overcome such limitations, and to extend the applicability
of information-theoretic tools to multivariate systems, is represented by
Partial Information Decomposition (pid) (Williams & Beer, 2010). The key idea
behind such method is the decomposition of the overall mi between a set of
source variables and a given target variable into non-negative constituents.
In particular, pid quantifies how much of the total information about the
target variable is encoded redundantly, synergistically or uniquely into given
subsets of variables. Redundancy quantifies information that is shared between
subsets of the partition, synergy describes the additional information that is
endowed to all subsets observed jointly but that is not available from
individual constituents of the partition, and uniqueness quantifies the
information that is lost when a given subset is not observed, removing the
amount of redundant and synergistic information associated to that subset. The
pid method requires partitioning the source system into all its possible
subsets and computing the information decomposition of all constituents with
respect to the target variable.

Despite its elegance, this measure is not without drawbacks. Indeed, there is
no consensus on the best way to define and compute pid, and several variants
have emerged, including (Barrett, 2014), who reformulate synergy and
redundancy for Gaussian systems (but that has been judged as poorly motivated
by (Venkatesh et al., 2023)), (Finn & Lizier, 2019), who use the algebraic
structure of information sharing, (Ay et al., 2019), who rely on cooperative
game theory, (Rosas et al., 2020), who build on concepts related to data
privacy and disclosure, (Kolchinsky, 2019), who use set theory, (van Enk,
2023), who deal with scalability issues by pooling probabilities, (Gutknecht
et al., 2023), who use a mereological formulation, and (Makkeh et al., 2021;
Ehrlich et al., 2023), who advocate for methods based on the exclusions of
probability mass. Nevertheless, the main limitations of pid persist in all
variants. Indeed, computational complexity grows extremely fast, precisely as
the Dedekind number of variables (which is more than $10^{31}$ for 9
variables). Moreover, pid computation relies on a partition of the system into
a set of sources and a unique target. This can be an artificial distinction
which limits usability and interpretability of the results. This latter
problem is partially addressed in (Varley et al., 2023), who introduce Partial
Entropy Decomposition (ped).

Motivated by these limitations, (Rosas et al., 2019) introduce the concept of
o-information, a measure which captures the synergy-redundancy dominance in
multivariate systems. In contrast to pid, this measure does not require the
system to be partitioned into sources and a target, and gracefully scales in
the number of its components (Martinez Mediano, 2022). Furthermore, recent
extensions such as o-information locality (Scagliarini et al., 2021) and
gradient computation (Scagliarini et al., 2023) allow a fine-grained analysis
of system behavior. However, o-information measures are accessible only in
restricted scenarios. Indeed, existing methods rely on estimation techniques
that requires either i) discrete distributions (or binning of continuous ones)
or ii) Gaussian distributions. In this work, we show — for the first time —
that such limitations can be lifted by using and extending recent methods to
estimate mi (Franzese et al., 2023a; Kong et al., 2023).

Our work is organized as follows: § 2 introduces the high-dimensional
interaction measures which we investigate in this work, while § 3 proposes
Score-based O-Information estimation (s$\Omega$i), our novel methodology which
allows scalable and flexible o-information estimation. § 4 validates
experimentally our proposed method, where we report a series of compelling
results on various synthetic systems, for which ground truth values are known
and accessible analytically. Furthermore, we consider a realistic endeavor by
revisiting previous studies (Venkatesh et al., 2023) that focus on the
analysis of brain activity in mice. Our method allows lifting previous
limiting assumptions, and allow synergy-redundancy characterizations that are
compatible with observations made by domain experts. Finally, we summarize our
findings in § 5.

##  2 High dimensional interaction measures

Consider the continuous multivariate random variable
$X=\\{X^{1},\dots,X^{N}\\}\sim p(x^{1},\dots,x^{N})$. We indicate the
collection of all but the $i_{th}$ random variable with the symbol
$X^{\setminus i}\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\\{X^{1},..,X^{i-1},X^{i+1},..,X^{N}\\}$. When necessary, we
indicate marginal and conditional distributions by properly specifying the
arguments of the distribution, e.g. $X^{i}\sim p(x^{i})$ or $X^{\setminus
i}\,|\,{X^{i}}\sim p(x^{1},\dots,x^{i-1},x^{i},\dots,x^{M}\,|\,x^{i})$.

A central quantity in this work is the Shannon entropy associated to a given
random variable $\mathcal{H}(X)\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\mathbb{E}\left[-\log p(X)\right]$ (Cover, 1999). Considering the
case of bi-variate (i.e. $N=2$) random variable $X$, entropy and conditional
entropy allow computation of the mutual information (mi) flow $\mathcal{I}$
between the two random variables $X^{1},X^{2}$:
$\mathcal{I}(X^{1};X^{2})=\mathcal{H}(X^{1})-\mathcal{H}(X^{1}|X^{2})$, where
$\mathcal{H}(X^{1}|X^{2})=\mathbb{E}\left[-\log p(X^{1}\,|\,X^{2})\right]$.
Importantly, such quantity can also be expressed as the Kullback-Leibler (kl)
divergence (Cover, 1999) between the joint and the product of marginal
distributions:
$\mathcal{I}(X^{1};X^{2})=\textsc{kl}\left[p(x^{1},x^{2})\;\|\;p(x^{1})p(x^{2})\right]$.
For the case of $N=3$, it is possible to define the mi as
$\mathcal{I}(X^{1};X^{2};X^{3})=\mathcal{I}(X^{1};X^{2})-\mathcal{I}(X^{1};X^{2}|X^{3})$,
where
$\mathcal{I}(X^{1};X^{2}|X^{3})=\mathcal{H}(X^{1}\,|\,X^{3})-\mathcal{H}({X^{1}\,|\,X^{2},X^{3}})$.
This quantity, also known as co-information or interaction information, can
counter-intuitively result in a negative value, and measures the difference
between synergistic and redundant interactions (Rosas et al., 2019).

Since, for $N>3$, interaction information becomes difficult to grasp (Williams
& Beer, 2010; Rosas et al., 2019), our goal in this work is to consider
extensions to mi, while preserving interpretability. In particular, a measure
of the interaction strengths in a system with $N>3$ can be obtained by
studying the average information between one variable and the rest of the
system:

$\mathcal{S}(X)\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\sum\limits_{i=1}^{N}\mathcal{I}(X^{i};X^{\setminus i}).$ (1)

This quantity, named s-information, can be decomposed into the redundant and
synergistic components of the considered multivariate system. In particular,
since $X^{\setminus i}=\\{X^{<i},X^{>i}\\}$, where
$X^{<i}=\\{X^{1},\dots,X^{i-1}\\}$ and $X^{>i}=\\{X^{i+1},\dots,X^{N}\\}$
(with $X^{>N}=\emptyset$), we can use the conditional mutual information laws
(Cover, 1999) and rewrite $\mathcal{S}(X)$ as:

$\mathcal{S}(X)=\sum\limits_{i=1}^{N}\mathcal{I}(X^{i};X^{>i})+\sum\limits_{i=1}^{N}\mathcal{I}(X^{i};X^{<i}\,|\,X^{>i}).$
(2)

The two positive series which constitute $\mathcal{S}(X)$ are equivalent to
the Total Correlation (tc) (Sun, 1975) and the Dual Total Correlation (dtc)
(Sun Han, 1980). Then, $\mathcal{S}(X)=\mathcal{T}(X)+\mathcal{D}(X)$, where
(proof in Appendix A)

$\displaystyle\mathcal{T}(X)=\sum_{i=1}^{N}\mathcal{H}(X^{i})-\mathcal{H}(X),$
(3)
$\displaystyle\mathcal{D}(X)=\mathcal{H}(X)-\sum_{i=1}^{N}\mathcal{H}(X^{i}|X^{\setminus
i}).$ (4)

tc is high in cases where, for each variable $X^{i}$, at least one of its
“children” (variables in $X^{>i}$) carries information about it. Importantly,
the number of children conveying information (whether 1, 2, or $N-1$) is
irrelevant. Since $\mathcal{T}(X)$ is permutation invariant, a high value
implies that for every ordering of the variables, and hence for all possible
combinations of children of a given variable, the average mutual information
between variables and their children remains high. This intuition, which
suggests redundancy, can similarly be obtained by considering the entropic
formulation. Indeed, whenever a system is composed of perfectly independent
variables ($X^{i}\perp X^{j},i\neq j$)
$\mathcal{H}(X)=\sum_{i=1}^{N}\mathcal{H}(X^{i})$ and consequently
$\mathcal{T}(X)=0$. On the other hand, a copy system ($X^{i}=X^{j},\forall
i,j$) achieves infinite $\mathcal{T}(X)$, as $\mathcal{H}(X)=-\infty$, since
the support of the joint distribution is on a lower than $N-$dimensional
space. tc also admits a representation in terms of kl divergences,
$\mathcal{T}(X)=\textsc{kl}\left[p(x)\;\|\;\prod\limits_{i=1}^{N}p(x^{i})\right]$,
which we will exploit later in our proposed methodology.

Similar considerations can be carried out for the dtc. Consider a single mi
term $I(X^{i};X^{<i}|X^{>i})$. The focus of this conditioning is about
quantifying how much additional information the variables $X^{<i}$ carry about
$X^{i}$ if we are also given access $X^{>i}$. Whenever the variables are
independent or redundant (the copy system), this value is identically zero.
However, whenever the aid of the extra measurements unlocks new bits of
information, which suggests a synergistic scenario, its value is positive.

Having recognized that the average interaction $\mathcal{S}(X)$ in a
multivariate system can be decomposed into measures of redundancy
$\mathcal{T}(X)$ and synergy $\mathcal{D}(X)$, we can introduce a new
information theoretic measure which quantifies the difference between the two
behaviours. This quantity, named o-information (Rosas et al., 2019), is
defined as

$\Omega(X)=\mathcal{T}(X)-\mathcal{D}(X).$ (5)

In summary, while s-information only quantifies the strength of interactions
in a system, o-information also determines the nature of these interactions,
being them redundant or synergistic. Intuitively, a redundancy-dominated
system is the most parsimonious explanation — in an Occam’s razor sense —
whenever $\Omega(X)>0$. Conversely, a negative value $\Omega(X)<0$ is
associated with a synergy-dominated system. o-information is a natural
generalization of mi for more than 3 variables: indeed, it is equal to the co-
information for $N=3$, and is a measure which preserves interpretability for
any positive $N$.

One important property of o-information is that it gracefully scales with the
number of random variables composing a system, as opposed to, e.g. the pid
measure, which has much worse scalability.

Since o-information measures the overall information dynamics among variables,
recent work focus on ways to study the individual influence of variables to
the high-order interactions, and capture the interaction structure of a
multivariate system (Scagliarini et al., 2023). The first order difference,
called the gradient of o-information, captures how much o-information changes
when adding or removing a given system variable $i$:

$\partial_{i}\Omega(X)=\Omega(X)-\Omega(X^{\setminus i}).$ (6)

A positive value implies that $X^{i}$ provides redundant information to the
system, while a negative one suggests that its interaction with other
variables is mainly synergistic.

##  3 Score-based o-information estimation

o-information and its gradient represent extremely useful information
theoretic measures to study multivariate systems. However, as it is clear from
Equations 3, 4 and 5, their estimation requires access to entropies,
conditional entropies and kl divergence measures. When strict assumptions
about the distribution of variables composing the system are possible, such as
discrete or Gaussian distributions, existing implementations of o-information
estimators have been used successfully in a number of application domains
(Varley et al., 2022; Sparacino et al., 2023; Stramaglia et al., 2021;
Chiarion et al., 2023). However, in more realistic cases where such
assumptions are not valid, there currently does not exist a method to estimate
the constituents of o-information in a reliable and scalable manner. In this
work, we present the first methodology allowing estimation of o-information
for more general scenarios. Our method unfolds according to the observation
that all quantities of interest can be expressed in terms of kl divergences,
and relies on a technique to estimate such divergences which scales gracefully
with the system size. Our key ingredient is the score function associated to
data distributions (Vincent, 2011; Song & Ermon, 2019) and the method we
present leverages recent advances in the field of mi estimation (Franzese et
al., 2023a; Kong et al., 2023).

###  3.1 Score-based divergence estimation

Consider the generic multivariate random variable $X$ with associated
distribution $p(x)$. Provided that certain minimal regularity assumptions are
met (Vincent, 2011), it is always possible to associate the distribution
$p(x)$ to its score function, defined as the gradient of its logarithm,
$\nabla\log p(x)$.

Recently, the community has showed tremendous interest (Song & Ermon, 2019;
Song et al., 2021) in a generalization of such concept, which involves
computing the score function of a noised version of the variable $X$, due to
the possibility of adopting such concept for generative modelling purposes.
Accordingly, in this work we define a noised version of the variable $X$ with
corresponding intensity indexed by $t\in[0,\infty)$. Then, the new variable is
constructed as $X_{t}=X+\sqrt{2t}W$, where $W$ is an $N-$dimensional Gaussian
random vector with zero mean and identity covariance matrix.

This new random variable can be associated to its time-varying score function
$\nabla\log p_{t}(x)$. In particular the analytic expression of $p_{t}(x)$ can
be obtained as the solution of the Partial Differential Equation (pde)
$\frac{\differential p_{t}(x)}{\differential t}=\Delta p_{t}(x)$, with initial
conditions given by $p_{0}(x)=p(x)$.

Next, we consider the kl divergence between two generic distributions and
define how it can be computed using score functions, a result which we will
use later for computing o-information.

######  Proposition 1.

(Franzese et al., 2023a; Kong et al., 2023) The kl divergence between two
generic distributions $p(x)$ and $q(x)$, defined as

$\textsc{kl}\left[p(x)\;\|\;q(x)\right]=\int
p(x)\log\frac{p(x)}{q(x)}\differential x,$

can be computed considering the time-varying score functions
$\nabla\log(p_{t})$ and $\nabla\log(q_{t})$, according to the following
expression:

$\textsc{kl}\left[p(x)\;\|\;q(x)\right]=\int
p_{t}(x)\norm{\nabla\log(\frac{p_{t}(x)}{q_{t}(x)})}^{2}\differential
x\differential t.$

Proof sketch. To avoid clutter, we drop the dependence on $x$ of the
distributions. Let’s define $r_{t}\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\int p_{t}\log\frac{p_{t}}{q_{t}}\differential x$.

Since it holds that
$r_{\infty}-\textsc{kl}\left[p\;\|\;q\right]=\int_{0}^{\infty}\frac{\differential
r_{t}}{\differential t}\differential t$, we need

$\int\frac{\differential r_{t}}{\differential t}\differential
t=\int\frac{\differential p_{t}}{\differential
t}\log(\frac{p_{t}}{q_{t}})+p_{t}\frac{\differential}{\differential
t}\log(\frac{p_{t}}{q_{t}})\differential x\differential t.$

Note that $\int p_{t}\frac{\differential}{\differential
t}\log(\frac{p_{t}}{q_{t}})\differential x\differential
t=\int\frac{\differential}{\differential t}p_{t}-\frac{p_{t}}{q_{t}}\Delta
q_{t}\differential x\differential t$, and
$\int\frac{\differential}{\differential t}p_{t}\differential x\differential
t=0$. Then, the expression above can be rewritten as $\int
p_{t}\Delta\log(\frac{p_{t}}{q_{t}})-\frac{p_{t}}{q_{t}}\Delta
q_{t}\differential x\differential t$. Integrating by parts we obtain
$\int-\nabla
p_{t}\nabla\log(\frac{p_{t}}{q_{t}})+\nabla(\frac{p_{t}}{q_{t}})\nabla
q_{t}\differential x\differential t$. Since $\nabla p_{t}=p_{t}\nabla\log
p_{t}$ and $\nabla(\frac{p_{t}}{q_{t}})\nabla q_{t}=p_{t}\nabla\log
q_{t}\nabla\log(\frac{p_{t}}{q_{t}})$, and $r_{\infty}=0$ (Franzese et al.,
2023b; Villani, 2009; Collet & Malrieu, 2008), the proposition follows. ∎

The result in Proposition 1 allows, in principle, the exact computation of kl
divergences, provided knowledge of the score functions
$\nabla\log{p_{t}},\nabla\log{q_{t}}$. Such knowledge is however out of reach
in practical cases, which is why in this work we consider a parametric
approximation of such vector fields, leading to a kl divergence estimator. In
particular, we leverage the methodology considered in (Song & Ermon, 2019;
Song et al., 2021) where the parametric score $s_{t}$ is obtained by
minimizing the so called denoising score-matching loss

$\int
p(x)p_{0t}(\tilde{x}\,|\,x)\norm{s_{t}(\tilde{x})-\nabla\log(p_{0t}(\tilde{x}\,|\,x))}^{2}\differential
x\differential\tilde{x}\differential t,$

where $p_{0t}(\tilde{x}\,|\,x)$ is the conditional distribution of the noised
random variable given initial conditions $X=x$, i.e. $p_{t}(\tilde{x})=\int
p_{0t}(\tilde{x}\,|\,x)p(x)\differential x$. Note that $p_{0t}$ has known
Gaussian distribution with mean $x$ and variance $2t$. This allows, together
with the knowledge of the score functions, the implementation of an estimator
for the kl divergence.

Informally, learning the score can be understood as learning to denoise the
variable $X_{t}$ to obtain $X$. Indeed, the score functions have analytic
expression $\nabla\log{p_{t}(x)}=\frac{\mathbb{E}[X\,|\,X_{t}=x]-x}{2t}$,
where the only unknown is $\mathbb{E}[X\,|\,X_{t}=x]$. An alternative, but
equivalent parametrization of the problem, consists in estimating the noise
$W$, given $X_{t}$. We use this approach in our work since is considered to be
more stable numerically (Ho et al., 2020).

###  3.2 Estimating o-information

Armed with Proposition 1, we can leverage score functions to estimate the
information-theoretic quantities introduced in § 2. Here we consider an
extension of the simple noising process described in § 3.1, where we allow i)
noising of only certain subsets of the variables or ii) deletion of subset of
variables. In practice, the first case corresponds to learning to denoise a
portion of the variables, given auxiliary information about the other
(noiseless) variables, e.g. to learn
$\mathbb{E}[X^{i}\,|\,X^{i}_{t}=\tilde{x}^{i},X^{\setminus i}=x^{\setminus
i}]$. Instead, the second case amounts to denoising problems akin to
$\mathbb{E}[X^{i}\,|\,X^{i}_{t}=\tilde{x}^{i}]$. In our implementation, we
follow the approach proposed in (Bounoua et al., 2023) (See Appendix B). Next,
we use such an intuition to derive a series of propositions that pave the way
to o-information computation.

In what follows, we use the compact notation
$\begin{bmatrix}(\cdot)^{i}\end{bmatrix}_{i=1}^{N}$, to indicate a
concatenation of $N$ elements in a column vector.

######  Proposition 2.

Given a multivariate random variable $X=\\{X^{1},\dots,X^{N}\\}\sim
p(x^{1},\dots,x^{N})$, and its corresponding noised version, the Total
Correlation $\mathcal{T}(X)$ is equal to:

$\int\frac{1}{4t^{2}}\mathbb{E}\norm{\mathbb{E}[X\,|\,X_{t}]-\begin{bmatrix}\mathbb{E}[X^{i}\,|\,X^{i}_{t}]\end{bmatrix}_{i=1}^{N}}^{2}\differential
t.$

Proof Sketch. Recall that
$\mathcal{T}(X)=\textsc{kl}\left[p(x)\;\|\;\prod\limits_{i=1}^{N}p(x^{i})\right]$.
Then, by virtue of Proposition 1, we have that $\mathcal{T}(X)$ equals

$\int p_{t}(x)\norm{\nabla\log
p_{t}(x)-\begin{bmatrix}\frac{\partial}{\partial x^{i}}\log
p_{t}(x^{i})\end{bmatrix}_{i=1}^{N}}^{2}\differential x\differential t.$

The terms $\frac{\partial}{\partial x^{i}}\log p_{t}(x^{i})$ correspond to
$\nicefrac{{1}}{{2t}}(\mathbb{E}[X^{i}\,|\,X^{i}_{t}=x^{i}]-x^{i})$. Then, the
proposition follows. ∎

######  Proposition 3.

Given a multivariate random variable $X=\\{X^{1},\dots,X^{N}\\}\sim
p(x^{1},\dots,x^{N})$, and its corresponding noised version, the s-information
$\mathcal{S}(X)$ is equal to:

$\int\frac{1}{4t^{2}}\mathbb{E}\norm{\begin{bmatrix}\mathbb{E}[X^{i}\,|\,X^{i}_{t}]\end{bmatrix}_{i=1}^{N}-\begin{bmatrix}\mathbb{E}[X^{i}\,|\,X^{i}_{t},X^{\setminus
i}]\end{bmatrix}_{i=1}^{N}}^{2}\differential t.$

Proof Sketch. In light of Equation 1, it holds that

$\mathcal{S}(X)=\sum_{i=1}^{N}\int p(x^{\setminus
i})\textsc{kl}\left[p(x^{i}\,|\,x^{\setminus
i})\;\|\;p(x^{i})\right]\differential x^{\setminus i},$

where the $i_{th}$ kl term of the sum is equal to (Proposition 1)

$\displaystyle\int p(x^{i}\,|\,x^{\setminus
i})p_{0t}(\tilde{x}^{i}\,|\,x^{i})$
$\displaystyle\norm{\frac{\partial}{\partial\tilde{x}^{i}}\log(\frac{p_{t}(\tilde{x}^{i})}{\hat{p}_{0t}(\tilde{x}^{i}\,|\,x^{\setminus
i})})}^{2}\differential\tilde{x}^{i}\differential x^{i}\differential t.$

Now, we can move the terms $p(x^{\setminus i})$ inside the kl computation
integrals and write the sum of the norms as the norm of a vector, which allows
computing s-information as

$\displaystyle\mathcal{S}(X)=\int p(x)p_{0t}(\tilde{x}\,|\,x)$
$\displaystyle\norm{\begin{bmatrix}\frac{\partial}{\partial\tilde{x}^{i}}\log
p_{t}(\tilde{x}^{i})\end{bmatrix}_{i=1}^{N}-\begin{bmatrix}\frac{\partial}{\partial\tilde{x}^{i}}\log
p_{t}(\tilde{x}^{i}\,|\,x^{\setminus
i})\end{bmatrix}_{i=1}^{N}}^{2}\differential\tilde{x}\differential
x\differential t,$

where ${p}_{t}(\tilde{x}^{i}\,|\,x^{\setminus
i})=\int{p}_{0t}(\tilde{x}^{i}\,|\,x^{i}){p}({x}^{i}\,|\,x^{\setminus
i})\differential x^{i}$.

Finally, the proposition follows since we can interpret the elements inside
the square norm in terms of denoisers, with
$\frac{\partial}{\partial\tilde{x}^{i}}\log
p_{t}(\tilde{x}^{i}\,|\,x^{\setminus
i})=\nicefrac{{1}}{{2t}}(\mathbb{E}[X^{i}\,|\,X^{i}_{t}=\tilde{x}^{i},X^{\setminus
i}=x^{\setminus i}]-\tilde{x}^{i})$∎

######  Proposition 4.

Given a multivariate random variable $X=\\{X^{1},\dots,X^{N}\\}\sim
p(x^{1},\dots,x^{N})$, and its corresponding noised version, the Dual Total
Correlation $\mathcal{D}(X)$ equals:

$\int\frac{1}{4t^{2}}\mathbb{E}\norm{\mathbb{E}[X\,|\,X_{t}]-\begin{bmatrix}\mathbb{E}[X^{i}\,|\,X^{i}_{t},X^{\setminus
i}]\end{bmatrix}_{i=1}^{N}}^{2}\differential t$

Proof Sketch. The starting point to obtain dtc is to recall that
$\mathcal{D}(X)=\mathcal{S}(X)-\mathcal{T}(X)$. Then, it is sufficient to
expand the square norms of $\mathcal{S}(X)$ and $\mathcal{T}(X)$ and combine
the different terms, to state that $\mathcal{D}(X)$ equals:

$\displaystyle\int p(x)p_{0t}(\tilde{x}\,|\,x)$ $\displaystyle\norm{\nabla\log
p_{t}(\tilde{x})-\begin{bmatrix}\frac{\partial}{\partial\tilde{x}^{i}}\log
p_{t}(\tilde{x}^{i}\,|\,x^{\setminus
i})\end{bmatrix}_{i=1}^{N}}^{2}\differential\tilde{x}\differential
x\differential t.$

This can be proven considering that i)
$\mathbb{E}\left[\mathbb{E}[X^{i}\,|\,X_{t}]]\mathbb{E}[X^{i}\,|\,X^{i}_{t}]\right]=\mathbb{E}[(\mathbb{E}[X^{i}\,|\,X^{i}_{t}])^{2}]$
ii)
$\mathbb{E}\left[\mathbb{E}[X^{i}\,|\,X^{i}_{t}]\mathbb{E}[X^{i}\,|\,X^{i}_{t},X^{\setminus
i}]\right]=\mathbb{E}\left[(\mathbb{E}[X^{i}\,|\,X^{i}_{t}])^{2}\right]$ and
iii)
$\mathbb{E}\left[\mathbb{E}[X^{i}\,|\,X_{t}]]\mathbb{E}[X^{i}\,|\,X^{i}_{t},X^{\setminus
i}]\right]=\mathbb{E}[(\mathbb{E}[X^{i}\,|\,X^{i}_{t}])^{2}]$. Then, the
proposition follows. ∎

Finally, to estimate o-information, it is sufficient to combine Proposition 2
and Proposition 4, and apply Equation 5. In practical terms, our method
requires access to denoisers for the three following scenarios: i) given
$X_{t}$ estimate $X$ ii) given $X_{t}^{i}$ estimate $X^{i}$ iii) given
$X_{t}^{i}$ and $X^{\setminus i}$ estimate $X^{i}$. To achieve this, we extend
the methodology proposed in (Bounoua et al., 2023), and amortize the three
different scenarios with a unique denoising network, which takes as input the
concatenation of noised and clean variables and outputs the corresponding
estimates (see Appendix B).

##  4 Experimental validation

We evaluate our method according to two strategies. First, we focus on a
synthetic setup that allows analytic computation of o-information and full
control on system scale. Then, we consider real data collected in a study of
brain activity in mice, to demonstrate how s$\Omega$i unlocks new avenues in
the application of information measures in real systems without the need for
restrictive assumptions.

###  4.1 Synthetic benchmark

We consider a canonical Gaussian system, whereby we control the number of
variables describing the system, their dimensionality, the inter-dependencies
between variables describing how they interact, and the strength of
interaction. Despite prior work such as (Venkatesh et al., 2023) also consider
the canonical Gaussian case, since it is restricted to such assumption by
construction, to the best of our knowledge, no other neural estimator capable
of estimating o-information for continuous variables exists. For this reason,
next, we construct an original baseline that also relies on neural estimation
of mi.

#### Baseline.

Recent work (Bai et al., 2023) describes a method to compute tc by noticing it
can be decomposed into a number of pairwise mi terms. Clearly, dtc can also be
decomposed into mi terms. Therefore, we extend (Bai et al., 2023) such that it
can be used as a baseline to compute o-information. The main limitation of
such a baseline is poor scalability: it requires training an individual model
for each mi term in which tc and dtc are decomposed in. We adopt the linear-
decomposition method (Bai et al., 2023), which results in $2(N-1)$ mi terms
(see Appendix C), and propose four variants to estimate mi based on (Belghazi
et al., 2018; Nguyen et al., 2007; Oord et al., 2018; Cheng et al., 2020). We
label this baseline approach according to the mi estimators: mine,nwj,
infonce, and club.

#### Experimental protocol

For each experiment, we use $100k$ samples for training the various neural
estimators, and $10k$ samples at inference time, to estimate o-information.
For our method s$\Omega$i, we use the vp-sde formulation (Song et al., 2021)
and learn a unique denoising network to estimate the various score terms. The
denoiser is a simple, stacked multilayer perceptron (mlp) with skip
connections, adapted to the input dimension. We apply importance sampling
(Huang et al., 2021; Song et al., 2021) at both training and inference time.
Finally, we use 10-sample Monte Carlo estimates for computing integrals. More
details about the implementation are included in Appendix C. For the baseline
variants, for each mi term we use an mlp that is sufficiently expressive given
the data dimension. All results are averaged over 5 seeds. Additional results
are included in Appendix D.

Our experiments unfold according to three inter-dependency scenarios, for
systems characterized by either redundancy, synergy or a mix of both
interactions.

#### Redundancy benchmark.

We consider $R=\mathcal{N}(0,\mathbb{I})$ as the redundant information
component in the system. All system variables are of the form
$X=\\{X^{1},\dots,X^{N}\\}=\\{R+\epsilon_{i},\dots,R+\epsilon_{N}\\}$, where
$\epsilon_{i}\sim\mathcal{N}(0,\sigma\mathbb{I})$ are mutually independent
random noise samples with standard deviation $\sigma$. We use $\sigma$ to
modulate the redundancy level: higher noise levels decrease the strength of
redundant interaction, and this has an impact on the value of o-information.

Next, we discuss results for a system with $N=10$ variables, organized as 3
redundant subsystem, each defined as described above. Figure 1 illustrates,
for various variable dimension, ranging from 5 to 20 dimensional Gaussians,
the ground-truth and the estimated o-information, for s$\Omega$i and the
various baselines. In this scenario, s$\Omega$i and baseline competitors
produce fairly accurate o-information estimates, when the dimensionality of
each random variable is small. When the dimension of systems variables grows,
however, the performance of the baseline methods degrades considerably. This
is due to the inherent limitations of the pairwise neural mi estimators, that
struggle with high dimensional data (Czy.z et al., 2023). Instead, the
performance of s$\Omega$i remains stable when increasing variable dimension,
and o-information estimates are accurate, even when interaction strength is
high.

(a) Dim = 5

(b) Dim = 10

(c) Dim = 15

(d) Dim = 20

Figure 1: Redundant system with 10 variables, organized into subsets of sizes
$\\{3,3,4\\}$ and increasing interaction strength.

#### Synergy benchmark.

In this case, we synthesize synergistic inter-dependency among system
variables by considering the following setup. For simplicity, consider three
random variables that behave as follows:

$\displaystyle X^{1}\sim\mathcal{N}(0,\mathbb{I}),\quad X^{2}=X^{1}+S$
$\displaystyle X^{3}=S+\epsilon,\quad\epsilon\sim\mathcal{N}(0,\sigma)\text{
and }X^{1}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X^{3}$

with $S\sim\mathcal{N}(0,\mathbb{I})$.

When $\sigma=0$, the synergy emerges through the Markov chain
$\\{X^{2},X^{3}\\}-X^{1}$, $\\{X^{1},X^{3}\\}-X^{2}$ and
$\\{X^{1},X^{2}\\}-X^{3}$, since no element alone is sufficient to recover the
remaining variables. We modulate $\sigma$ to achieve different synergistic
strengths. More generally, we simulate $N$ synergistic variables as:
$X^{1}\sim\mathcal{N}(0,\mathbb{I})$, $X^{2}=X^{1}+S_{1}+...+S_{N-2}$ and
$X^{i}=S_{i-2}+\epsilon_{i-2}\forall i\in\\{3,..,N\\}$.

Results in Fig. 2 show that s$\Omega$i achieves consistent results in all
scenarios, whereas the baselines behave poorly. Indeed, a synergy-only setting
is challenging, as it’s dominated by high dtc values required to capture high-
order interactions, on which the baselines based on pairwise mi estimator
fail.

(a) Dim = 5

(b) Dim = 10

(c) Dim = 15

(d) Dim = 20

Figure 2: Synergistic system with 10 variables, organized into subsets of
sizes $\\{3,3,4\\}$ and increasing interaction strength.

#### Mixed benchmark.

In general, systems components are characterized by a mix of redundant and
synergistic interactions. Then, we synthesize such a system by creating
subgroups dominated by redundancy and synergy, respectively, following the
procedures defined above.

Results in Fig. 3, demonstrate that our method s$\Omega$i stands out as the
best estimator in this challenging scenario. Baseline methods produce poor
estimates, especially when the synergistic interaction is dominant. Note that
s$\Omega$i reports a negative o-information whenever the system is synergy-
dominant and also succeeds in capturing interaction strengths, when the system
equilibrium changes in favor redundant interactions, by estimating correctly a
positive o-information.

(a) Dim = 5

(b) Dim = 10

(c) Dim = 15

(d) Dim = 20

Figure 3: Mixed-interaction system with 10 variables, organized into 2
redundancy-dominant subsets of size $\\{3,4\\}$ variables and one synergy-
dominant subset with $3$ variables. o-information is modulated by fixing the
synergy inter-dependency and increasing the redundancy.

#### Discussion.

We attribute the superior performance of s$\Omega$i, compared to the
baselines, to several factors. Score-based estimators have shown to be
extremely successful in fitting complex distributions, for example in the
context of generative modeling (Song & Ermon, 2020; Song et al., 2021).
Moreover, our technique relies on Proposition 1, whereby the difference of
score functions has been shown to produce an accurate estimate of kl
divergences, due to canceling effects of estimation errors (Franzese et al.,
2023a). Note also that the baselines we adopt in our work use mi estimators
that produce a bound only. Moreover, using individual models to estimate
several mi terms can naturally suffer from cumulative bias, which is avoided
in our case by amortizing computation with a single neural network.

#### Gradient of O-information

While o-information provides global information about dominance of either
synergy or redundancy, the contribution of individual variables to either
effects is not available. Next, we rely on the gradient of o-information to
study individual system components, as introduced in § 2. Indeed, our method
s$\Omega$i can be easily extended to output such gradients, by estimating
additional score functions, as described in Appendix B. In Figure 4, we
illustrate gradients of o-information applied to the mixed benchmark scenario
discussed above. While o-information of the whole system can be positive due
to the redundancy strength of some subgroup of variables, we notice that three
variables report a negative gradient, which is indicative of their synergistic
interaction. In Figure 4, ground truth gradient values are showed using a
diamond marker. Our estimator, despite suffering from some bias, correctly
attributes the role and interaction type of each system constituent.

(a) Dim = 5

(b) Dim = 10

(c) Dim = 5

(d) Dim = 10

Figure 4: Gradient of o-information for the mixed benchmark, for a system of 6
variables, and a system of 10 variables, and different dimension of variables.

###  4.2 Application to a real system

Multivariate analysis is a powerful tool for the field of neuroscience, as it
allows scientists to analyze activity patterns of different brain regions.
Understanding how the brain processes and transmits information during
different stimulus requires analysing the underlying inter-dependencies
between different brain regions. To show that s$\Omega$i is an effective tool
also in practical use cases, we now consider the Visual Behavior project,
which used the Allen Brain Observatory to collect a highly standardized
dataset consisting of recordings of neural activity in mice that have learned
to perform a visually guided task (Allen-Institute, 2022).

A visual change-detection task experiment was conducted on 80 mice using six
neuropixels probes tasked to report the activity of different regions of the
visual cortex. During the recordings, a set of 8 natural scenes were presented
in 250 ms flashes, at intervals of 750 ms. The same image was shown during
several flashes before a change to a new image. The mouse had to perform an
action to receive a water reward when the image changed. Ultimately, the
purpose of this experiment is to investigate how the different brain region of
the mice react to different types of stimulus, such as detecting a new image
(change) or not (no change).

In this work, we follow the prepossessing procedure described by (Venkatesh et
al., 2023), where in each experimental session, good quality units from each
area are chosen (See Appendix C). For each trial, the recorded spikes are
binned in 50 ms intervals, starting from the stimulus flash. We consider two
types of flashes: change and no change. For both cases, s$\Omega$i is used for
each time bin to estimate O-information (o-information). The reported
estimation is done using 10 Monte Carlo integration steps and averaged over
multiple seeds. We first consider three visual cortex regions VISp, VISl and
VISal, as done in (Venkatesh et al., 2023). We then extend the experiment to
six brain regions by including VISrl, VISam and VISpm.

We show our results in Figure 5, where the distribution of o-information
values are reported as box-plots for each bin. We remark that values of
o-information are higher in cases of change stimulus, and lower for the no
change stimulus. This suggests that higher amount of redundant information in
the visual cortex regions is transmitted in case of a flash with new scene.
Interestingly, when considering six areas of the visual cortex, our
observations remain valid, suggesting that the measured behaviour is common to
these other brain areas as well. Our results are aligned with (Venkatesh et
al., 2023). However, prior work rely on the pid measure, which requires the
brain regions to be artificially organized into two areas and a target
variable, due to scalability issues affecting pid. Our work confirms that
s$\Omega$i does not have such a limitation and allows a single estimation
procedure to obtain the same conclusions.

(a) 3 areas

(b) 6 areas

Figure 5: o-information estimate in the visual cortex region activity after
two types of stimulus flash across 72 trial sessions. Top: Analysis using
three brain region areas, Bottom: Extended analysis using six brain region
areas. The step size is set to $2ms$ which results in 25 dimensional data for
each bin per area. Different step sizes led to the same behavior (see Appendix
D).

##  5 Conclusion

We addressed the problem of analyzing multivariate systems, whereby the
essence of complexity does not only lie in the nature of the individual system
components, but also in the structure of their inter-dependencies. Indeed, the
analysis of high-order interaction among variables has emerged as an important
tool to deepen our understanding of such complex systems, with application
domains including machine learning, neuroscience, climate modeling, and many
more.

Recently, the scientific community has spent considerable effort on extending
information theory to allow the study of complex, multivariate systems
according to notions of uniqueness, redundancy and synergy. While no consensus
exists yet, on a information measure that can fully and reliably characterize
high-order interactions, in this work we focused on o-information, which has
desirable properties such as interpretability and scalability in number of
variables. The current state of the art is however at a roadblock. No
technique exists, to the best of our knowledge, to compute information
theoretic measures for arbitrary, continuous variables, which are common in
real world scenarios. Then, the endeavour of our work was to present a method
to lift such limitations, and endow practitioners and scientists with a
flexible and reliable tool to study complex systems associated to natural
phenomena.

In this paper, we proposed s$\Omega$i, a novel technique that leverages recent
neural estimators of mutual information and uses score functions of joint and
conditional distributions to compute divergences. We showed that s$\Omega$i
can compute o-information by training a single parametric function, which is
efficient and flexible. We validated our technique with a comprehensive
experimental protocol, both in synthetic and realistic settings. We
demonstrated that s$\Omega$i is accurate and robust across different system
configurations and complexities. We also applied s$\Omega$i to a case study of
mice brain activity, where we obtained plausible and interpretable results,
and showcased the scalability of s$\Omega$i to handle larger systems than
previously possible. We believe that our work contributes to a substantial
advancement of information measures computation and their applications to
real-world, complex systems.

## Impact Statement

This paper presents work to improve current methods to compute information
measures of complex systems, modeled as ensembles of multiple random
variables. Such information measures have been recently brought to the
attention of the scientific community, for their potential in explaining the
high-order interactions between systems part, and specifically to understand
information redundancy, uniqueness and synergy. Applications of such measures
range from multi-modal machine learning, neuroscience, climate modeling and
many more. There are many potential societal consequences of our work, none
which we feel must be specifically highlighted here.

## References

  * Allen-Institute (2022) Allen-Institute.  Visual behavior neuropixels dataset overview.  2022\.  URL https://portal.brain-map.org/explore/circuits/visual-behavior-neuropixels.
  * Ay et al. (2019) Ay, N., Polani, D., and Virgo, N.  Information decomposition based on cooperative game theory.  _ArXiv_ , abs/1910.05979, 2019.  URL https://api.semanticscholar.org/CorpusID:204512236.
  * Bai et al. (2023) Bai, K., Cheng, P., Hao, W., Henao, R., and Carin, L.  Estimating total correlation with mutual information estimators.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 2147–2164. PMLR, 2023.
  * Barrett (2014) Barrett, A. B.  An exploration of synergistic and redundant information sharing in static and dynamical gaussian systems.  _CoRR_ , abs/1411.2832, 2014.  URL http://arxiv.org/abs/1411.2832.
  * Belghazi et al. (2018) Belghazi, M. I., Baratin, A., Rajeshwar, S., Ozair, S., Bengio, Y., Courville, A., and Hjelm, D.  Mutual information neural estimation.  In _Proceedings of the 35th International Conference on Machine Learning_ , 2018.
  * Bounoua et al. (2023) Bounoua, M., Franzese, G., and Michiardi, P.  Multi-modal latent diffusion.  _arXiv preprint arXiv:2306.04445_ , 2023.
  * Cheng et al. (2020) Cheng, P., Hao, W., Dai, S., Liu, J., Gan, Z., and Carin, L.  Club: A contrastive log-ratio upper bound of mutual information.  06 2020.
  * Chiarion et al. (2023) Chiarion, G., Sparacino, L., Antonacci, Y., Faes, L., and Mesin, L.  Connectivity analysis in eeg data: A tutorial review of the state of the art and emerging trends.  _Bioengineering_ , 10, 2023.  URL https://api.semanticscholar.org/CorpusID:257644289.
  * Collet & Malrieu (2008) Collet, J.-F. and Malrieu, F.  Logarithmic sobolev inequalities for inhomogeneous markov semigroups.  _ESAIM: Probability and Statistics_ , 12:492–504, 2008.
  * Cover (1999) Cover, T. M.  _Elements of information theory_.  John Wiley & Sons, 1999.
  * Czy.z et al. (2023) Czy.z, P., Grabowski, F., Vogt, J. E., Beerenwinkel, N., and Marx, A.  Beyond normal: On the evaluation of mutual information estimators.  _ArXiv_ , abs/2306.11078, 2023.  URL https://api.semanticscholar.org/CorpusID:259202453.
  * Dosi & Roventini (2019) Dosi, G. and Roventini, A.  More is different… and complex! the case for agent-based macroeconomics.  _Journal of Evolutionary Economics_ , 29:1–37, 2019.
  * Ehrlich et al. (2023) Ehrlich, D. A., Schick-Poland, K., Makkeh, A., Lanfermann, F., Wollstadt, P., and Wibral, M.  Partial information decomposition for continuous variables based on shared exclusions: Analytical formulation and estimation.  _arXiv preprint arXiv:2311.06373_ , 2023.
  * Finn & Lizier (2019) Finn, C. and Lizier, J. T.  Generalised measures of multivariate information content.  _Entropy_ , 22, 2019.  URL https://api.semanticscholar.org/CorpusID:202888469.
  * Franzese et al. (2023a) Franzese, G., Bounoua, M., and Michiardi, P.  Minde: Mutual information neural diffusion estimation.  2023a.
  * Franzese et al. (2023b) Franzese, G., Rossi, S., Yang, L., Finamore, A., Rossi, D., Filippone, M., and Michiardi, P.  How much is enough? a study on diffusion times in score-based generative models.  _Entropy_ , 2023b.
  * Ganmor et al. (2011) Ganmor, E., Segev, R., and Schneidman, E.  Sparse low-order interaction network underlies a highly correlated and learnable neural population code.  _Proceedings of the National Academy of sciences_ , 108(23):9679–9684, 2011.
  * Gat & Tishby (1998) Gat, I. and Tishby, N.  Synergy and redundancy among brain cells of behaving monkeys.  _Advances in neural information processing systems_ , 11, 1998.
  * Gutknecht et al. (2023) Gutknecht, A. J., Makkeh, A., and Wibral, M.  From babel to boole: The logical organization of information decompositions.  2023\.
  * Ho et al. (2020) Ho, J., Jain, A., and Abbeel, P.  Denoising diffusion probabilistic models.  In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H. (eds.), _Advances in Neural Information Processing Systems_ , volume 33, pp. 6840–6851. Curran Associates, Inc., 2020.
  * Huang et al. (2021) Huang, C.-W., Lim, J. H., and Courville, A. C.  A variational perspective on diffusion-based generative models and score matching.  _Advances in Neural Information Processing Systems_ , 34:22863–22876, 2021.
  * Kingma & Ba (2015) Kingma, D. and Ba, J.  Adam: A method for stochastic optimization.  In _International Conference on Learning Representations (ICLR)_ , 2015.
  * Kolchinsky (2019) Kolchinsky, A.  A novel approach to multivariate redundancy and synergy.  _CoRR_ , abs/1908.08642, 2019.  URL http://arxiv.org/abs/1908.08642.
  * Kong et al. (2023) Kong, X., Liu, O., Li, H., Yogatama, D., and Steeg, G. V.  Interpretable diffusion via information decomposition.  2023\.
  * Latham & Nirenberg (2005) Latham, P. E. and Nirenberg, S.  Synergy, redundancy, and independence in population codes, revisited.  _Journal of Neuroscience_ , 25(21):5195–5206, 2005.
  * MacKay (2003) MacKay, D. J.  _Information theory, inference and learning algorithms_.  Cambridge university press, 2003.
  * Makkeh et al. (2021) Makkeh, A., Gutknecht, A. J., and Wibral, M.  Introducing a differentiable measure of pointwise shared information.  _Physical Review E_ , 103(3):032149, 2021.
  * Martinez Mediano (2022) Martinez Mediano, P. A.  Integrated information theory in complex neural systems.  2022\.
  * Nguyen et al. (2007) Nguyen, X., Wainwright, M. J., and Jordan, M.  Estimating divergence functionals and the likelihood ratio by penalized convex risk minimization.  In _Advances in Neural Information Processing Systems_ , 2007.
  * Oord et al. (2018) Oord, A. v. d., Li, Y., and Vinyals, O.  Representation learning with contrastive predictive coding.  _Advances in neural information processing systems_ , 2018.
  * Rosas et al. (2019) Rosas, F. E., Mediano, P. A. M., Gastpar, M., and Jensen, H. J.  Quantifying high-order interdependencies via multivariate extensions of the mutual information.  _Physical review. E_ , 100 3-1:032305, 2019.  URL https://api.semanticscholar.org/CorpusID:67855406.
  * Rosas et al. (2020) Rosas, F. E., Mediano, P. A. M., Rassouli, B., and Barrett, A.  An operational information decomposition via synergistic disclosure.  _Journal of Physics A: Mathematical and Theoretical_ , 53, 2020.  URL https://api.semanticscholar.org/CorpusID:210932609.
  * Runge et al. (2019) Runge, J., Bathiany, S., Bollt, E., Camps-Valls, G., Coumou, D., Deyle, E., Glymour, C., Kretschmer, M., Mahecha, M. D., Muñoz-Marí, J., et al.  Inferring causation from time series in earth system sciences.  _Nature communications_ , 10(1):2553, 2019.
  * Scagliarini et al. (2021) Scagliarini, T., Marinazzo, D., Guo, Y., Stramaglia, S., and Rosas, F. E.  Quantifying high-order interdependencies on individual patterns via the local o-information: Theory and applications to music analysis.  _Physical Review Research_ , 2021.  URL https://api.semanticscholar.org/CorpusID:237303787.
  * Scagliarini et al. (2023) Scagliarini, T., Nuzzi, D., Antonacci, Y., Faes, L., Rosas, F., Marinazzo, D., and Stramaglia, S.  Gradients of o-information: Low-order descriptors of high-order dependencies.  _Physical Review Research_ , 5, 01 2023.  doi: 10.1103/PhysRevResearch.5.013025.
  * Shannon (1948) Shannon, C. E.  A mathematical theory of communication.  _The Bell system technical journal_ , 27(3):379–423, 1948.
  * Song & Ermon (2019) Song, Y. and Ermon, S.  Generative modeling by estimating gradients of the data distribution.  In Wallach, H., Larochelle, H., Beygelzimer, A., d’ Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.
  * Song & Ermon (2020) Song, Y. and Ermon, S.  Improved techniques for training score-based generative models.  In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H. (eds.), _Advances in Neural Information Processing Systems_ , volume 33, pp. 12438–12448. Curran Associates, Inc., 2020.
  * Song et al. (2021) Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., and Poole, B.  Score-based generative modeling through stochastic differential equations.  In _International Conference on Learning Representations_ , 2021.
  * Sparacino et al. (2023) Sparacino, L., Faes, L., Mijatović, G., Parla, G., Re, V. L., Miraglia, R., de Ville de Goyet, J., and Sparacia, G.  Statistical approaches to identify pairwise and high-order brain functional connectivity signatures on a single-subject basis.  _Life_ , 13, 2023.  URL https://api.semanticscholar.org/CorpusID:264314627.
  * Stramaglia et al. (2021) Stramaglia, S., Scagliarini, T., Daniels, B. C., and Marinazzo, D.  Quantifying dynamical high-order interdependencies from the o-information: An application to neural spiking dynamics.  _Frontiers in Physiology_ , 11, 2021.  ISSN 1664-042X.  doi: 10.3389/fphys.2020.595736.  URL https://www.frontiersin.org/articles/10.3389/fphys.2020.595736.
  * Sun (1975) Sun, T.  Linear dependence structure of the entropy space.  _Inf Control_ , 29(4):337–68, 1975.
  * Sun Han (1980) Sun Han, T.  Multiple mutual informations and multiple interactions in frequency data.  _Information and Control_ , 46(1):26–45, 1980.  ISSN 0019-9958.  doi: https://doi.org/10.1016/S0019-9958(80)90478-7.  URL https://www.sciencedirect.com/science/article/pii/S0019995880904787.
  * Tax et al. (2017) Tax, T. M., Mediano, P. A., and Shanahan, M.  The partial information decomposition of generative neural network models.  _Entropy_ , 19(9), 2017.  ISSN 1099-4300.  doi: 10.3390/e19090474.  URL https://www.mdpi.com/1099-4300/19/9/474.
  * van Enk (2023) van Enk, S. J.  Pooling probability distributions and partial information decomposition.  _Physical review. E_ , 107 5-1:054133, 2023.  URL https://api.semanticscholar.org/CorpusID:256615444.
  * Varley et al. (2022) Varley, T. F., Pope, M., Faskowitz, J., and Sporns, O.  Multivariate information theory uncovers synergistic subsystems of the human cerebral cortex.  _Communications Biology_ , 6, 2022.  URL https://api.semanticscholar.org/CorpusID:249642639.
  * Varley et al. (2023) Varley, T. F., Pope, M., Puxeddu, M. G., Faskowitz, J., and Sporns, O.  Partial entropy decomposition reveals higher-order information structures in human brain activity.  _Proceedings of the National Academy of Sciences of the United States of America_ , 120, 2023.  URL https://api.semanticscholar.org/CorpusID:255825886.
  * Venkatesh et al. (2023) Venkatesh, P., Bennett, C., Gale, S., Ramirez, T. K., Heller, G., Durand, S., Olsen, S. R., and Mihalas, S.  Gaussian partial information decomposition: Bias correction and application to high-dimensional data.  In _Thirty-seventh Conference on Neural Information Processing Systems_ , 2023.  URL https://openreview.net/forum?id=1PnSOKQKvq.
  * Villani (2009) Villani, C.  _Optimal transport: old and new_ , volume 338.  Springer, 2009.
  * Vincent (2011) Vincent, P.  A connection between score matching and denoising autoencoders.  _Neural Computation_ , 23(7):1661–1674, 2011.
  * Williams & Beer (2010) Williams, P. L. and Beer, R. D.  Nonnegative decomposition of multivariate information, 2010.

## Score-based o-information Estimation — Supplementary material

##  Appendix A Proofs

We here prove the equivalences about tc and dtc. Starting from tc

$\displaystyle\sum_{i=1}^{N}\mathcal{H}(X^{i})-\mathcal{H}(X)=\sum_{i=1}^{N}\mathcal{H}(X^{i})-\sum_{i=1}^{N}\mathcal{H}(X^{i}\,|\,X^{>i})=\sum\limits_{i=1}^{N-1}\mathcal{I}(X^{i};X^{>i})=\mathcal{T}(X)$

Concerning dtc

$\displaystyle\mathcal{H}(X)-\sum_{i=1}^{N}\mathcal{H}(X^{i}|X^{\setminus
i})=\mathcal{H}(X^{1})+\mathcal{H}(X^{\setminus
1}\,|\,X^{1})-\mathcal{H}(X^{1}|X^{\setminus
1})-\sum_{i=2}^{N}\mathcal{H}(X^{i}|X^{\setminus i})=$
$\displaystyle\mathcal{I}(X^{1};X^{\setminus 1})+\mathcal{H}(X^{\setminus
1}\,|\,X^{1})-\sum_{i=2}^{N}\mathcal{H}(X^{i}|X^{\setminus i})=$
$\displaystyle\mathcal{I}(X^{1};X^{\setminus
1})+\mathcal{H}(X^{2}\,|\,X^{1})+\mathcal{H}(X^{\setminus
1,2}\,|\,X^{1},X^{2})-\mathcal{H}(X^{2}|X^{\setminus
2})-\sum_{i=3}^{N}\mathcal{H}(X^{i}|X^{\setminus i})=$
$\displaystyle\mathcal{I}(X^{1};X^{\setminus
1})+\mathcal{I}(X^{2};X^{>2}|X^{1})+\mathcal{H}(X^{\setminus
1,2}\,|\,X^{1},X^{2})-\sum_{i=3}^{N}\mathcal{H}(X^{i}|X^{\setminus i})=\dots$
$\displaystyle\sum\limits_{i=1}^{N-1}\mathcal{I}(X^{i};X^{>i}\,|\,X^{<i})=\mathcal{D}(X)$

where for the last equality it suffices to consider trivial reordering
arguments,
$\sum\limits_{i=2}^{N}\mathcal{I}(X^{i};X^{<i}\,|\,X^{>i})=\sum\limits_{i=1}^{N-1}\mathcal{I}(X^{i};X^{>i}\,|\,X^{<i})$.

##  Appendix B Details of s$\Omega$i

In the section we provide additional implementation details about s$\Omega$i.

###  B.1 Computing o-information

In § 3.1, we presented how tc and dtc can be estimated using denoising score
functions. Our estimators requires different score functions which can be
obtained by learning different denoisers. More particularly, tc requires the
joint denoiser $\mathbb{E}[X\,|\,X_{t}]$ and the marginals
$\mathbb{E}[X^{i}\,|\,X^{i}_{t}]$ for $i\in\\{1,\dots,N\\}$. dtc estimation is
obtained using the joint and the following conditional terms
$\mathbb{E}[X^{i}\,|\,X^{i}_{t},X^{\setminus i}]$ for $i\in\\{1,\dots,N\\}$.
Our formulation in § 3.1 is general and can be applied to a wide range of
denoising score learning techniques. For the implementation of s$\Omega$i, we
adopt VP-Stochastic Differential Equation (sde) framework (Song & Ermon,
2019). The latter perturbs the data using an sde parameterized by a drift
$f_{t}$ and a diffusion coefficient $g_{t}$.

#### Muti-variate denoising score network.

We extend the work from (Bounoua et al., 2023) to amortize the learning of all
the required terms using a single denoising score network. The denoising score
network $\epsilon_{\theta}$ accepts as input the concatenation of the
variables each perturbed at different times. The second input is a vector of
size $N$ which describes the state of each variable and allows a
parametrization of different denoising score functions.

The joint term corresponds to the case where all the variables are perturbed
with the same intensity $t$ and all the elements of the vector
$\tau=[t,\dots,t]$ are set equivalently to $t$. The conditional terms
correspond to the case where only the conditioned variable $i$ is perturbed
with intensity $t$ whereas the remaining conditioning variables $\setminus
i_{\text{th}}$ are kept unperturbed at $t=0$. Consequently the parameter
describing this case is of the form $[0,\dots,t,\dots,0]$.

While (Bounoua et al., 2023) framework is not able to learn the marginal
denoising score, it’s possible via an additional parameterization to include
this configuration. This corresponds to the case where the marginal variable
$i$ is perturbed with intensity $t$ while all the other variables are made
uninformative. The non marginal variables $\setminus i_{\text{th}}$ are
replaced with pure noise corresponding to a maximal perturbation at $t=T$.
Consequently the parameter describing this case is of the form
$[T,\dots,t,\dots,T]$.

#### Training.

The training is carried out through a randomized procedure. At each training
step, we select randomly a set of the denoising score functions required for
the o-information estimation (joint, conditional or marginals). These
denoising scores function are learned by the single network following
Algorithm 1. In total, estimating o-information requires calling $2N+1$
denoising score functions which we learn using a single denoising network.

1

Data: $X=\\{X^{i}\\}_{i=1}^{N}$

$t\sim\mathcal{U}[0,T]$

// Importance sampling schemes (Huang et al., 2021; Song et al., 2021) can be
adopted to reduce variance

if  _Joint_ then

$X_{t}\sim p_{t}$

// Obtain noisy version of all the variables using vpsde (Song & Ermon, 2019)
with drift $f_{t}$ and diffusion coefficient $g_{t}$.

$s_{t}(X_{t})=\epsilon_{\theta}([X_{t}^{1},\dots,X_{t}^{N}],\tau=[t,\dots,t,\dots,t])$

Return $\nabla_{\theta}\norm{s_{t}(X_{t})-\nabla\log{p_{t}(X_{t}|X)}}$

// Denoising score matching of all the variables

if  _Conditional_ then

$X^{i}_{t}\sim p_{t}$

// Obtain noisy version of the variable $i$ while the remaining variables are
kept unperturbed at ($t=0$)

$s_{t}(X^{i}_{t}|X^{\setminus
i})=\epsilon_{\theta}\left([X^{1},\dots,X^{i-1},X^{i}_{t},X^{i+1},\dots,X^{N}],\tau=[0,\dots,t,\dots,0]\right)$

Return $\nabla_{\theta}\norm{s_{t}(X^{i}_{t}|X^{\setminus
i})-\nabla\log{p_{t}(X^{i}_{t}|X^{i})}}$

// Denoising score matching of the conditioning variable $i$

if  _Marginal_ then

$X^{i}_{t}\sim p_{t}$

$X^{\setminus i}_{T}\leftarrow p_{T}=\mathcal{N}(0,\mathbb{I})$

// Obtain noisy version of the variable $i$ while the remaining variables are
replaced with pure noise ($t=T$).

$s_{t}(X^{i}_{t})=\epsilon_{\theta}([X_{T}^{1},\dots,X_{T}^{i-1},X_{t}^{i},X_{T}^{i+1},\dots,X_{T}^{N}],\tau=[T,\dots,t,\dots,T])$

Return
$\nabla_{\theta}\norm{s_{t}(X^{i}_{t})-\nabla\log{p_{t}(X^{i}_{t}|X^{i})}}$

// Denoising score matching of the marginal variable $i$

Algorithm 1 s$\Omega$i Training step

#### Inference.

Once all the denoising score functions are learned, it’s possible to estimate
tc and dtc via a Monte Carlo estimation of the integral over $t$ in
Proposition 2 and Proposition 4 . The outer integration w.r.t. to the time
instant is possible by sampling $t\sim\mathcal{U}(0,T)$, and then using the
estimation $\int_{0}^{T}(\cdot)\differential
t=T\mathbb{E}_{t\sim\mathcal{U}(0,T)}[(\cdot)]$. In practice we adopt 10 steps
for the computation of the expectation. The procedure to estimate
o-information is described in algorithm 2. First, samples from $x\sim p(x)$
are considered, then sampling the time $t\sim\mathcal{U}[0,T]$. A perturbed
version of the variables $X_{t}$ is computed using the Variance preserving SDE
(vpsde). The joint, conditional and marginal denoising scores are computed
leveraging the single denoising score network. This is possible by choosing
different perturbation times and manipulating the vector $\tau$ as described
earlier. Computing the difference of the denoising scores functions (see
Proposition 2 and Proposition 2 ) allows the computation of tc and dtc
respectively. Please note that it is possible to implement importance sampling
schemes to reduce the variance, along the lines of what described by Huang et
al. (2021).

Data: $X=\\{X^{i}\\}_{i=1}^{N}$

$t\sim\mathcal{U}[0,T]$

// Importance sampling scheme can also be adopted

$X_{t}\sim p_{t}$

// Obtain the noisy version of all the variables using vpsde (Song & Ermon,
2019) with drift $f_{t}$ and diffusion coefficient $g_{t}$.

$s_{t}(X_{t})\leftarrow\epsilon_{\theta}([X_{t}^{1},\dots,X_{t}^{N}],\tau=[t,\dots,t,\dots,t])$
// Compute the joint score

for  _i = 1 to N // Compute the conditional and marginal terms_

__ do

$s_{t}(X^{i}_{t}|X^{\setminus
i})\leftarrow\epsilon_{\theta}\left([X^{1},\dots,X^{i-1},X^{i}_{t},X^{i+1},\dots,X^{N}],\tau=[0,\dots,t,\dots,0]\right)$

$s_{t}(X^{i}_{t})\leftarrow\epsilon_{\theta}([X_{T}^{1},\dots,X_{T}^{i-1},X_{t}^{i},X_{T}^{i+1},\dots,X_{T}^{N}],\tau=[T,\dots,t,\dots,T])$

// Similarly to Algorithm 1 the non marginal variables are replaced with pure
noise $X_{T}^{\setminus i}\sim\mathcal{N}(0,\mathbb{I})$

end for

$\hat{\mathcal{T}}(X)\leftarrow\frac{g^{2}_{t}}{2}\norm{s_{t}(X_{t})-\left[s_{t}(X^{i}_{t})\right]_{i=1}^{N}}^{2}$
// See Proposition 2

$\hat{\mathcal{D}}(X)\leftarrow\frac{g^{2}_{t}}{2}\norm{s_{t}(X_{t})-\left[s_{t}(X^{i}_{t}|X^{\setminus
i})\right]_{i=1}^{N}}^{2}$ // See Proposition 4

$\hat{\Omega}(X)\leftarrow\hat{\mathcal{T}}(X)-\hat{\mathcal{D}}(X)$

Return $\hat{\Omega}(X)$

Algorithm 2 s$\Omega$i inference time

###  B.2 Computing gradient of o-information

To compute the gradient of o-information recall that
$\partial_{i}\Omega(X)=\Omega(X)-\Omega(X^{\setminus i})$. The first order
gradient of o-information requires the estimation of o-information of all the
subsystems of size $N-1$.

$\displaystyle\Omega(X^{\setminus i})=\mathcal{T}(X^{\setminus
i})-\mathcal{T}(X^{\setminus i})$ (7) $\displaystyle=\sum_{j=1,j\neq
i}^{N}\mathcal{H}(X^{j})-\mathcal{H}(X^{\setminus i})$ (8)
$\displaystyle-(\mathcal{H}(X^{\setminus i})-\sum_{j=1,j\neq
i}^{N}\mathcal{H}(X^{j}|X^{\setminus\\{i,j\\}}))$ (9)

Many denoising score functions in Equation 9 were also used to estimate the
global o-information. To learn the additional necessary terms to compute
$\Omega(X^{\setminus i})$, the randomized set of scores adopted during the
training step (see § B.1) is extended to account for the new requirements.
Please note that we still use a single denoising network that considers all
the terms necessary to compute o-information and its gradient. A large number
of learned denoising score functions is a potential reason for the bias
observed in our experiment Figure 4. A highly flexible architecture capable of
fitting large number of scores may be needed to infer gradient of
o-information.

##  Appendix C Experimental settings

###  C.1 Canonical multivariate Gaussian system

In this section we provide additional details about the construction of the
synthetic benchmark § 4.1.

#### Redundancy benchmark.

All the variable of the system are composed of a redundant component and
unique information specific to each variable.

We modulate the redundant inter-dependency strength by setting different
values for $\sigma$. We consider a standardized system where all the variables
mean is $0$ and standard deviation equal to $\mathbb{I}$. This results in the
following covariance matrix:

$\begin{bmatrix}\mathbb{I}&\rho\mathbb{I}&\vdots&\rho\mathbb{I}\\\
\rho\mathbb{I}&\mathbb{I}&\dots&\rho\mathbb{I}\\\
\vdots&\vdots&\ddots&\rho\mathbb{I}\\\
\rho\mathbb{I}&\rho\mathbb{I}&\dots&\mathbb{I}\\\ \end{bmatrix}$ (10)

With $\rho=\frac{1}{1+\sigma^{2}}$ which modulates the interactions strength
in the system.

#### Synergy benchmark.

We consider a standardized system where all the variables mean is $0$ and
standard deviation equal to $\mathbb{I}$. This results in the following
covariance matrix :

$\begin{bmatrix}\mathbb{I}&\frac{1}{\sqrt{N-1}}\mathbb{I}&0&\dots&0\\\
\frac{1}{\sqrt{N-1}}\mathbb{I}&\mathbb{I}&\frac{\rho}{\sqrt{N-1}}\mathbb{I}&\dots&\frac{\rho}{\sqrt{N-1}}\mathbb{I}\\\
0&\frac{\rho}{\sqrt{N-1}}\mathbb{I}&\mathbb{I}&\dots&0\\\
0&\vdots&0&\ddots&0\\\ 0&\frac{\rho}{\sqrt{N-1}}&0&\dots&\mathbb{I}\\\
\end{bmatrix}$ (11)

Where $\rho=\frac{1}{\sqrt{1+\sigma^{2}}}$ modulates the interactions strength
in the system.

#### Mixed benchmark.

The covariance matrix is easy to obtain as the mixed benchmark is made of
independent subsystems.

#### Ground Truth.

Having access to the covariance matrix of the system, computing entropy in
close form for Gaussian distribution is possible. For
$X\sim\mathcal{N}(\mu,\sigma)$ :

$\mathcal{H}(X)=\frac{1}{2}\log(2\pi\sigma^{2})+\frac{1}{2}\ $ (12)

For a multivariate Gaussian distribution
$X^{d}\sim\mathcal{N}_{d}(\mu,\Sigma)$ :

$\mathcal{H}(X)=\frac{D}{2}(1+\log(2\pi))+\frac{1}{2}\log det(\Sigma)\ $ (13)

###  C.2 s$\Omega$i implementation details

We provide code-base for s$\Omega$i implementation at
111https://github.com/MustaphaBounoua/soi . The training of s$\Omega$i is
carried out using Adam optimizer (Kingma & Ba, 2015). We use Exponential
moving average (EMA) with a momentum parameter $m=0.999$. Importance sampling
(Huang et al., 2021) (222https://github.com/CW-Huang/sdeflow-light) at train
and test-time. The hyper-parameters are presented in Table 1. To estimate the
gradient of o-information (Figure 4) the model width is double the one
presented in Table 1 to account for the additional necessary terms to learn.
Concerning the experiments in Figure 5 , we use the same architecture used for
the canonical examples and follow the same procedure to choose the model
capacity( see Table 1 for the hyper-parameters details).

Table 1: s$\Omega$i network training details. $Dim$ of the task correspond the sum of the dimensions of all variables of the system. For the neural data application we report the number of training iterations (.,.) corresponding the ”change” case and ”No change” case. The number of iteration used for the ”No change” is higher since the dataset contains more ”no change” flashes compared to ”change” flashes. | Width | Time embed | Batch size | Lr | Iterations | Number of params
---|---|---|---|---|---|---
($Dim\leq 50$) | 128 | 128 | 256 | 1e-2 | 195k | 320k
($Dim\leq 100$) | 192 | 192 | 256 | 1e-2 | 195k | 747k
($Dim\geq 100$) | 256 | 256 | 256 | 1e-2 | 195k | 1003k
Neural application |  |  |  |  |  |
($Dim\leq 30$) | 128 | 128 | 256 | 1e-2 | (100k,160k) | 320k
($Dim\leq 75$) | 192 | 192 | 256 | 1e-2 | (100k,160k) | 737k
($Dim\leq 150$) | 256 | 256 | 256 | 1e-2 | (100k,160k) | 1300k
($Dim\geq 150$) | 384 | 384 | 256 | 1e-2 | (100k,160k) | 3000k

###  C.3 Baselines

(Bai et al., 2023) decomposes tc into $N-1$ mi terms which are estimated using
pairwise neural mi estimator. Similarly by leveraging Equation 16 dtc can also
be retrieved by estimating $N-1$ additional mi terms.

$\displaystyle\mathcal{T}(X)=\sum\limits_{i=1}^{N-1}\mathcal{I}(X^{i};X^{>i})$
(14)
$\displaystyle\mathcal{D}(X)=\mathcal{S}(X)-\mathcal{T}(X)=\sum\limits_{i=1}^{N}\mathcal{I}(X^{i};X^{\setminus
i})-\mathcal{T}(X)$ (15)
$\displaystyle\mathcal{D}(X)=\sum\limits_{i=2}^{N}\mathcal{I}(X^{i};X^{\setminus
i})-\sum\limits_{i=2}^{N-1}\mathcal{I}(X^{i};X^{>i})$ (16)

Our implementation in based on the the official codebase
333https://github.com/Linear95/TC-estimation of (Bai et al., 2023). We use the
same architecture and hyper parameters from (Bai et al., 2023):
$\text{LR}=1e-3$, $\text{Batch size}=64$. We use an mlp architecture for all
the variants of the baseline with 3 linear layers with varying width. For each
mi term, the capacity of the neural network is aligned to the input dimension.
Adam optimizer (Kingma & Ba, 2015) is used for training. We increase the width
of the hidden layer to accommodate the data dimension. For the variant of the
baseline implemented with mine, we used smaller layer size as large capacity
led to divergence during training. To ensure the best performance, we train
each mi estimator model for $80k$ steps for a number of variables $N=10$ and
$40k$ for number of variables $N=6$. In the different experiments, we reported
the performance results averaged over 5 seeds and dropped the baseline in case
of divergence during training.

###  C.4 The Visual Behavior Neuropixels

Hereafter we describe the different pre-processing steps applied on the Visual
Behavior Neuropixels in § 4.2. We follow the same procedure described by
(Venkatesh et al., 2023). The selected mice are the ones with both familiar
and novel sessions and a minimum number of 20 units in each of the six brain
regions: VISp, VISl, VISal ,VISrl, VISam and VISpm. Only the units of good
quality are kept. The selection criteria was based on an SNR at least 1, and
with fewer than 1 inter-spike interval violations. The non-change flashes
correspond to the ones where the image does not change and happen between 4
and 10 flashes after the trial start. Trials corresponding to a change are
naturally the ones when the image has changed. Only flashes that occurred
while the animal was engaged ( based on the reward information) is kept, while
the ones corresponding to an omission, or after an omission, and flashes
during which the animal licked, were all removed.

The trials were aligned to the start of each stimulus flash, and the 250ms
recordings were divided into 5 bins of 50 ms duration averaged over the units
of the same region. We use different step sizes to count the spikes which
resulted in different dimensional representation but resulted in the same
intuition (See Figure 11,Figure 10 and Figure 9). Please note that unlike
(Venkatesh et al., 2023), we don’t use PCA to reduce the dimension of the
data, and count the number of spikes per unit by averaging the activity over
the units of the same region indexed by time.

##  Appendix D Additional results

(a) Dim=5

(b) Dim=10

(c) Dim=15

(d) Dim=20

Figure 6:  Redundant system with 6 variables, organized into subsets of sizes
$\\{3,3\\}$ and increasing interaction strength.

(a) Dim=5

(b) Dim=10

(c) Dim=15

(d) Dim=20

Figure 7:  Synergistic system with 6 variables, organized into subsets of
sizes $\\{3,3\\}$ and increasing interaction strength.

(a) Dim=5

(b) Dim=10

(c) Dim=15

(d) Dim=20

Figure 8:  Mixed-interaction system with 6 variables, organized into a
redundancy-dominant subsets of size $3$ variables and one synergy-dominant
subset with $3$ variables. o-information is modulated by fixing the synergy
inter-dependency and increasing the redundancy.

(a) 3 areas

(b) 6 areas

(c) 3 areas

(d) 6 areas

Figure 9: o-information and s-information estimate in the visual cortex region
activity after two types of stimulus flash across 72 trial sessions. Left:
Analysis using three brain region areas, Right: Extended analysis using six
brain region areas. The step size is set to $1ms$ which results in 50
dimensional data for each bin per area.

(a) 3 areas

(b) 6 areas

Figure 10: s-information estimate in the visual cortex region activity after
two types of stimulus flash across 72 trial sessions. Left: Analysis using
three brain region areas, Right: Extended analysis using six brain region
areas. The step size is set to $2ms$ which results in 25 dimensional data for
each bin per area.

(a) 3 areas

(b) 6 areas

(c) 3 areas

(d) 6 areas

Figure 11: o-information and s-information estimate in the visual cortex
region activity after two types of stimulus flash across 72 trial sessions.
Left: Analysis using three brain region areas, Right: Extended analysis using
six brain region areas. The step size is set to $5ms$ which results in 10
dimensional data for each bin per area.