# Neural dSCA: demixing multimodal interaction among brain areas during
naturalistic experiments

Yu Takagi
University of Tokyo, Japan
<EMAIL_ADDRESS>
Laurence T. Hunt
University of Oxford, UK
<EMAIL_ADDRESS>
Ryu Ohata
University of Tokyo, Japan
<EMAIL_ADDRESS>
Hiroshi Imamizu
University of Tokyo/ATR, Japan
<EMAIL_ADDRESS>
Jun-ichiro Hirayama
AIST, Japan
<EMAIL_ADDRESS>
Corresponding author.

###### Abstract

Multi-regional interaction among neuronal populations underlies the brain’s
processing of rich sensory information in our daily lives. Recent neuroscience
and neuroimaging studies have increasingly used naturalistic stimuli and
experimental design to identify such realistic sensory computation in the
brain. However, existing methods for cross-areal interaction analysis with
dimensionality reduction, such as reduced-rank regression and canonical
correlation analysis, have limited applicability and interpretability in
naturalistic settings because they usually do not appropriately ‘demix’ neural
interactions into those associated with different types of task parameters or
stimulus features (e.g., visual or audio). In this paper, we develop a new
method for cross-areal interaction analysis that uses the rich task or
stimulus parameters to reveal how and what types of information are shared by
different neural populations. The proposed neural demixed shared component
analysis combines existing dimensionality reduction methods with a practical
neural network implementation of functional analysis of variance with latent
variables, thereby efficiently demixing nonlinear effects of continuous and
multimodal stimuli. We also propose a simplifying alternative under the
assumptions of linear effects and unimodal stimuli. To demonstrate our
methods, we analyzed two human neuroimaging datasets of participants watching
naturalistic videos of movies and dance movements. The results demonstrate
that our methods provide new insights into multi-regional interaction in the
brain during naturalistic sensory inputs, which cannot be captured by
conventional techniques.

##  1 Introduction

Naturalistic stimuli have been increasingly used in experimental neuroscience
in both human [1, 2] and animal studies [3] to investigate richer and more
realistic stimulus-induced neural activities than those in traditional
experimental design (Figure 1a). Recent advancements in data-driven feature
extraction methods, including deep neural network [4, 5], also contribute to
this trend by drastically reducing the annotation costs for complex stimuli.
One of the biggest challenges of using naturalistic stimuli is that they are
usually non-categorical and high-dimensional, which is in stark contrast to
traditional stimuli that are typically designed to be categorical and low-
dimensional.

In addition to the increasing complexity of experimental stimuli, in recent
neuroscience studies, there has been an increasing emphasis on investigating
complex interactions between neural populations across brain areas [6, 7, 8,
9, 10] (Figure 1b) rather than simple stimulus-response relationships of
individual neuron (or voxel) groups. This is important because the brain is a
network system that consists of a tremendous number of neurons that interact
with each other. The ongoing revolution in both measurement techniques and the
scale of data provides neuroscientists with a fascinating opportunity to
advance the understanding of cross-areal interactions in both neurophysiology
[10] and human neuroimaging [6, 7]. In particular, several authors have
emphasized the importance of analyzing ‘multivariate’ (multi-to-multi)
interactions between brain areas beyond ‘univariate’ interactions identified
among individual neurons/voxels or their simple averages across areas. Several
statistical techniques for dimensionality reduction have been used [11] to
explore multivariate cross-areal interactions in an interpretable, and
computationally and statistically efficient manner.

Figure 1: Basic ideas and properties of the proposed NdSCA. a In traditional
experiments (left), task parameters are usually pre-defined, categorical, and
low-dimensional, whereas in naturalistic experiments (right), an arbitrary
number of task parameters can, in principle, be extracted from continuous,
high-dimensional, and usually multimodal (e.g., images or audio) naturalistic
stimuli. b Comparison among different types of interactions. Areas A, B, and C
are communicating in a univariate (top) or multivariate (bottom) manner. c
Mixed (top) and demixed (bottom) interactions in matrix format. Coloured areas
indicate significant interaction. Only the lower-triangular part is visualized
as we primarily assume symmetric interactions

However, the use of naturalistic stimuli complicates the analysis of cross-
areal interactions, which poses several technical challenges. Existing methods
that use canonical correlation analysis (CCA) [10] or reduced-rank regression
(RRR) do not take external variables into account, which results in cross-
areal connections that are ‘mixed’ in terms of their association with task
parameters, such as stimulus features and other contextual factors (Figure
1c). The uncontrolled, high-dimensional, and multi-modal nature of
naturalistic stimuli further makes the scenario extremely difficult because
previous attempts in task-related neural interaction analysis [12, 13, 14, 15,
16] commonly assumed traditional experimental designs and thus suffered from
conceptual or computational problems. In fact, how to efficiently ‘demix’
cross-areal interactions into those associated with different task parameters
of naturalistic experiments is an open problem. Novel technical advances are
strongly desired so that one can fully take advantage of large-scale datasets
and the rich information in naturalistic stimuli.

In this paper, we propose neural demixed shared component analysis
(NdSCA)111All code will be published after acceptance., as well as a
simplified alternative, as a novel technique for cross-areal neural
interaction analysis using dimensionality reduction. The method extends
demixed shared component analysis (dSCA) [17], proposed previously for
traditional experimental design, in a non-trivial manner (see Table 1 for
comparison with other methods). In particular, we propose using a recently
proposed neural-network implementation of functional analysis of variance
(ANOVA) with latent variables [18] to efficiently demix the complicated
nonlinear effects of naturalistic stimuli into those associated with different
types of task parameters. The continuous and multimodal nature of naturalistic
stimuli cannot be handled by the conventional ANOVA on which dSCA relies. Our
NdSCA combines this novel demixing technique with linear dimensionality
reduction so that interactions among brain regions and their association with
various task parameters are detected and visualized in an interpretable, and
both statistically and computationally efficient manner.

This paper is organized as follows. First, we briefly review previous methods
for cross-areal interaction analysis, in addition to methods related to the
idea of demixing (Section 2), and contrast their availability in traditional
and naturalistic scenarios in neuroscience experiments. Then, we introduce the
underlying ideas and specific algorithms of NdSCA and its simplified version
(Section 3). Finally, we validate our methods using both simulation and
applications to two human functional magnetic resonance imaging (fMRI)
datasets (Section 4), and then conclude with a brief discussion (Section 5).

|  | FC | PPI | RRR | dSCA | NdSCA-simple | NdSCA
---|---|---|---|---|---|---|---
Interaction among areas | demixed | - | ✓ | - | ✓ | ✓ | ✓
multivariate | - | - | ✓ | ✓ | ✓ | ✓
Covariate (task parameter) | multiple | - | ✓ | - | - | ✓ | ✓
multimodal | - | ✓ | - | ✓ | - | ✓
continuous | - | ✓ | - | - | ✓ | ✓
nonlinear | - | - | - | - | - | ✓
Table 1: Comparison among the methods for detecting interaction among areas.

##  2 Related work

#### Detecting information sharing across different brain areas.

In previous studies, information sharing among different brain areas was
investigated for different scales from pairs of neurons [19] to interaction
among multi-voxels [16, 20]. Particularly in human neuroimaging, a simple
univariate correlation between large brain regions of interest (ROIs) is often
calculated, which is called functional connectivity (FC). In line with the
increased scale of datasets, both in the field of neurophysiology and human
neuroimaging, researchers have started to investigate multivariate information
sharing. Because simple regression often does not work in a high-dimensional
setting, they have used statistical techniques such as principal component
analysis (PCA) [8, 16], RRR [9], or CCA [10] to effectively reduce the
dimensionality. However, these methods do not provide results that are demixed
with respect to experimental task parameters. This is problematic because,
although the results of several studies have suggested that task-related
information sharing might have richer information than that of the resting
state [21], the content of the shared information is unclear.

#### Demixing task parameters

It is well known that even a single neuron has mixed selectivity, i.e., the
neuron’s firing rate responds to more than one task parameter [22]. This is
true even after dimensionality reduction [23]. Mixed selectivity is also
commonly observed in a large-scale voxel in human fMRI [1]. It is a critical
problem for researchers who are interested in cognitive processing during
complex experiments because the analysis and interpretation then become
complicated in order to dissociate different types of computations run
simultaneously in the brain. To overcome this problem, several approaches have
been proposed. Among others, demixed PCA (dPCA) [14, 15] has successfully been
used in neurophysiology, which was recently extended to dSCA to capture
information sharing among areas [17]. The key idea is to ‘marginalize’
multivariate information in a single area in association with a specific task
parameter of interest. These methods allow experimenters to combine rich data
with equally rich task design, identifying low-dimensional components that
vary along axes defined by features of the experimental task. In human
neuroimaging, psychophysiological interaction (PPI) [12] and beta series
correlation (BSC) [13], which are closely related mathematically [24], have
specifically been used for a long time but they are designed for traditional
experiments. In fact, PPI and BSC primary consider univariate interactions and
also not applicable to large-dimensional task parameters with a reasonable
computational cost.

##  3 Neural demixed shared component analysis (NdSCA)

###  3.1 Cross-areal interaction analysis with dimensionality reduction

Our goal is to investigate task or stimulus-related interactions among brain
areas induced by naturalistic stimuli in an efficient and interpretable
manner. Assume that brain activity is measured simultaneously in two areas.
For each area, we observe the multivariate time series of firing rates or
blood oxygen level-dependent signals sampled at T discrete time points. We
then provide a conventional FC measure using the correlation of average
activities between two areas, but this is suboptimal if the interaction is not
through averages of the areas.

Linear dimensionality reduction techniques are particularly useful for
evaluating the multivariate interaction of neuron/voxel groups between brain
areas. A popular method in the neuroscience literature is RRR. Let
$\mathbf{X}$ and $\mathbf{Y}$ denote given data matrices of two areas, with
sizes $M_{X}\times T$ and $M_{Y}\times T$, respectively, where $M_{X}$ and
$M_{Y}$ denote the number of neurons or voxels in respective areas and T
denotes the time length. Then, standard least-squares RRR minimizes the
following:

$L_{\mathrm{RRR}}=\left\|\mathbf{Y}-\mathbf{WX}\right\|^{2},$ (1)

where $\|.\|$ represents the Frobenius norm; $\mathbf{W}$ is a coefficient
matrix of size $M_{X}\times M_{Y}$ and its rank is constrained to not be
greater than $P$ $(<\min(M_{X},M_{Y}))$. The solution $\mathbf{W}_{RRR}$ is
obtained as

$\mathbf{W}_{RRR}=\mathbf{W}_{OLS}\mathbf{VV}^{T},$ (2)

where $\mathbf{W}_{OLS}$ is the ordinary least-squares solution and the
columns of the $M_{X}\times P$ matrix $\mathbf{V}$ contain the top $P$
principal components of the optimal linear predictor
$\mathbf{\hat{Y}}_{OLS}:=\mathbf{W}_{OLS}\mathbf{X}$ [15]. Note that RRR
reduces to CCA if $\mathbf{Y}$ is whitened [25]. Given an RRR solution, we
define the interaction strength between the two areas as the total explained
variance (EV) of the target activity.

Although RRR might achieve higher sensitivity than simple correlation, RRR
still fails to capture task-related information sharing, that is, results are
mixed in terms of task parameters. Therefore, in the following, we consider
task-related information sharing explicitly.

###  3.2 Demixing effects of naturalistic stimuli via functional ANOVA

Suppose that we are given multiple task parameters of interest, denoted by
$c_{k}$ ($k=1,2,\ldots,K$), where $K$ is the number of different types of task
parameters, in addition to brain activity measurements. First, we consider
traditional experiments such that $c_{k}$ are discrete (univariate) and K is
relatively small. Then, their induced effects on brain activity may be modeled
using the well-known framework of (multivariate) ANOVA, which underlies the
theory of dPCA and its extension for cross-area interaction analysis, that is,
dSCA. The key idea is that ANOVA eventually decomposes a given data matrix
$\mathbf{Y}$ into multiple matrices, each corresponding to the main or
interaction effects by the discrete task parameters; for instance, if $K=2$,
two-way ANOVA identifies the decomposition, given by

$\mathbf{Y}=\mathbf{Y}_{1}+\mathbf{Y}_{2}+\mathbf{Y}_{12}+\mathbf{E},$ (3)

where $\mathbf{Y}_{0}$ corresponds to the overall mean; $\mathbf{Y}_{1}$ and
$\mathbf{Y}_{2}$ correspond to the main effects of $c_{1}$ and $c_{2}$,
respectively; $\mathbf{Y}_{12}$ corresponds to the interaction effects; and
$\mathbf{E}$ corresponds to residuals that are not explained by task
parameters.222Note that in each matrix, except for $\mathbf{E}$, the columns
are identical if their associated task-parameter values are the same; for
example, $\mathbf{Y}_{0}$ replicates the overall mean vector across columns,
and the columns in $\mathbf{Y}_{1}$ with the same value of $c_{1}$ are
identical, each given by the corresponding $c_{1}$-conditional empirical mean.
The procedure to obtain the ANOVA decomposition in addition to each resultant
term was called marginalization in [15, 17]. Given the decomposition, dSCA
between source $\mathbf{X}$ and target $\mathbf{Y}$ is formalized as solving
RRR for each marginalization of $\mathbf{Y}$, with the $\mathbf{Y}$ in Eq.(3)
replaced by the corresponding term in the ANOVA decomposition; dSCA reduces to
dPCA if the source and target areas are the same. Intuitively, each RRR then
identifies a rank-reduced connection from the source to target areas such that
the connection specifically explains the target variability related to the
particular subset of task parameters.

Generalizing the idea of dSCA in naturalistic settings is, however, not
straightforward. The problem is that the idea of marginalization strongly
relies on conventional ANOVA, and thus cannot directly be applicable if a task
parameter $c_{k}$ is continuous. Moreover, features obtained from naturalistic
stimuli are usually rather high-dimensional, even in each individual modality,
such as visual or audio. In practice, modeling the effect of high-dimensional
task parameters easily leads to high-computational complexity if not combined
with any suitable parametric architecture and associated algorithms. Another
potential issue is that all the task parameters of interest would not
necessarily be observed (extracted), which is more likely to occur in
naturalistic experiments than traditional experiments as the stimuli are not
strongly controlled.

To overcome all these issues, we propose extending the basic idea of dSCA
using neural decomposition (ND) [18], which is a recently proposed practical
instantiation of functional ANOVA [26] with a conditional variational
autoencoder (CVAE) [27]. Functional ANOVA generalizes the marginalization of
dSCA in a principled manner, and the use of the deep neural network technique
enhances the practical utility in high-dimensional settings. We refer to the
resultant analysis framework as NdSCA. In practice, the algorithm simply
replaces the marginalization part of dSCA with the generalized part using ND
(Algorithm 1). Therefore, in the following, we focus on presenting the basic
idea of ND and its interpretation from a naturalistic experiment viewpoint.

Algorithm 1 Neural demixed shared component analysis (NdSCA) between two brain
areas

0: Matrices $(\mathbf{X},\mathbf{Y})$ of source and target brain activity;
those of $K$ task parameter vectors $\mathbf{c}_{k}$

0: Rank-reduced coefficient matrices $\mathbf{W}_{RRR}$ for all relevant
marginalizations

1: Run ND [18] to decompose $\mathbf{Y}$ into marginalizations
$\mathbf{Y}_{\phi}$ corresponding to a functional ANOVA

2: for each relevant marginalization $\mathbf{Y}_{\phi}$ (e.g.
$\phi\in\\{z,c,zc\\}$ for the model (4)) do

3: Solve an RRR problem (1) from $\mathbf{X}$ to $\mathbf{Y}_{\phi}$ with
Eq.(2)

4: end for

Now, we consider a naturalistic setting in which each task parameter, denoted
by ${\bf c}_{k}$, may be continuous and vector-valued. For example, ${\bf
c}_{1}$ and ${\bf c}_{2}$ represent visual and audio features of the
naturalistic stimulus extracted at each time frame of the brain activity
measured. To model the potential effects from unobserved but relevant task
parameters, which are more likely to exist in naturalistic settings than
traditional settings, we further introduce an additional random vector z of
latent variables. To ease exposition, assume for simplicity that the number of
observed task-parameter vectors is $K=1$. Then, we can generally write the
task-induced effect on brain activity $\mathbf{y}$ as $F({\bf z,c})$, where
$F$ denotes any appropriate nonlinear mapping. The idea of functional ANOVA
and thus ND is to decompose $F({\bf z,c})$ into marginal (main) and
interaction effects of ${\bf c}$ and ${\bf z}$, such that

$\mathit{F}(\bf{z,c})=\mathit{F}_{\mathrm{0}}+\mathit{F}_{\mathrm{z}}(z)+\mathit{F}_{\mathrm{c}}(c)+\mathit{F}_{\mathrm{zc}}(z,c)$
(4)

where $F_{0}$ denotes the intercept and other terms; $F_{c}(\mathbf{c})$ and
$F_{z}(\mathbf{z})$ represent the main effects of ${\bf c}$ and ${\bf z}$,
respectively; and $F_{\mathrm{z}c}(\bf z,c)$ represents their interaction
effect. Generalization with more than one task-parameter vector
$\mathbf{c}_{k}$ is straightforward, with their main and interaction effects
appropriately included up to a reasonable order. Note that a decomposition
such as Eq.(4) is clearly non-unique if not further constrained; functional
ANOVA ensures the identifiability of the model by introducing constraints on
each term such that its integration over every input argument is zero [18],
which naturally generalizes the typical sum-to-zero constraint in classical
ANOVA.

Although the non-parametric form of functional ANOVA above is theoretically
relevant, it is not practically useful and thus needs to be further
implemented with parametric architectures, which must easily incorporate
possibly high-dimensional $\bf{c}_{k}$ and latent input $\bf z$ in our
context. From this viewpoint, the recent neural network-based realization, ND,
of functional ANOVA is particularly suitable as it offers flexible and
computationally efficient implementation based on a CVAE [28], explicitly
modeling and performing inference on latent variables $\bf z$. Briefly, a CVAE
learns a parametric generative model $p_{\theta}(\bf{y}|\bf{z,c})$ defined by
the decoder network $\bf{y}=F_{\theta}(\bf{z,c})$ and a prior $p(\bf{z})$,
while approximating the posterior $p(\bf{z}|\bf{y,c})$ parametrically with an
encoder network $q_{\phi}(\bf{z}|\bf{y})$. Both the parameters $\phi$ and
$\theta$ are trained by approximately maximizing the variational lower bound
of the conditional log-likelihood (i.e., $\log p_{\theta}(\bf{y}|\bf{c})$ for
a single datum), typically using stochastic gradient ascent. In ND, the
additive form and specific input dependences of each term in Eq.(4) are
explicitly modeled in the decoder network, and penalty and augmented
Lagrangian methods are developed to deal with the integral constraints (see
[18] for details).

Given the decoder network learned, we can obtain an additive matrix
decomposition similar to Eq.(4) by first estimating latent vector $\bf z$ for
every column of the target matrix $\mathbf{Y}$ and then computing terms such
as $F_{c}(\mathbf{z,c})$ at each instance of $\bf z$ and $\bf c$. Then, if we
perform RRR as above by regarding each of these terms as a generalized
marginalization of $\mathbf{Y}$, NdSCA obtains rank-reduced cross-areal
interaction and low-dimensional components of activities in the areas. We
quantify the interaction so that it is sensitive only to a particular (main or
interaction) effect of observed and unobserved task parameters. Note that the
final RRR step of NdSCA provides a linear relationship among voxels/neurons
and thus the contribution of each voxel/neuron can be easily visualized and
interpreted. This interpretability is crucial for neuroscientists; previous
approaches based on decoding [20] or representational similarity analysis [16]
did not provide it.

###  3.3 Simplified alternative

The novel marginalization procedure of NdSCA can be applied to various
scenarios of neuroscientific data analysis, even in non-naturalistic settings,
because of the flexibility of ND and the underlying functional ANOVA model.
However, the method can be overly complicated if we consider more restricted
and rather controlled (more traditional) scenarios in which the effect of non-
observed task parameters (latent variables $\bf z$) can be ignored. In
particular, a simple alternative would be sufficient if only a single type
(modality) of continuous stimuli was used and the linearity of stimulus
effects could be reasonably assumed.

Specifically, we consider a single multivariate task parameter, which is
continuous and possibly high-dimensional. Let $\mathbf{C}$ denote the
corresponding data matrix, with the instances of the task parameter vectors as
its columns. Then, if no latent variables are necessarily introduced, the
linear effect on a single-region activity is modeled simply by
$\mathbf{Y}=\mathbf{WC}+\mathbf{E}$ without resorting to more general
functional ANOVA modeling. We can therefore use multivariate linear regression
or similar techniques to identify the task-specific effect in $\bf{Y}$, that
is, $\bf{Y}_{\bf{C}}:=\bf{WC}$; in this case, we specifically use the same
RRR/CCA technique as above to improve the identification.

In fact, this simple alternative, which we refer to as NdSCA-simple, is novel
and useful in its own right. This also offers a reasonable baseline for
validating the performance of NdSCA in our simulation and real-data analysis
below. Note, however, that the method is suboptimal if the simplifying
assumptions are violated, for example, when we need to simultaneously
incorporate multiple types (modalities) of task parameter vectors, as is
typical in naturalistic experiments.

##  4 Results

#### Synthetic data

To confirm that demixing is necessary for capturing task-parameter-specific
information sharing and NdSCA can deal with naturalistic experiments for that
purpose, we generated simulated neuronal populations that is analogous to the
naturalistic experiment.

Suppose that we simultaneously recorded population neural activities in three
brain areas, A, B, and C, during a traditional experiment with two task
parameters of interest: Audio and Image (see Figure 1a right for an
illustration). Both parameters are high-dimensional and continuous, such as
features extracted using an artificial neural network. We set the dimension of
the task parameters and neurons to 10 and 40, respectively. We synthesized the
data so that the interactions (co-activations) of neural activity between
areas A, B, and C varied in response to Audio and those between A and C varied
in response to Image, as indicated by the two colors in Figure 2a. We
simulated $3\times 10^{4}$ observations, and split them into training,
validation, and test data (1000/1000/1000). See the Supplementary Material for
the details.

For comparison, we applied FC, RRR, PPI, and our NdSCA to the dataset. For
each method, we calculated the EVs of predictions from one region to another.
To calculate the strength of interactions, for example, between areas A and B,
we averaged the EVs of predictions from A to B and B to A. We calculated the
statistical significance of interactions by running permutation testing 1,000
times. For PPI and NdSCA, we also calculated the statistical significance for
each area in terms of the representation of the specific task parameter using
ND, and calculated the strength of interactions only if both source and target
areas had information about the task parameter. To train ND, we used the Adam
optimizer [27]. We trained the model at a maximum of 50 epochs with a batch
size of 128. We used early stopping with a patience of 10 epochs. We set the
number of latent units $\bf z$ to 1. We used a warm-up process that weighted
the reconstruction loss more than the Kullback–Leibler loss [29]. For FC and
PPI, we applied PCA to reduce the dimensions of neurons to one, whereas for
RRR and NdSCA, we used multi-dimensional data without PCA. We also applied PCA
to the task parameters to reduce the dimensions to three before applying PPI
and NdSCA.

Figure 2: Simulations demonstrate that only NdSCA detected complex task-
parameter specific interaction in the naturalistic experiment. a Areas A, B,
and C communicate Audio information (top; indicated by red lines), whereas
areas A and C communicate Image information (bottom; indicated by blue line).
b Results obtained from standard FC, RRR, PPI, and NdSCA ordered from left to
right. Only NdSCA detected demixed interaction. n.s. corresponds to a no
significant interaction.

We found that FC and PPI could not detect any information sharing because they
do not use an underlying multivariate structure (Figure 2b; $P>0.05$).
Although RRR detected information sharing among areas ($P<0.05$ for lines in
Figure 2b), it did not indicate what information was shared across areas.
Strikingly, only NdSCA revealed the demixed underlying information sharing
structure ($P<0.05$ for coloured lines in Figure 2b). This is because only
NdSCA could manage multivariate information with complex task parameters.

#### Naturalistic experiment – movie watching fMRI

Next, we applied NdSCA to a human fMRI dataset of participants watching movies
[30] to validate the advantage of our methods. The dataset consisted of 86
participants watching 10 full-length movies. We used six participants watching
‘The Shawshank Redemption’ that lasted for 8,181 s (Figure 3a). The authors
also provided time-aligned Word (transcribed from audio) and Face (appear/not
appear) annotations for each movie. We used a binary feature vector for Face
and three-dimensional latent semantic matrices for Word that were converted
using fastText ([31]; see the Supplementary Material). We split all scans into
training, validation, and test data (with 4,800, 1,200 and 2,181 scans,
respectively).

We investigated information sharing among small-scale ROIs and also among
large-scale functional networks. Specifically, we used 10 ROIs obtained from
automated anatomical labeling ([32]; Figure 3b, see Supplementary Material),
which are related to faces, words and associative functions, and eight
functional networks in [33] (Figure 3e). Before applying ND, we applied PCA to
each ROI/network to reduce the dimensionality to either 20 (as source) or 50
(as taarget). NdSCA was run based on a functional ANOVA model for the two
types of task parameters (i.e. face and word) and one latent vector
$\mathbf{z}$, while their interaction terms were not included to simplify the
analysis. We calculated EVs for each participant, and averaged across
participants for visualization and statistical testing.

Figures 3d and 3g show that RRR detected information sharing among almost all
ROIs and networks ($P<0.05$ for coloured pairs in Figures 3d and 3g; FWE
corrected). However, as shown in the simulation analysis, they did not
indicate what information they shared. While PPI detected task-related
interaction to some extent, e.g., Face-related interaction between fusiform
face area (FFA) and V1, PPI failed to detect Word-related information sharing.
Strikingly, only NdSCA successfully revealed Word- and Face-related
information sharing. The ROIs and networks that were known for visual
information processing (e.g., FFA for the ROIs; Visual 1 and Visual 2
networks) exhibited Face-related information sharing with each other, whereas
areas known for audio and semantics (e.g., A1 and temporal areas for the ROIs;
medial temporal/frontoparietal networks) Word-related information sharing
($P<0.05$ for coloured pairs in Figures 3d and 3g; FWE corrected). Notably,
NdSCA had much higher sensitivity than NdSCA-simple, which might suggest that
the underlying relationship between brain activity and the task parameters is
nonlinear.

We also investigated the demixed interaction associated with the effect of
latent variables $\mathbf{z}$ (representing non-observed but relevant task or
stimulus parameters) by applying NdSCA with the $F_{z}(\mathbf{z})$ term in
Eq.(4) reconstructed for every estimated instance of z. The result appeared to
be the same as that obtained by non-demixed RRR (Supplementary Figure 1),
although we did not investigate their differences in detail as this is beyond
the scope of this paper. The apparent similarity, at least, indicates that a
large number of variations in cross-regional interaction can be nonlinearly
explained by a relatively small number of factors that are not related to Face
or Word. This is in fact reasonable because naturalistic movies contain many
other features and also cross-areal interaction is strongly affected by the
intrinsic interaction, reflecting spontaneous brain activity, even during
tasks [34, 35].

Figure 3: NdSCA successfully demixed different types of computation. a
Schematic of the task. Human participants watched 2-hours of a movie. We
extracted a binary face vector (appear/not appear) and multi-dimensional
latent semantic matrices as features. b We used 10 ROIs. See the Supplementary
Material. c ROIs were divided into visual, audio/temporal, and others. d
Results obtained from RRR, PPI, NdSCA-simple, and NdSCA are ordered from left
to right. For PPI and NdSCAs, the results of Face (top) and Word (bottom) are
displayed separately. e-g We also conducted similar analysis for eight
functional networks. Non-significant pairs are coloured black. n.s.
corresponds to a no significant interaction.

#### Naturalistic experiment – Dance clip fMRI

Finally, we applied NdSCA to a human fMRI dataset of participants watching
dance clips (Figure 4a). Participants watched 1,163 dance clips that were
performed by 30 dancers from 10 genres for 10–50 seconds [36, 37], which
resulted in nearly five hours and approximately 15,000 scans (see the
Supplementary Material).333We collected the fMRI data used in this experiment.
We split all scans into training, validation, and test data (with 10,874,
1,200 and 2,756 scans). Note that the dance clips in the test dataset
consisted of music and choreography that did not appear in both training and
validation datasets. We used the same ROIs and networks as in the previous
analysis.

We used automatically extracted audio and motion features from the clips. For
audio features, we used the publicly available audio processing toolbox
Librosa [38] to extract a one-dimensional envelope, mel-frequency cepstral
coefficients, and chroma, resulting in a 33-dimensional feature vector. For
the motion features, we used 17 COCO-format three-dimensional (3D) human joint
locations, which resulted in a 51-dimensional feature vector [37]. We applied
PCA to reduce the dimensionality of each feature to three. Again, NdSCA with
the two type of feature (task parameter vectors) and one latent vector was
used with no interaction terms among them for simplicity. See the
Supplementary Material for the details.

Figures 4c and 4e show that, again, only NdSCA successfully detected task-
parameter-specific information sharing ($P<0.05$ for coloured pairs in Figures
4c and 4e; FWE corrected). It shows that audio-related information sharing was
propagated across the brain, whereas motion-related information sharing was
more strongly observed around the visual ROIs/networks.

Figure 4: NdSCA demixed information sharing into meaningful task parameters,
whereas the conventional method failed. a Schematic of the task. Human
participants watched 5-hours of dance clips. We extracted multi-dimensional 3D
joint locations (motion) and audio features. b We used the same ROIs as in the
previous analysis. c Results obtained from RRR, PPI, NdSCA-simple, and NdSCA
are ordered from left to right. d-e Same analyses for eight functional
networks. Non-significant pairs are coloured black. n.s. corresponds to a no
significant interaction.

Given that only NdSCA provided task-specific information sharing, we further
investigated whether NdSCA could also reveal the different information sharing
patterns between different conditions. Although previous studies investigated
different patterns of interaction between different conditions [7], there
measures are mixed in terms of task-parameters, thus make it difficult to
interpret the results. In this dataset, there was a control condition under
which participants watched dance clips played backward (Fig 5a). We
hypothesized that pattern of interaction among areas in the backward-play
condition was substantially changed from those in the forward-play condition.
We also hypothesized that the effects were not equal across different task
parameters, audio and motion.

Figure 5: NdSCA reveals strong double-dissociation of effects of backward
playing between audio- and motion-related information shraing a Comparison
between forward and backward play. b-c NdSCA (right columns) indicated that
interactions for audio and motion were differently manipulated by backward
playing, whereas RRR could not demonstrate this (left column) for the ROIs
($\bf b$) and networks ($\bf c$). Non-significant pairs are coloured white.

Figures 5b and 5c show that NdSCA confirmed our hypothesis: manipulation
effects on audio and motion-related information sharing were drastically
different, rather than demonstrating a mere overall decrease/increase of
information sharing ($P<0.05$ for coloured pairs in Figures 5b and 5c; FWE
corrected). Specifically, audio-related information sharing changed globally,
whereas motion-related information sharing did not change consistently across
the ROIs/networks. Interestingly, audio-related interaction around A1 did not
suffer even in the condition of backward playing (indicated by the black
arrow). All the results above could never be observed by conventional
connectivity methods, including RRR, which suggests that our proposed method
had an advantage over existing methods [7, 21].

##  5 Conclusion and discussion

Here we proposed NdSCA, a new technique for demixing complex task parameter
specific information sharing among brain areas during naturalistic
experiments. With both simulations and naturalistic fMRI data analyses, we
demonstrated that NdSCA can provide a novel neuroscience insight by demixing
patterns of cross-areal information sharing into those associated with
different types of task or stimulus parameters (features). Previously, non-
demixed types of techniques have commonly been used in cross-areal interaction
analysis and thus often yielded results that were not easily interpreted.

However, our method has several limitations. First, in return for an
expression power, NdSCA needs relatively larger data compared to simple linear
methods. However, given that the size of the publicly available dataset is
getting larger, we believe that NdSCA can contribute to the trend. Second, the
current implementation of NdSCA only captures linear information sharing among
areas. Although it provides interpretability, future study should consider
non-linear information sharing among areas.

Future research topics include: (i) applying NdSCA to time-resolved data
obtained by, e.g. microelectrodes or magnetoencephalography; (ii)
investigating the behavioural relevance of the shared component; (iii)
extending NdSCA to deal with more than two areas, that can be easily
implemented thanks to the flexible architecture of NdSCA.

## Acknowledgments and Disclosure of Funding

YT was supported by JSPS overseas fellowship. LH was supported by a Sir Henry
Dale Fellowship from Wellcome and the Royal Society (208789/Z/17/Z). HI was
supported by JSPS KAKENHI (19H05725 and 21H00959) and AMED (JP18dm0307008). JH
was supported by a project, JPNP20006, commissioned by the New Energy and
Industrial Technology Development Organization (NEDO), and by JSPS KAKENHI
(18KK0284 and 21K12055).

## References

  * [1] Alexander G. Huth, Shinji Nishimoto, An T. Vu, and Jack L. Gallant.  A Continuous Semantic Space Describes the Representation of Thousands of Object and Action Categories across the Human Brain.  Neuron, 76(6):1210–1224, 2012.
  * [2] Uri Hasson, Yuval Nir, Ifat Levy, Galit Fuhrmann, and Rafael Malach.  Intersubject Synchronization of Cortical Activity during Natural Vision.  Science, 303(5664):1634–1640, 2004.
  * [3] Carsen Stringer, Marius Pachitariu, Nicholas Steinmetz, Matteo Carandini, and Kenneth D. Harris.  High-dimensional geometry of population responses in visual cortex.  Nature, 571(7765):361–365, 2019.
  * [4] Umut Güçlü and Marcel A.J. van Gerven.  Deep neural networks reveal a gradient in the complexity of neural representations across the ventral stream.  Journal of Neuroscience, 35(27):10005–10014, 2015.
  * [5] Alexander Mathis, Pranav Mamidanna, Kevin M. Cury, Taiga Abe, Venkatesh N. Murthy, Mackenzie Weygandt Mathis, and Matthias Bethge.  DeepLabCut: markerless pose estimation of user-defined body parts with deep learning.  Nature Neuroscience, 21(9):1281–1289, 2018.
  * [6] David C Van Essen, Stephen M Smith, Deanna M Barch, Timothy EJ Behrens, Essa Yacoub, Kamil Ugurbil, Wu-Minn HCP Consortium, et al.  The wu-minn human connectome project: an overview.  Neuroimage, 80:62–79, 2013.
  * [7] Erez Simony, Christopher J Honey, Janice Chen, Olga Lositsky, Yaara Yeshurun, Ami Wiesel, and Uri Hasson.  Dynamic reconfiguration of the default mode network during narrative comprehension.  Nature communications, 7(1):1–13, 2016.
  * [8] Matthew T. Kaufman, Mark M. Churchland, Stephen I. Ryu, and Krishna V. Shenoy.  Cortical activity in the null space: Permitting preparation without movement.  Nature Neuroscience, 17(3):440–448, 2014.
  * [9] João D. Semedo, Amin Zandvakili, Christian K. Machens, Byron M. Yu, and Adam Kohn.  Cortical Areas Interact through a Communication Subspace.  Neuron, 102(1):249–259.e4, 2019.
  * [10] Nicholas A. Steinmetz, Peter Zatka-Haas, Matteo Carandini, and Kenneth D. Harris.  Distributed coding of choice, action and engagement across the mouse brain.  Nature, 576(7786):266–273, 2019.
  * [11] John P. Cunningham and Byron M. Yu.  Dimensionality reduction for large-scale neural recordings.  Nature Neuroscience, 17(11):1500–1509, 2014.
  * [12] K J Friston, C Buechel, G R Fink, J Morris, E Rolls, and Raymond J Dolan.  Psychophysiological and modulatory interactions in neuroimaging.  Neuroimage, 6(3):218–229, 1997.
  * [13] Jesse Rissman, Adam Gazzaley, and Mark D’Esposito.  Measuring functional connectivity during distinct stages of a cognitive task.  NeuroImage, 23(2):752–763, 2004.
  * [14] Wieland Brendel, Ranulfo Romo, and Christian K. MacHens.  Demixed principal component analysis.  Advances in Neural Information Processing Systems 24: 25th Annual Conference on Neural Information Processing Systems 2011, NIPS 2011, pages 1–9, 2011.
  * [15] Dmitry Kobak, Wieland Brendel, Christos Constantinidis, Claudia E. Feierstein, Adam Kepecs, Zachary F. Mainen, Xue Lian Qi, Ranulfo Romo, Naoshige Uchida, and Christian K. Machens.  Demixed principal component analysis of neural population data.  eLife, 5(APRIL2016):1–36, 2016.
  * [16] Takuya Ito, Kaustubh R. Kulkarni, Douglas H. Schultz, Ravi D. Mill, Richard H. Chen, Levi I. Solomyak, and Michael W. Cole.  Cognitive task information is transferred between brain regions via resting-state network topology.  Nature Communications, 8(1):1–13, 2017.
  * [17] Yu Takagi, Steven Kennerley, Jun-ichiro Hirayama, and Laurence Hunt.  Demixed shared component analysis of neural population data from multiple brain areas.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 6235–6244. Curran Associates, Inc., 2020.
  * [18] Kaspar Märtens and Christopher Yau.  Neural Decomposition: Functional ANOVA with Variational Autoencoders.  In International Conference on Artificial Intelligence and Statistics, pages 2917–2927. PMLR, 2020.
  * [19] L. G. Nowak, M. H.J. Munk, A. C. James, P. Girard, and J. Bullier.  Cross-correlation study of the temporal interactions between areas V1 and V2 of the macaque monkey.  Journal of Neurophysiology, 81(3):1057–1074, 1999.
  * [20] Marc N. Coutanche and Sharon L. Thompson-Schil.  Informational connectivity: Identifying synchronized discriminability of multi-voxel patterns across the brain.  Frontiers in Human Neuroscience, 7(JAN):1–14, 2013.
  * [21] Emily S Finn and Peter A Bandettini.  Movie-watching outperforms rest for functional connectivity-based prediction of behavior.  NeuroImage, 235:117963, 2021.
  * [22] Mattia Rigotti, Omri Barak, Melissa R. Warden, Xiao Jing Wang, Nathaniel D. Daw, Earl K. Miller, and Stefano Fusi.  The importance of mixed selectivity in complex cognitive tasks.  Nature, 497(7451):585–590, 2013.
  * [23] Valerio Mante, David Sussillo, Krishna V. Shenoy, and William T. Newsome.  Context-dependent computation by recurrent dynamics in prefrontal cortex.  Nature, 503(7474):78–84, 2013.
  * [24] Xin Di, Zhiguo Zhang, and Bharat B. Biswal.  Understanding psychophysiological interaction and its relations to beta series correlation.  Brain Imaging and Behavior, 15(2):958–973, 2021.
  * [25] Fernando De Torre.  A Least-Squares Framework for Component Analysis.  Ieee Transactions on Pattern Analysis and Machine Intelligence, 34(6):1041–1055, 2012.
  * [26] Ilya M Sobol.  Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates.  Mathematics and computers in simulation, 55(1-3):271–280, 2001\.
  * [27] Diederik P. Kingma and Max Welling.  Auto-encoding variational bayes.  2nd International Conference on Learning Representations, ICLR 2014 - Conference Track Proceedings, (Ml):1–14, 2014.
  * [28] Kihyuk Sohn, Xinchen Yan, and Honglak Lee.  Learning structured output representation using deep conditional generative models.  Advances in Neural Information Processing Systems, 2015-January:3483–3491, 2015.
  * [29] Lars Maaløe, Casper Kaae Sønderby, Søren Kaae Sønderby, and Ole Winther.  Auxiliary Deep Generative Models.  In Maria Florina Balcan and Kilian Q Weinberger, editors, Proceedings of The 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pages 1445–1453, New York, New York, USA, 2016. PMLR.
  * [30] Sarah Aliko, Jiawen Huang, Florin Gheorghiu, Stefanie Meliss, and Jeremy I. Skipper.  A naturalistic neuroimaging database for understanding the brain using ecological stimuli.  Scientific Data, 7(1):1–21, 2020.
  * [31] Piotr Bojanowski, Edouard Grave, Armand Joulin, and Tomas Mikolov.  Enriching Word Vectors with Subword Information.  Transactions of the Association for Computational Linguistics, 5:135–146, 2017.
  * [32] N Tzourio-Mazoyer, B Landeau, D Papathanassiou, F Crivello, O Etard, N Delcroix, B Mazoyer, and M Joliot.  Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain.  NeuroImage, 15(1):273–89, jan 2002.
  * [33] Emily S Finn, Xilin Shen, Dustin Scheinost, Monica D Rosenberg, Jessica Huang, Marvin M Chun, Xenophon Papademetris, and R Todd Constable.  Functional connectome fingerprinting: identifying individuals using patterns of brain connectivity.  Nature Neuroscience, 18(October):1–11, 2015.
  * [34] I. Tavor, O Parker Jones, R B Mars, S M Smith, T E Behrens, and S Jbabdi.  Task-free MRI predicts individual differences in brain activity during task performance.  Science (New York, N.Y.), 352(6282):216–20, 2016.
  * [35] Michael W. Cole, Takuya Ito, Danielle S. Bassett, and Douglas H. Schultz.  Activity flow over resting-state networks shapes cognitive task activations.  Nature Neuroscience, 19(12):1718–1726, 2016.
  * [36] Shuhei Tsuchida, Satoru Fukayama, Masahiro Hamasaki, and Masataka Goto.  AIST Dance Video Database: Multi-Genre, Multi-Dancer, and Multi-Camera Database for Dance Information Processing.  In ISMIR, pages 501–510, 2019.
  * [37] Ruilong Li, Shan Yang, David A Ross, and Angjoo Kanazawa.  Learn to Dance with AIST++: Music Conditioned 3D Dance Generation.  arXiv preprint arXiv:2101.08779, 2021.
  * [38] Brian McFee, Colin Raffel, Dawen Liang, Daniel PW Ellis, Matt McVicar, Eric Battenberg, and Oriol Nieto.  librosa: Audio and music signal analysis in python.  In Proceedings of the 14th python in science conference, volume 8, pages 18–25. Citeseer, 2015.