# Optimal Transport for Latent Integration with An Application to
Heterogeneous Neuronal Activity Data

Yubai Yuan  Department of Statistics, The Pennsylvania State University  Babak
Shahbaba  Department of Statistics, University of California Irvine  Norbert
Fortin  Department of Neurobiology and Behavior, University of California
Irvine  Keiland Cooper  Department of Neurobiology and Behavior, University of
California Irvine  Qing Nie  Department of Mathematics, Department of
Developmental and Cell Biology, Department of Biomedical Engineering,
University of California Irvine  Annie Qu 111Corresponding author:
<EMAIL_ADDRESS>Department of Statistics, University of California Irvine

###### Abstract

Detecting dynamic patterns of task-specific responses shared across
heterogeneous datasets is an essential and challenging problem in many
scientific applications in medical science and neuroscience. In our motivating
example of rodent electrophysiological data, identifying the dynamical
patterns in neuronal activity associated with ongoing cognitive demands and
behavior is key to uncovering the neural mechanisms of memory. One of the
greatest challenges in investigating a cross-subject biological process is
that the systematic heterogeneity across individuals could significantly
undermine the power of existing machine learning methods to identify the
underlying biological dynamics. In addition, many technically challenging
neurobiological experiments are conducted on only a handful of subjects where
rich longitudinal data are available for each subject. The low sample sizes of
such experiments could further reduce the power to detect common dynamic
patterns among subjects. In this paper, we propose a novel heterogeneous data
integration framework based on optimal transport to extract shared patterns in
complex biological processes. The key advantages of the proposed method are
that it can increase discriminating power in identifying common patterns by
reducing heterogeneity unrelated to the signal by aligning the extracted
latent spatiotemporal information across subjects. Our approach is effective
even with a small number of subjects, and does not require auxiliary matching
information for the alignment. In particular, our method can align
longitudinal data across heterogeneous subjects in a common latent space to
capture the dynamics of shared patterns while utilizing temporal dependency
within subjects. Our numerical studies on both simulation settings and
neuronal activity data indicate that the proposed data integration approach
improves prediction accuracy compared to existing machine learning methods.

Key words: Distribution alignment, Electrophysiological study, Gromov-
Wasserstein barycenter, Latent embedding, Optimal transport, Semi-supervised
learning

##  1 Introduction

New emerging and rapidly expanding technologies for conducting high throughput
scientific studies have allowed the simultaneous recording of activity from
large ensembles of neurons during complex behavior [3, 16]. These technical
advances are expected to lead to paradigm-shifting breakthroughs in our
fundamental understanding of information processing in the brain, as well as
to a novel platform to investigate sources of deficits and potential
treatments in neurological disorders. However, at this time, the main limiting
factor toward achieving this goal is that existing statistical and
computational methods are insufficient to address the new challenges
associated with extracting and integrating the essential information contained
in these large and complex multidimensional datasets. Among these challenges,
by far the most difficult task is how to properly integrate data across
multiple subjects from a highly heterogeneous population [17, 26]. This paper
provides a fundamental solution and a rigorous framework to address this
issue.

Our proposed method is motivated by, and applied to, a recent
electrophysiological neuroscience study, which was designed to elucidate the
neuronal mechanisms underlying memory for sequences of events, an ability
essential to our normal daily life function. Although it is well established
that this ability depends on the integrity of the hippocampus, a region near
the center of our brain, it remains unclear how the activity of ensembles of
hippocampal neurons can trigger forms of memory. Since the structure and
function of the hippocampus is well observed across mammals, this type of
information is primarily obtained by recording neural activity in the
hippocampus of rodents due to the invasiveness of the approach. In particular,
considerable rodent research has shown that hippocampal neurons can encode
sequences of locations, but there has been no direct evidence that the coding
property can be extended to the non-spatial memories which are typically
investigated in human studies.

The experiment recorded neural activity in the CA1 region of the hippocampus
as rats performed a nonspatial sequential memory task such as testing a
sequence of 5 odors (see Figure 1) where neuron functional activity was
recorded by neuronal electrophysiological data [29]. Although statistical
machine learning methods can be used to extract and visualize latent
representations of sequential neuronal activities associated with different
tasks in each subject, it is still a great challenge to integrate subject-
specific latent representations across different subjects. Yet this data
integration is a crucial step, as the primary objective of such studies is to
identify fundamental coding properties which are common across subjects [2,
28].

Figure 1: Neural activity was recorded from hippocampal region CA1 as animals
performed a complex non-spatial sequence memory task. (a) The task involves
repeated presentations of sequences of non-spatial events (odor stimuli) to
subjects (rodents). Using an automated odor delivery system (top left),
sequences of five odors were presented in the same odor port. (b) Example
ensemble activity from representative subject during one sequence
presentation.

To integrate cross-subject neuronal data, there are several neuronal data
integration methods designed to directly extract the cross-subject latent
spaces, i.e., model the subject-wise neuronal signal as the reconstruction
from a set of reference neuronal dynamics, which are estimated in the common
latent space via Gaussian process [15], canonical correlation analysis [37],
and principal component analysis [23]. These methods assume that the subject-
wise signals are linear transformations of the common latent dynamics, which
implies that the heterogeneity among subjects can be captured via linear
transformations between their neuronal data. This assumption can be
restrictive and oversimplified when the individual heterogeneity in neuronal
activities is complex and strong.

In this paper, we propose to aggregate information by finding the common
relations between neuron activity patterns and experimental outcomes among
different subjects. On the other hand, this type of data is known to contain
noise, redundancies, and subject-specific idiosyncrasies, which poses great
difficulties in developing methods to identify the underlying common neural
patterns among subjects. One of the most significant challenges is the
intrinsic heterogeneity within neuron activities [12, 24]. Specifically, the
underlying mapping between neuron activity pattern and a specific outcome can
be differentially expressed among different subjects due to systematic
differences across subjects caused by random deviation from common patterns.
Also, this kind of experiment is typically conducted on only a handful of
subjects where high-dimensional neuron signal data are available from each
subject [12, 24]. The systematic heterogeneity across individuals and low
sample size significantly undermines the capacity of traditional machine
learning methods to identify cross-subjects patterns. To alleviate this
problem, a new data integration method is highly needed to properly account
for heterogeneity in order to discern the underlying common patterns across
subjects. Specifically, we propose a novel integrated latent alignment (ILA)
framework by incorporating deep representation learning [7, 18, 6] and optimal
transport (OT) [33, 1] to integrate heterogeneous data.

Compared with existing data integration methods, optimal transport can relax
the linearity assumption to reconstruct the nonlinear transformation between
individual spike sequences and cross-subject spike patterns, and therefore one
can integrate heterogeneous spike data with more flexible distributions from
subjectwise perturbations. In addition, the optimal-transport-based barycenter
is able to integrate multiple distributions of neuronal encodings to obtain
shared geometric structures [27], which enables one to detect common
clustering patterns from heterogeneous neural encodings of different subjects’
spike data, and preserve common spike patterns via filtering out subjects’
heterogeneity when aligning spike data.

The proposed ILA framework consists of three stages: First, we use deep
learning to project the temporal neural data to a low-dimensional latent space
for each subject, which allows us to identify the macro neuron patterns from
each subject. Notice that the representations of the neuron activity are
projected into different latent spaces for different subjects, or the
compressed neural features from each subject. Therefore, in the second stage,
we develop a novel data integration method based on optimal transport [1, 14,
31]. Specifically, we incorporate cross-subject heterogeneity by aligning both
the geometric structure of the latent neural features and the feature-response
patterns from different subjects. The cross-subject alignment allows us to
aggregate subject-wise data in a shared latent space. In the last stage, we
adopt either supervised learning or unsupervised learning methods in the
aligned latent space to capture the common neural patterns across subjects.

The proposed ILA has several significant advantages. First, it provides a
principle strategy to integrate data from heterogeneous subjects and increase
statistical power in identifying the common patterns. In addition, the
integration is effective even when the number of subjects is small. Secondly,
the proposed method does not require restrictive assumptions on the mappings
for the latent feature spaces from heterogeneous subjects. Instead, the ILA
utilizes the smoothness in association between neural activity features and
responses to align the patterns directly, which allows for perturbations
between responses and features across subjects. Finally, the proposed method
can incorporate heterogeneous temporal latent features, and task-specific
machine learning algorithms. Both the numerical studies on simulations and the
electrophysiological data show that the proposed data integration method
improves the performance of classification tasks compared with methods using
individual data only.

This paper is organized as follows: Section 2 introduces the background of the
electrophysiological study and optimal transport. Section 3 introduces the
proposed data integration method for supervised learning and the corresponding
optimization algorithms. Section 4 demonstrates simulation studies, and
Section 5 illustrates an application of the proposed method to an
electrophysiological study. The last section provides conclusions and
discussion.

##  2 Background and Notations

In this section, we provide a brief introduction to our motivating data from
the electrophysiological study, and then introduce the optimal transport
techniques related to our method. The electrophysiological data were collected
as rats performed a non-spatial sequence memory task [4]. Briefly, the task
involves repeated presentations of a sequence of odors (odors A,B,C,D, and E)
in the same odor port, and requires subjects to determine whether each odor is
presented “in sequence” (InSeq; e.g., ABC…) or “out of sequence” (OutSeq;
e.g., ABD…) to receive a water reward (Figure 1a). Note that the task was
self-paced in that each odor presentation was initiated by the rat placing its
nose in the port (median interval between consecutive odors $\sim$5 s). In
each session, the same sequence was presented multiple times (Figure 1a,
right), with approximately half the presentations including all stimuli InSeq
(Figure 1a,top) and the other half including one stimulus OutSeq (Figure
1a,bottom). Incorrect InSeq/OutSeq judgments resulted in termination of the
sequence. Neural spiking activity was recorded using 22 tetrodes (bundles of
four electrodes) in rats tested on the well-trained sequence (ABCDE). For each
odor presentation (trial), the data typically features spike counts from
$\sim$50-90 neurons, and trial identifiers (e.g., odor presented,
InSeq/OutSeq, response correct/incorrect). Figure 1b illustrates examples of
ensemble activity (putative principal neurons and interneurons) from a
representative subject during one sequence presentation. More detailed
description of the data is in the Section 11 of Supplemental Materials.

In summary, this dataset is well-suited to develop and test the required
methodology to integrate neural activity data across subjects. First, the
dataset was shown to contain complex and distributed forms of coding that were
undetectable when examining the activity of individual neurons independently,
but emerged when considering the activity of each ensemble of neurons as a
whole. Second, the size of the dataset is in a tractable range (50-90 neurons
per subject); the ensembles are large enough to detect emergent coding
properties, but small enough to allow the iterative development and testing of
analytical and computational solutions. Third, the task structure provides
points of convergence across subjects (e.g., common sequences of stimuli and
associated neural states) as well as points of divergence (e.g., errors by a
specific animal).

Here, we use optimal transport (OT), which has been widely applied in machine
learning fields (e.g., domain adaptation and transfer learning), to provide a
principled way of aligning distributions while preserving the geometric
structures. The classical OT based on the Wasserstein distance requires that
data distributions be defined on the same space, which may not be satisfied in
many real applications where data could arise from different domains of
distributions. To overcome this limitation, the Gromov–Wasserstein distance
[25] and fused Gromov–Wasserstein distance [32] are introduced to extend the
concept of OT to cross-domain distribution alignment.

Gromov-Wasserstein distance  Consider two datasets
$\\{(X^{(s)},Y^{(s)}),P^{(s)}\\}\;(s=1,2)$ from two different metric-measure
spaces where $(X^{(s)},Y^{(s)})$ denotes a set of the covariates and the
responses for the $s$th dataset. Notice that $X^{(1)}$ and $X^{(2)}$ can be
from different domains with different data types and dimensions of features.
The $P^{(s)}$ denotes the probability mass distribution on the samples from
the $s$th dataset to represent relative importance for the corresponding
samples. If no prior information is known on the $P^{(s)}$, one can set the
empirical distribution version of
$P^{(s)}=\frac{1}{N_{s}}\mathbbm{1}^{\top}_{N_{s}}$ with mass of each sample
point being $\frac{1}{N_{s}}$, where
$\mathbbm{1}_{N_{s}}=(1,1,\cdots,1)_{N_{s}}$ and $N_{s}$ is the sample size of
the $s$th dataset. To quantify the discrepancy between distributions on
different domain spaces, [25] proposes the Gromov-Wasserstein distance on
$\\{X^{(s)},P^{(s)}\\},(s=1,2)$ to measure the geometric structural similarity
of two distributions. Specifically, one constructs the distance matrix
$D^{(s)}=(d^{(s)}_{ij}),(s=1,2)$ to represent similarities among samples
within the $s$th dataset where $d^{(s)}_{ij}=d^{(s)}(X^{(s)}_{i},X^{(s)}_{j})$
is the distance between $X^{(s)}_{i}$ and $X^{(s)}_{j}$ corresponding to a
given metric $\bm{X}\times\bm{X}\rightarrow R^{+}$ in the $s$th space. The
Gromov-Wasserstein distance is defined as

$\text{GW}(D^{(1)},D^{(2)},P^{(1)},P^{(2)})=\min_{\Omega\in\bm{\Omega}_{\\{P^{(1)},P^{(2)}\\}}}\sum_{i,j,k,l}L\left(D^{(1)}_{i,k},D^{(2)}_{j,\ell}\right)\cdot\Omega_{i,j}\cdot\Omega_{k,\ell},$

where $D^{(1)}_{i,k}$ and $D^{(2)}{j,l}$ denote the $(i,k)$ element and
$(j,l)$ element in $D^{(1)}$ and $D^{(2)}$, and $\Omega=(\Omega_{ij})$ is a
transport coupling matrix between two datasets satisfying

$\Omega\in\bm{\Omega}_{\\{P^{(1)},P^{(2)}\\}}=\left\\{\Omega\in\left(\mathbb{R}_{+}\right)^{N_{1}\times
N_{2}}:\Omega\mathbbm{1}^{\top}_{N_{2}}=P^{(1)},\Omega^{\top}\mathbbm{1}^{\top}_{N_{1}}=P^{(2)}\right\\}.$

Here $L$ is a loss function to account for the element-wise discrepancy
between the distance matrices. Typical choices of $L$ include the quadratic
loss $L(a;b)=(a-b)^{2}$ and Kullback-Leibler divergence
$L(a;b)=a\log(a/b)-a+b$.

Fused Gromov-Wasserstein distance To further measure the discrepancy of joint
distribution $(X^{(s)},Y^{(s)})$ between two datasets, the fused Gromov-
Wasserstein distance is a generalization of the aforementioned Gromov-
Wasserstein distance. To incorporate the response information $Y^{(s)}$ into
the alignment between two datasets, we can first construct
$M=\big{(}d(Y^{(1)}_{i},Y^{(2)}_{j})\big{)}_{ij}$ which is a $N_{1}\times
N_{2}$ distance matrix, and $d(\cdot,\cdot)$ stands for a similarity
measurement between responses within the same metric space. Then the fused
Gromov-Wasserstein distance is defined as:

$\displaystyle\text{FGW}(D^{(1)},D^{(2)},M,P^{(1)},P^{(2)})=\min_{\Omega\in\bm{\Omega}_{\\{P^{(1)},P^{(2)}\\}}}\sum_{i,j,k,l}\big{\\{}(1-\alpha)M_{ij}+\alpha
L\left(D^{(1)}_{i,k},D^{(2)}_{j,\ell}\right)\big{\\}}\\!\cdot\\!\Omega_{i,j}\\!\cdot\\!\Omega_{k,\ell},$

where $0\leq\alpha\leq 1$. The FGW distance searches for a coupling $\Omega$
between two datasets that minimizes a weighted cost combining the mass
transportation on response from $Y^{(1)}_{i}$ to $Y^{(2)}_{j}$, and on
covariate pairs from $D^{(1)}_{i,k}$ to $D^{(2)}_{j,l}$. The $i$th column in
the coupling matrix $\Omega$ decides the proportions of probability mass in
the $i$th sample point from the first dataset to be transported to each sample
point in the second dataset. Notice that when $\alpha=1$, the FGW distance
degenerates to the GW distance.

Optimal transport for data integration The optimal transport can be applied to
cross-dataset alignment. We consider two different covariate spaces
$(X^{(1)},d^{{(1)}},\mu^{(1)})$ and $(X^{(2)},d^{{(2)}},\mu^{(2)})$ where
$d^{(1)},d^{(2)}$ and $\mu^{(1)},\mu^{(2)}$ denote the metrics and probability
measures on each covariate space. The optimal transport can be used to align
the geometric structures of $\mu^{(1)}$ and $\mu^{(2)}$. Specifically, we
randomly sample pairs of covariates $\\{(x_{i},x_{j})\\}$ from
${X}^{(1)}\times{X}^{(1)}$ and pairs of covariates $\\{(y_{k},y_{l})\\}$ from
${X}^{(2)}\times{X}^{(2)}$, respectively. Given that the pairwise metrics
$\\{d^{(1)}(x_{i},x_{j})\\}$ and $\\{d^{(2)}(y_{k},y_{l})\\}$ capture the
corresponding geometric information of $\mu^{(1)}$ and $\mu^{(2)}$, we aim to
seek an alignment plan between $\mu^{(1)}$ and $\mu^{(2)}$ to minimize the
discrepancy between the two metrics
$L((x_{i},x_{j}),(y_{k},y_{l}))=|d^{(1)}(x_{i},x_{j})-d^{(2)}(y_{k},y_{l})|$
at the population level. Based on the set of joint distributions
$\bm{\Omega}=\\{\Omega\in
P({X}^{(1)}\times{X}^{(2)})\mid\Omega({X}^{(1)}\times
A)=\mu^{(2)}(A),\Omega(B\times{X}^{(2)})=\mu^{(1)}(B)\\}$, the optimal
alignment plan to minimize geometric structure discrepancy between $\mu^{(1)}$
and $\mu^{(2)}$ is defined as:

$\displaystyle\Omega=\underset{\Omega\in\bm{\Omega}}{\operatorname{argmin}}\int_{({X}^{(1)}\times{X}^{(2)})\times({X}^{(1)}\times{X}^{(2)})}\\!\\!\\!\\!\\!\\!L((x_{i},x_{j}),(y_{k},y_{l}))d\Omega\left(\mathbf{x}_{i},\mathbf{y}_{k}\right)d\Omega\left(\mathbf{x}_{j},\mathbf{y}_{l}\right).$
(1)

Notice that the metric discrepancy $L((x_{i},x_{j}),(y_{k},y_{l}))$ is defined
on the product space, therefore we use the product of alignment plan
$\Omega\left({x}_{i},{y}_{k}\right)\times\Omega\left({x}_{j},{y}_{l}\right)$
to serve as the alignment between sample pairs $(x_{i},x_{j})$ and
$(y_{k},y_{l})$. The alignment plan obtained in (1) is well-defined and
induces the Gromov–Wasserstein distance which generalizes the Wasserstein
distance. We extend the above intuition into the proposed methods in Section
3.

##  3 Methodology

In this section, we propose a new data integration framework ILA for both
supervised and unsupervised learning tasks. The proposed ILA framework
consists of three stages: data compression, cross-subject alignment, and
common pattern learning, which is illustrated in Figure 2. The data
compression in stage 1 of the ILA framework directly adopts existing dimension
reduction methods such as principal component analysis [23] or the autoencoder
[21]. Autoencoders have been used to compress neural activity data in other
applications [35, 30, 10, 29] and empirical studies have demonstrated their
superior performance in capturing non-linear features from high-dimensional
data compared with other unsupervised dimensionality reduction methods such as
PCA [5], Isomap, and locally linear embedding [34]. Therefore, we utilize
autoencoder in stage 1, and focus on the key ingredients stage 2 and stage 3
of our method that utilizes the idea of incorporating the heterogeneity
domains in datasets by aligning them into a common data space. We first
introduce a cross-domain data integration method for the supervised learning
task in Section 3.1, and then introduce an integration method for the cross-
domain temporal patterns in Section 3.2.

Figure 2: The proposed integrated latent alignment framework (ILA). In stage
1, we utilize an autoencoder to compress the electrophysiological data for
each subject. In stage 2, we use optimal transport to align the compressed
data from different subjects in latent spaces. In the last stage, a supervised
learning algorithm is implemented in the aligned latent space to extract
common patterns.

###  3.1 Integrated supervised learning on multiple heterogeneous datasets

In this subsection, we consider multiple datasets
$\\{(X^{(s)},Y^{(s)})\\}_{s=1}^{S}$ defined in $S$ different domains. The
joint distributions $\mathbb{P}^{(s)}(\bm{X}^{(s)},\bm{Y}^{(s)})$ on the $s$th
dataset can be distinct from one another due to the heterogeneity among
marginal distributions $\mathbb{P}^{(s)}(\bm{X}^{(s)}),\;s=1,2,\cdots,S$. To
enable the cross-domain alignment, we assume that the marginal distributions
$\mathbb{P}^{(s)}(\bm{X}^{(s)}),\;s=1,2,\cdots,S$ share the same geometric
structure in that for two datasets $s$ and $s^{\prime}$: there exists a
probabilistic mapping
$\Gamma^{s^{\prime},s}:\bm{X}^{(s^{\prime})}\times\bm{Y}^{(s^{\prime})}\to\bm{X}^{(s)}\times\bm{Y}^{(s)}$
satisfying

$\displaystyle\mathbb{P}^{(s)}\big{(}\bm{X}^{(s)},\bm{Y}^{(s)}\big{)}\approx\mathbb{P}^{(s^{\prime})}\big{(}\Gamma^{s^{\prime},s}(\bm{X}^{(s^{\prime})},\bm{Y}^{(s^{\prime})})\big{)},1\leq
s\neq s^{\prime}\leq S.$ (2)

Notice that the measure-preserving mapping $\Gamma^{s^{\prime},s}$ corresponds
to the alignment plan $\Omega$ in (1) between two spaces
$\bm{X}^{(s^{\prime})}\times\bm{Y}^{(s^{\prime})}$ and
$\bm{X}^{(s)}\times\bm{Y}^{(s)}$. In addition, the geometry information of
$\mathbb{P}^{(s)}(\bm{X}^{(s)})$ can be represented as a set of pairwise
distances $\\{d^{(s)}(x_{i},x_{j})\\}$ where $x_{i},x_{j}$ are samples from
$\mathbb{P}^{(s)}(\bm{X}^{(s)})$. Then the distribution of
$\\{d^{(s)}(x_{i},x_{j})\\}$ represents the concentration of samples from
$\bm{X}^{(s)}$ at different regions. Therefore, a common geometric structure
between $\bm{X}^{(s)}$ and $\bm{X}^{(s^{\prime})}$ implies that the mapping
$\Gamma^{s^{\prime},s}$ preserves distances between the two domains as

$\displaystyle d^{(s^{\prime})}\big{(}x_{i},x_{j}\big{)}\approx
d^{(s)}\big{(}\Gamma^{s^{\prime},s}(x_{i}),\Gamma^{s,s^{\prime}}(x_{j})\big{)},1\leq
s\neq s^{\prime}\leq S.$

The cross-domain relation (2) guarantees the existence of a shared
distribution of response conditional on predictors across $S$ datasets as
follows:

$\displaystyle\mathbb{P}^{(s)}\big{(}\bm{Y}^{(s)}|\bm{X}^{(s)}\big{)}\approx\mathbb{P}^{(s^{\prime})}\big{(}\bm{Y}^{(s^{\prime})}|\Gamma^{s^{\prime},s}(\bm{X}^{(s^{\prime})})\big{)},1\leq
s\neq s^{\prime}\leq S.$

This enables us to increase the effective training sample size in the $s$
dataset via integrating labeled data from other datasets $s^{\prime}$
associated with $s$, whereby the performance of the machine learning model of
the $s$th domain can be improved.

Without loss of generality, we consider the classification task as follows.
For the $s$th dataset $1\leq s\leq S$, let the training dataset and the set of
unlabeled sample points be denoted as $(X_{l}^{(s)},Y_{l}^{(s)})$ and
$X_{p}^{(s)}$, respectively. We denote the sample points in the training set
and unlabeled dataset as $(X_{l,i}^{(s)},Y_{l,i}^{(s)})$ and $X_{p,i}^{(s)}$
such that
$X_{l}^{(s)}=\\{X_{l,i}^{(s)}\\}_{i=1}^{n_{l}^{(s)}},Y_{l}^{(s)}=\\{Y_{l,i}^{(s)}\\}_{i=1}^{n^{(s)}_{l}}$
and $X_{p}^{(s)}=\\{X_{p,i}^{(s)}\\}_{i=1}^{n^{(s)}_{p}}$, where $n^{(s)}_{l}$
and $n^{(s)}_{p}$ are the number of training samples and unlabeled samples in
the $s$th dataset, respectively. In the following, we consider the data
integration for the classification task on the $s$th dataset. The data
integration on the other datasets follows the same procedure. The proposed
method is summarized as follows:

Step 1: Individual supervised learning Given a specific classifier $f(\cdot)$,
we train the classifier on $(X_{l}^{(s)},Y_{l}^{(s)})$ separately as
$f_{s}(\cdot)$, and then predict on the unlabeled set $X^{(s)}_{p}$ to obtain
predicted labels $\hat{Y}^{(s)}_{p}=f_{s}(X^{(s)}_{p})$ for each dataset
$s=1,\cdots,S$.

Step 2: One-to-many cross-domain alignment We consecutively align training
sets from other $S-1$ datasets
$(X_{l}^{(s^{\prime})},Y_{l}^{(s^{\prime})}),\;s^{\prime}\neq s,1\leq
s^{\prime}\leq S$ to the region of unlabeled data $X^{(s)}_{p}$ in the $s$th
dataset. Specifically, we adopt the fused Gromov-Wasserstein optimal transport
as follows:

$\displaystyle\Omega^{s,s^{\prime}}=\operatorname*{argmin}_{\bm{\Omega}}\sum_{i,j,k,l}\big{\\{}(1-\alpha)M^{s,s^{\prime}}_{ij}+\alpha
L(D^{(s)}_{i,k},D^{(s^{\prime})}_{j,\ell})\big{\\}}\Omega_{i,j}\Omega_{k,\ell},$
(3)
$\displaystyle\bm{\Omega}=\left\\{\Omega\in\left(\mathbb{R}_{+}\right)^{n_{p}^{(s)}\times
n_{l}^{(s^{\prime})}}:\Omega\mathbbm{1}_{n_{l}^{(s^{\prime})}}=\frac{1}{n_{p}^{(s)}}\mathbbm{1}_{n_{p}^{(s)}},\Omega^{\top}\mathbbm{1}_{n_{p}^{(s)}}=\frac{1}{n_{l}^{(s^{\prime})}}\mathbbm{1}_{n_{l}^{(s^{\prime})}}\right\\},$

where
$M^{s,s^{\prime}}=\big{(}d(\hat{Y}^{(s)}_{p,i},{Y}_{l,j}^{(s^{\prime})})\big{)}$,
$D^{(s)}=\big{(}d^{(s)}(X^{(s)}_{p,i},X^{(s)}_{p,j})\big{)}$, and
$D^{(s^{\prime})}=\big{(}d^{(s^{\prime})}(X^{(s^{\prime})}_{l,i},X^{(s^{\prime})}_{l,j})\big{)}$
are distance matrices, $d(\cdot)$ is the similarity measure of the labels, and
$d^{(s)}(\cdot)$ denotes the distance between predictors on the $s$th domain.
The $i$th column in the coupling matrix $\Omega^{s,s^{\prime}}$ decides the
proportions of probability mass in the training sample point
$(X^{(s^{\prime})}_{l,i},Y^{(s^{\prime})}_{l,i})$ to be transported to each
unlabeled sample point in $\\{X^{(s)}_{p,i}\\}$. We choose L2 norm
$L(\cdot)=\|\cdot\|^{2}$ in (3).

Given that $X_{l}^{(s^{\prime})}$ and $X_{p}^{(s)}$ are not in the same
domain, $\Omega^{s,s^{\prime}}$ intends to match the internal geometric
structures between $X_{l}^{(s)}$ and $X_{p}^{(s^{\prime})}$, which are
measured by distance matrices $D^{(s)}$ and $D^{(s^{\prime})}$, respectively.
On the other hand, $\Omega^{s,s^{\prime}}$ also aims to align the label
distribution between datasets $s$ and $s^{\prime}$ though the term
$(1-\alpha)M^{s,s^{\prime}}_{ij}$. Minimizing (3) allows us to match the joint
distribution $\mathbb{P}(X_{l}^{(s^{\prime})},Y_{l}^{(s^{\prime})})$ to
$\mathbb{P}(X_{p}^{(s)},Y_{p}^{(s)})$ so that (2) can be held where
$\Omega^{s,s^{\prime}}$ serves as the underlying cross-domain transformation
$\Gamma^{s,s^{\prime}}$, and hence we can integrate training data from
heterogeneous domains to capture the common label-predictor pattern. Because
$X_{p}^{(s)}$ are not labeled, i.e., the corresponding $Y_{p}^{(s)}$ are not
observed, we approximate $\mathbb{P}(X_{p}^{(s)},Y_{p}^{(s)})$ by
$\mathbb{P}(X_{p}^{(s)},\hat{Y}_{p}^{(s)})$ where, $\hat{Y}_{p}^{(s)}$ are
obtained in Step 1.

Step 3: Training data mapping We map the training set
$(X_{l}^{(s^{\prime})},Y_{l}^{(s^{\prime})})\;(s^{\prime}\neq s)$ from other
$S-1$ datasets to the domain of the $s$th dataset. Specifically, we
reconstruct the $j$th training sample in the $s^{\prime}$th dataset
$(X_{l,j}^{(s^{\prime})},Y_{l,j}^{(s^{\prime})})$ on the $s$th domain as
$(\hat{X}_{l,j}^{(s^{\prime})},Y_{l,j}^{(s^{\prime})})$, where

$\displaystyle\hat{X}_{l,j}^{(s^{\prime})}=\underset{X}{\operatorname{argmin}}\sum_{i}\Omega^{s,s^{\prime}}_{ij}d^{(s)}\left(X,X^{(s)}_{p,i}\right),$
(4)

and $\Omega^{s,s^{\prime}}$ is obtained from Step 2. The barycentric mapping
(4) corresponds to a weighted average of unlabeled data $X^{(s)}_{p}$ in the
$s$th domain. Intuitively, the barycentric mapping reconstructs the training
sample from the $s$the dataset in the domain of the $s^{\prime}$th dataset. We
denote the reconstructed samples as
$\hat{X}^{(s^{\prime})}_{l}=\\{\hat{X}^{(s^{\prime})}_{l,j}\\}$. Based on Step
2 and assumption (2), $(\hat{X}^{(s^{\prime})}_{l},Y^{(s^{\prime})}_{l})$
follows the same distribution in the $s$ domain as
$(X^{(s)}_{p},Y^{(s)}_{p})$.

Step 4: Retrain the model with augmented training data We obtain reconstructed
training data from other $S-1$ datasets as
$(\tilde{X}^{(s)}_{l},\tilde{Y}^{(s)}_{l})=\underset{s^{\prime}\neq
s}{\cup}(\hat{X}^{(s^{\prime})}_{l},Y^{(s^{\prime})}_{l})$, and train the
classifier $f(\cdot)$ on the augmented training dataset
$(X^{(s)}_{l}\cup\tilde{X}^{(s)}_{l},Y^{(s)}_{l}\cup\tilde{Y}^{(s)}_{l})$ as
$f_{s}^{aug}(\cdot)$. Then we update the prediction on unlabeled data
$X^{(s)}_{p}$ as $\hat{Y}^{(s)}_{p}=f_{s}^{aug}(X^{(s)}_{p})$.

In summary, the proposed data integration method is designed for the semi-
supervised learning scenario, where the training data from each dataset
consists of a small proportion of labeled samples
$\\{X^{(s)}_{l},Y^{(s)}_{l}\\}_{s=1}^{S}$ and a large number of unlabeled
samples $\\{X^{(s)}_{p}\\}_{s=1}^{S}$. The objective is to predict outcomes
for the unlabeled samples.

###  3.2 Temporal integration of multiple heterogeneous dataset

In this subsection, we introduce the proposed data integration for supervised
learning on multiple longitudinal heterogeneous datasets. Specifically, we
consider $S$ longitudinal heterogeneous datasets
$\\{\big{(}X^{(s)}(t),Y^{(s)}(t)\big{)}\\}_{s=1}^{S}$ measured at time points
$t=1,\cdots,T$, where covariates $X^{(s)}(t)$, $s=1,\cdots,S$ change over
time, but are defined on different domains.

To perform data integration, we denote the temporal alignment mappings
$\Gamma^{s^{\prime},s}$ between domains $s^{\prime}$ and $s$ as:

$\displaystyle\mathbb{P}^{(s)}\big{(}\bm{X}^{(s)}(t),\bm{Y}^{(s)}(t)\big{)}\approx\mathbb{P}^{(s^{\prime})}\Big{(}\Gamma^{s^{\prime},s}\big{(}\bm{X}^{(s^{\prime})}(t),\bm{Y}^{(s^{\prime})}(t)\big{)}\Big{)},1\leq
s\neq s^{\prime}\leq S,\;1\leq t\leq T.$ (5)

We consider the case where the alignment mappings
$\\{\Gamma^{s^{\prime},s}\\}$ are time invariant, and hence (5) implies that
the common temporal pattern is identifiable among the $S$ different datasets,
i.e.,

$\displaystyle\mathbb{P}^{(s)}\big{(}\bm{X}^{(s)}(t),\bm{Y}^{(s)}(t)\big{)}\approx\mathbb{P}\Big{(}\Gamma^{s}\big{(}\bm{Z}(t),\bm{Y}(t)\big{)}\Big{)},1\leq
s\leq S,\;1\leq t\leq T,$ (6)

where $\\{\bm{Z}(t)\\}_{t=1}^{T}$ and $\\{\bm{Y}(t)\\}_{t=1}^{T}$ denote the
covariate-associated and response-associated dynamic process on the aligned
data domains, respectively. In other words, the common temporal pattern on the
conditional distribution of response given that covariates is captured by
$\\{(\bm{Z}(t),\bm{Y}(t)\\}_{t=1}^{T}$. Then the $s$th mapping
$\\{\Gamma^{s}\\}_{s=1}^{S}$ projects the common temporal patterns to domains
of the $s$th dataset. By inferring $\\{(\bm{Z}(t),\bm{Y}(t)\\}_{t=1}^{T}$ and
$\\{\Gamma^{s}\\}_{s=1}^{S}$, we can integrate cross-domain temporal data into
a common latent space to improve the performance of the supervised learning
task as the effective training sample size increases.

In the following, we denote the labeled training dataset as
$\big{(}X^{(s)}_{l}(t),Y^{(s)}_{l}(t)\big{)}$ and the unlabeled dataset as
$X^{(s)}_{p}(t)$, such that
$X^{(s)}_{l}(t)=\\{X^{(s)}_{l}(t)_{i}\\}_{i=1}^{n_{l}^{(s)}}$,
$Y^{(s)}_{l}(t)=\\{Y^{(s)}_{l}(t)_{i}\\}_{i=1}^{n_{l}^{(s)}})$ and
$X^{(s)}_{p}(t)=\\{X^{(s)}_{p}(t)_{i}\\}_{i=1}^{n_{p}^{(s)}}$, $1\leq s\leq
S,1\leq t\leq T$. The number of total samples in the $s$th dataset is
$n^{(s)}=n_{l}^{(s)}+n_{p}^{(s)}$. Let
$X^{(s)}(t)=(X_{l}^{(s)}(t),X_{p}^{(s)}(t))$ be the set of covariates in the
$s$th domain. The proposed method is summarized as follows:

Step 1: Individual supervised learning Given a specific classifier $f(\cdot)$,
we train the classifier on $(X_{l}^{(s)}(t),Y_{l}^{(s)}(t))$ separately as
$f^{(s)}_{t}(\cdot)$ for each dataset $s$ and time point $t$, $1\leq s\leq
S,1\leq t\leq T$. Then obtain predicted labels through unlabeled data
$X^{(s)}_{p}(t)$ via
$\hat{Y}^{(s)}_{p}(t)=f_{t}^{(s)}\big{(}X^{(s)}_{p}(t)\big{)}$ for each
dataset and time point.

Step 2: Cross-domain temporal alignment We align the time-evolving datasets
from $S$ different domains via the proposed temporal Fused Gromov-Wasserstains
(FGW) barycenter. Denoting the augmented labels for the $s$th dataset as
$Y^{(s)}(t)=\big{(}Y^{(s)}_{l}(t),\hat{Y}^{(s)}_{p}(t)\big{)},1\leq t\leq T$,
we estimate the alignment couplings $\\{\Omega^{(s)}\\}_{s=1}^{S}$ and the
common dynamic process $\\{C(t)\\}_{t=1}^{T}=\\{D(t),Y(t)\\}_{t=1}^{T}$ in the
following aligned space:

$\displaystyle\\{C(t)\\}_{t=1}^{T},$
$\displaystyle\\{\Omega^{(s)}\\}_{s=1}^{S}=\operatorname*{argmin}_{\begin{subarray}{c}C(t),t=1,\cdots,T\\\
\Omega^{(s)}\in\bm{\Omega}^{(s)},s=1,\cdots,S\end{subarray}}\sum_{t=1}^{T}\sum_{s=1}^{S}w_{s}\;\text{FGW}\Big{(}C(t),\Omega^{(s)}\Big{)},$
(7) $\displaystyle\text{where FGW}\Big{(}C(t),\Omega^{(s)}\Big{)}$
$\displaystyle\\!=\\!\alpha\\!\sum_{ij}\big{(}M^{(s)}(t)\big{)}_{ij}\Omega^{(s)}_{ij}\\!+\\!(1\\!-\\!\alpha)\\!\sum_{i,j,u,v}\\!\|\big{(}D^{(s)}(t)\big{)}_{ij}\\!-\\!\big{(}D(t)\big{)}_{uv}\|^{2}\Omega^{(s)}_{iu}\Omega^{(s)}_{jv},$
$\displaystyle\bm{\Omega}^{(s)}$
$\displaystyle=\left\\{\Omega\in\left(\mathbb{R}_{+}\right)^{n^{(s)}\times
n}:\Omega\mathbbm{1}_{n}=\frac{1}{n^{(s)}}\mathbbm{1}_{n^{(s)}},\Omega^{\top}\mathbbm{1}_{n^{(s)}}=\frac{1}{n}\mathbbm{1}_{n}\right\\},$

where $M^{(s)}(t)=\Big{(}d\big{(}Y^{(s)}(t)_{i},Y(t)_{j}\big{)}\Big{)}$,
$D^{(s)}(t)=\Big{(}d^{(s)}\big{(}X^{(s)}(t)_{i},\;X^{(s)}(t)_{j}\big{)}\Big{)},\;1\leq
t\leq T$, $\alpha\in[0,1]$ is the trade-off parameter, and $n$ is the number
of template samples in the common latent space. The similarity matrix $D(t)\in
R_{+}^{n\times n}$ represents the pairwise distances aligned with $D^{(s)}(t)$
from different datasets at time $t$, and consists of the common covariates
dynamics $Z(t)$ in (6). Specifically, given the positions of $n$ template
samples $\\{Z(t)_{i}\\}_{i=1}^{n}$ in the aligned latent space at time $t$, we
have $D(t)_{ij}=\tilde{d}\big{(}Z(t)_{i},Z(t)_{j}\big{)}$, where
$\tilde{d}(\cdot)$ is the distance between $Z(t)_{i}$ and $Z(t)$ in the latent
space. In addition, the transport coupling
$\\{\Omega^{(s)}\in[0,1)^{(n_{l}^{(s)}+n_{p}^{(s)})\times n}\\}$ serves as the
domain transformation $\\{\Gamma^{s}\\}$ in (6) which assigns the probability
mass on $\\{X^{(s)}(t)_{i},Y^{(s)}(t)_{i}\\}_{i=1}^{n_{l}^{(s)}+n_{p}^{(s)}}$
for template samples $\\{Z(t)_{i},Y(t)_{i}\\}_{i=1}^{n}$ in the latent space
with $1\leq t\leq T$. Finally, $w_{s}\geq 0$ in (7) denotes the weights of the
$s$th dataset in the alignment as $\sum_{s=1}^{S}w_{s}=1$, and we set
$w_{s}=\frac{1}{S}$ by default.

Step 3: Data integration on aligned space We embed the aligned similarity
matrix $D(t)$ to a $k$-dimensional latent space at each time point $t$.
Specifically, we infer $Z(t)=\\{Z(t)_{i}\in R^{k}\\}_{i=1}^{n}$ via

$\displaystyle
Z(t)=\operatorname*{argmin}_{\bm{Z}}\sum_{i,j}^{n}(D(t)_{ij}-\|\bm{Z}_{i}-\bm{Z}_{j}\|_{2})^{2}.$

Then we map all $s$ datasets to latent space $\\{Z(t)\\}_{t=1}^{T}$, i.e., the
training sample $\big{(}X^{(s)}_{l}(t)_{i},Y^{(s)}_{l}(t)_{i}\big{)}$ and the
$j$th testing sample $X^{(s)}_{p}(t)_{j}$ in the $s$th dataset are
reconstructed in the latent space as
$\big{(}\hat{X}^{(s)}_{l}(t)_{i},Y^{(s)}_{l}(t)_{i}\big{)}$ and
$\hat{X}^{(s)}_{p}(t)_{j}$ respectively, where

$\displaystyle\hat{X}_{l}^{(s)}(t)_{i}=\underset{x}{\operatorname{argmin}}\sum_{u}\Omega^{s}_{iu}d\big{(}x,Z(t)_{u}\big{)},\;\hat{X}_{p}^{(s)}(t)_{j}=\underset{x}{\operatorname{argmin}}\sum_{u}\Omega^{s}_{ju}d\big{(}x,Z(t)_{u}\big{)},$

with $1\leq s\leq S,1\leq t\leq T$, and $\Omega^{s}$ is obtained from Step 2.
In conjunction with data integration, inferring $\\{Z(t),Y(t)\\}$ in a low-
dimensional latent space also allows us to identify and model the governing
temporal pattern on conditional distribution of the response after integrating
heterogeneous data. This strategy can be useful for other downstream tasks and
applications, such as clustering and pattern exploration.

Step 4: Retrain the model with augmented training data For the $s$th dataset
at time $t$, we aggregate training data from all $S$ datasets reconstructed in
the common latent space as

$\big{(}\tilde{X}^{(s)}_{l}(t),\tilde{Y}^{(s)}_{l}(t)\big{)}=\underset{1\leq
s\leq S}{\cup}\big{(}\hat{X}_{l}^{(s)}(t),Y^{(s)}_{l}(t)\big{)},$

and train the classifier $f(\cdot)$ on the augmented training dataset
$\big{(}\tilde{X}^{(s)}_{l}(t),\tilde{Y}^{(s)}_{l}(t)\big{)}$ as
$\tilde{f}_{t}^{(s)}(\cdot)$. Consequently, we update the prediction on
unlabeled data $X^{(s)}_{p}(t)$ as
$\hat{Y}^{(s)}_{p}(t)=\tilde{f}_{t}^{(s)}\big{(}\hat{X}^{(s)}_{p}(t)\big{)}$.

Remark: our proposed method relaxes the assumption on the bijection mapping
between joint distribution $(X,Y)$ cross-datasets. Instead we assume that
there exist common clustering structures among distributions of covariates $X$
from different datasets. Specifically, for covariates $X^{(s)}$ and
$X^{(s^{\prime})}$ from the $s$th and $s^{\prime}$th dataset, there exists a
cross-dataset mapping $\Gamma:\bm{X}^{(s)}\to\bm{X}^{(s^{\prime})}$ which
approximately preserves the metric information as

$\displaystyle d^{\left(s\right)}\left(X^{(s)}_{i},X^{(s)}_{j}\right)\approx
d^{(s^{\prime})}\left(\Gamma\left(X^{(s)}_{i}\right),\Gamma\left(X^{(s)}_{j}\right)\right),$
(8)

where $d^{(s)}(\cdot)$ and $d^{(s^{\prime})}(\cdot)$ are the metrics on
dataset $s$ and $s^{\prime}$, respectively. Under the context of neuroscience
applications where $\\{X^{(s)}\\}$ stand for neural activities in different
subjects, condition (1) is a reasonable assumption based on recent findings
that the correlations among neurons generally lead to strong clustering
organizations in neural activity patterns on the population level [36, 11].
These underlying common clustering structures among neural activities imply
the existence of metric preserve $\Gamma(\cdot)$. On the other hand, the
marginal distribution $\mathbb{P}^{(s)}(\Gamma(\bm{X}^{(s)}))$ can still be
different from $\mathbb{P}^{(s^{\prime})}(\bm{X}^{(s^{\prime})})$ due to
subjects’ heterogeneity.

In addition, the proposed method assumes that correspondence between response
$Y$ and clusters of covariates $X$ can be preserved across different datasets.
Specifically, we assume that the condition distributions of response $Y$ are
invariant as

$\displaystyle\mathbb{P}^{(s)}(\bm{Y}\mid\bm{X}^{(s)})\approx\mathbb{P}^{(s^{\prime})}(\bm{Y}\mid\Gamma(\bm{X}^{(s)})).$
(9)

Notice that (9) is the covariate shift assumption, which is generally made in
most domain adaption problems [13]. Under the context of neuroscience
applications, existing studies show that the clustering structures of neural
activities generally serve as common neural codes carrying external sensory
information and neural processing statuses information[8, 22, 19], which
support the assumption (9) for our real application. The detailed optimization
procedure and hyperparameter tuning are summarized in the Section 8 of
Supplemental Materials.

##  4 Simulation Study

In this section, we investigate the performance of the proposed data
integration approach and temporal integration method for several
classification problems, and compare the performance to a baseline model that
uses the data from individuals. Here we choose two classifiers for the
numerical comparisons. The first is the quadratic discriminative analysis
(QDA), and the second is the label propagation algorithm (LP); the former is a
supervised learning classifier, whereas the latter is a semi-supervised
learning method utilizing the structure information on covariates.

###  4.1 Simulations of data integration for classification

In this simulation setting, we consider six different datasets where the label
of data points $Y$ in each dataset can belong to one of three classes. The
data generation is as follows. For each class $Y=k,(k=1,2,3)$, the two-
dimensional covariates $X$ are sampled from the mixture of Gaussian
distributions as $X\mid(Y=k)\sim 0.5\times N(\mu^{(1)}_{k},I)+0.5\times
N(\mu^{(2)}_{k},I)$ where $I=\big{(}(1,0)^{T},(0,1)^{T}\big{)}$ is the
covariance matrix, and $\mu^{(1)}_{k},\mu^{(2)}_{k}$ denote the class-specific
means. Specifically, we set

$\displaystyle\mu^{(1)}_{1}=(-a,a),\mu^{(2)}_{1}=(a,a),\mu^{(1)}_{2}=(2a,0),\mu^{(2)}_{2}=(a,-a),\mu^{(1)}_{3}=(-a,-a),\mu^{(2)}_{3}=(-2a,0),$
(10)

where $a$ represents the separability of $X$ among three classes. For each
dataset, we randomly generate $N$ samples, where $Y$ follows a multinomial
distribution with three classes of equal probability, and $X$ are sampled from
the previous mixture Gaussian distributions. Then we obtain data points
$(\Gamma_{l}X,Y)$ as the $l$th dataset, where $\Gamma_{l}(l=1\cdots,6)$ are
the rotation matrices $\Gamma_{l}=\big{(}(\text{cos $\theta_{l}$},-\text{sin
$\theta_{l}$})^{T},(\text{sin $\theta_{l}$},\text{cos
$\theta_{l}$}^{T}\big{)}$, and $\theta_{l}=20\times(l-1)$. Notice that the
covariates $X$ from each dataset have six clusters corresponding to different
labels $Y$; we index these clusters as $c=1,\cdots,6$. Next we split data
points from each dataset into a training set with
$N_{train}=\sum_{s=1}^{6}n_{l}^{(s)}$ data points, and a testing set with
$N_{test}=\sum_{s=1}^{6}n_{p}^{(s)}$ data points, where
$N_{train}+N_{test}=N$. For the first dataset, we sample $90\%$ of the
training data from clusters $c=\\{1,3,6\\}$ and $10\%$ from the other three
clusters. Then the remaining data points serve as testing data. The generation
of the training and testing data for other five datasets follows the same
procedure, while the training data are sampled from different clusters.
Specifically, the cluster indexes of training data for the second to sixth
dataset are $\\{1,4,6\\},\\{2,3,5\\},\\{2,3,6\\}$,$\\{1,4,5\\},\\{2,4,5\\}$,
respectively. The above data rotation and training data sampling scheme
increase the heterogeneity of the distribution of $X$ across different
datasets.

In the following, we compare the classification performance of the data
integration procedure in Section 3.1 with that of the individual learning
scheme. For individual learning, we independently train a classifier on the
training data, and then make predictions on the testing data from each
dataset. While for the data integration scheme, we train the classifier for
the $l$th dataset by using the training data from both the $l$th dataset and
$M$ datasets. The number of datasets for integration varies as $M=1,3,5$. We
then calculate the misclassification rate averaged from the 6 datasets using
two different schemes. The comparisons are conducted under different sizes of
training set $N_{train}$ and testing set $N_{test}$. In addition, we measure
the separability intensity of the classification as pairwise distances among
cluster centers as $\Delta:=\sum_{1\leq k_{1}\leq k_{2}\leq 3}\sum_{1\leq
i\leq j\leq 2}\|\mu^{i}_{k_{1}}-\mu^{j}_{k_{2}}\|_{2}$, and denote the
intensity as $\Delta_{0}$, $\Delta_{1}$, $\Delta_{2}$ when $a$ in (10) varies
among $0.7,1,1.5$. A larger $a$ indicates a strong separation of
classifications.

Signal intensity $N_{train}/N_{test}$ LP QDA Ind. Prop. Ind. Prop. $M=1$ $M=3$
$M=5$ $M=1$ $M=3$ $M=5$ $\Delta_{0}$ 90/120 0.33(0.02) 0.30(0.02) 0.29(0.02)
0.29(0.02) 0.32(0.02) 0.29(0.02) 0.28(0.02) 0.28(0.02) 150/150 0.33(0.02)
0.30(0.02) 0.29(0.02) 0.28(0.02) 0.33(0.02) 0.29(0.02) 0.28(0.02) 0.27(0.02)
240/180 0.36 (0.01) 0.33(0.02) 0.31(0.02) 0.30(0.02) 0.37(0.01) 0.30(0.02)
0.29(0.02) 0.28(0.02) $\Delta_{1}=1.4{\Delta_{0}}$ 90/120 0.25(0.02)
0.20(0.02) 0.18(0.02) 0.18(0.02) 0.23(0.02) 0.19(0.02) 0.17(0.02) 0.17(0.02)
150/150 0.26(0.01) 0.20(0.02) 0.18(0.02) 0.18(0.02) 0.24(0.02) 0.18(0.02)
0.16(0.02) 0.16(0.01) 240/180 0.32 (0.02) 0.25(0.03) 0.23(0.03) 0.22(0.02)
0.32(0.02) 0.22 (0.02) 0.21(0.02) 0.20(0.02) ${\Delta_{2}}=2.2{\Delta_{0}}$
90/120 0.17(0.02) 0.09(0.02) 0.07(0.02) 0.06(0.02) 0.14(0.03) 0.07(0.02)
0.07(0.02) 0.07(0.01) 150/150 0.20(0.02) 0.11(0.04) 0.08(0.02) 0.08(0.01)
0.16(0.02) 0.07(0.01) 0.07(0.02) 0.06(0.01) 240/180 0.29 (0.01) 0.21(0.06)
0.19(0.05) 0.17(0.05) 0.29(0.02) 0.14 (0.04) 0.13(0.03) 0.12(0.02)

Table 1: Misclassification rate over 30 repeated experiments from the
individual learning scheme (Ind.) and data integration scheme (Prop.), where
the data integration utilizes $M$ datasets. The misclassification rate
reported in Table 1 is the average over 6 different datasets.

Table 1 shows that the proposed data integration method can substantially
reduce the misclassification rate for both supervised learning classifier QDA
and semi-supervised learning classifier LP by $9\%$ to $62\%$ over different
simulation schemes. Notice that the improvement from the data integration
increases as the number of incorporated datasets and the size of the training
set within $N_{{train}}$ increase. This improvement pattern implies that the
proposed alignment scheme captures the underlying common conditional
distribution of label $Y$ given covariate $X$, which boosts the effective
sample size of training data on each dataset. In addition, the improvement of
our method is more significant when the separability of $X$ across different
classes is strong, implying that the improvement also depends on the signal
intensity of the shared response-covariate patterns within each individual
dataset. In addition, we perform an additional simulation study where the
labels are unbalanced in training sets. The results show that the proposed
method achieves lower misclassification rate compared with the individual
learning method. The detailed simulation setup and results are in the Section
12 of Supplemental Materials.

###  4.2 Simulations for temporal integration for multiple datasets

In this subsection, we consider simulation studies with the two-dimensional
covariates dynamically changing over time. We compare the classification
performance of integrating the temporal data procedure in Section 3.2 with
that of the individual learning approach. For individual learning, we
independently train a classifier on the training data, and then make a
prediction on the testing data from each dataset at each time point. For our
temporal integration scheme, we train the classifier on the common latent
space integrating aligned training data from different datasets and time
points. We then calculate the misclassification rate averaged from the six
temporal datasets using these two approaches.

Specifically, we consider data points with three classes $Y=k(k=1,2,3)$. For
each dataset at time point $t=1$, we generate the covariates in the same way
as Section 4.1, and use $X^{(1)}$ to denote the observed $X$ at time $t=1$. To
mimic the temporal pattern in $X$ over time points $t=1,\cdots,T$, we propose
a sequential update on $X^{(t)}$ based on $X^{(t-1)}$ as
$X^{(t)}=\bm{b}+\bm{A}X^{(t-1)}$, where $\bm{b}$ introduces a global drift of
$X$ from time point $t-1$ to $t$, and $\bm{A}$ provides a transformation on
$X^{(t-1)}$. The stochastic drift term is $\bm{b}=(b_{1},b_{2})$, where
$b_{1},b_{2}$ are sampled from uniform distribution $\text{Unif}(0,0.4)$. The
transformation matrix $\bm{A}$:

$\displaystyle\bm{A}=\begin{pmatrix}1+\lambda&0\\\
0&1+\lambda\end{pmatrix}\begin{pmatrix}\text{cos}\theta&-\text{sin}\theta\\\
\text{sin}\theta&\text{cos}\theta\end{pmatrix},$

where $\lambda$ and $\theta$ are randomly sampled from $\text{Unif}(-0.2,0.2)$
and $\text{Unif}(-0.2\pi,0.2\pi)$, respectively. Notice that the labels of the
data points do not change during the process. Therefore, we construct a data
sequence $(X^{(1)},X^{(2)},\cdots,X^{(T)},Y)$ for each dataset. Then we split
data points from each dataset into $30$ training samples and $90$ testing
samples following the same procedure discussed in Section 4.1. In the
following, we perform numerical experiments under different lengths of time
$T=4,8,12,20$, using different numbers of integrated datasets $M=1,3,5$.

Length of time Num of datasets Classifier LP QDA Ind. Prop. Ind. Prop. $T=4$
$M=1$ 0.146(0.02) 0.141(0.02) 0.156(0.01) 0.179(0.01) $M=3$ 0.176(0.01)
0.133(0.01) 0.184(0.01) 0.139(0.01) $M=5$ 0.169 (0.01) 0.125(0.01) 0.174(0.01)
0.118(0.01) $T=8$ $M=1$ 0.142(0.02) 0.145(0.02) 0.156(0.01) 0.153(0.01) $M=3$
0.182(0.01) 0.128(0.01) 0.184(0.01) 0.121(0.01) $M=5$ 0.172 (0.02) 0.121(0.01)
0.174(0.01) 0.113(0.01) $T=12$ $M=1$ 0.142(0.01) 0.145(0.02) 0.156(0.01)
0.137(0.01) $M=3$ 0.188(0.01) 0.135(0.01) 0.184(0.01) 0.115(0.01) $M=5$ 0.176
(0.01) 0.127(0.01) 0.174(0.01) 0.108(0.01) $T=20$ $M=1$ 0.132(0.00)
0.064(0.00) 0.087(0.02) 0.084(0.01) $M=3$ 0.161(0.00) 0.108(0.01) 0.117(0.02)
0.093(0.01) $M=5$ 0.157 (0.00) 0.104(0.01) 0.183(0.01) 0.095(0.01)

Table 2: Misclassification rates from the individual learning approach (Ind.)
and the proposed temporal data integration scheme (Prop.) under different
lengths of time $T$, where the proposed method utilizes $M$ datasets. QDA and
LP stand for quadratic discriminative analysis and label propagation
algorithm, respectively.

Table 2 shows that the proposed temporal data integration method can reduce
the misclassification rate for both classifiers QDA and LP compared with the
individual learning even when the length of time is short $(T=4)$ or the
number of incorporated datasets is small $(M=1)$. These results demonstrate
the effectiveness of temporal alignment in reducing the systematic
heterogeneity even with a small sample size. For both the QDA and LP
classifiers, the improvement from temporal integration increases as $T$
increases given that $M$ is fixed. On the other hand, the improvement also
increases as $M$ increases given that $T$ is fixed. In addition, when $T=4$,
the classification improvement from our method over the individual learning
method increases about $2.5\%$ when $M$ increases from 3 to 5. On the other
hand, when $T=20$, the improvement increases about $13\%$ when $M$ increases
from 3 to 5 on averaging two classifiers. This result implies that the
proposed method can utilize common temporal patterns on the conditional
distribution of response to further enhance the power of detecting the
patterns shared across different datasets.

Notice that by applying the proposed method in Section 3.2, we obtain the
sequential latent space to capture the common dynamic patterns of observed
covariates and response-covariate association. To measure the performance of
temporal alignment over all time windows, we investigate the misclassification
rates of each time point, and calculate the average misclassification rates on
the trajectory. A lower average misclassification rate indicates that the
aligned temporal pattern can better capture the overall dynamic of response-
covariate association.

###  4.3 Simulations for Synthetic Neuronal Spike Data

In this subsection, we consider the simulation experiments to mimic the
application of the proposed data integration method to real neuronal spike
data, where we simulate the neuronal spike data by following the widely used
Poisson generalized linear modeling and Latent variable modeling [20].
Specifically, we consider total $M=6$ datasets, where each consists of spike
counts over $T=20$ time points. The neuronal spikes of the $m$th dataset at
the $t$th time point are governed by dynamic low-dimensional latent status
$X_{t}^{(m)}\in R^{k},m=1,\cdots,6,t=1,\cdots,20$, where $k$ denotes the
number of latent factors. The spike counts for the $j$th neuron at the $i$th
trial are generated as follows: $S_{j}^{(i)}(t)\mid
X_{t}^{(m)}\sim\text{Poiss}\big{(}\sigma(W_{ij}X_{t}^{(m)})\big{)},i=1,\cdots,n,\;j=1,\cdots,p$,
where $W_{ij}\in R^{p}$ denotes the loading vector, and
$\sigma(\cdot)=\frac{\exp(\cdot)}{1+\exp(\cdot)}$ is the link function. Here
we consider the number of trials $n=60$ and number of neurons $p=60$. The
latent status $\\{X_{t}^{(m)}\\}$ are sampled from mixture of Gaussian
distributions with multiple clusters corresponding to different external
stimuli. The detailed generation procedure for $\\{X_{t}^{(m)}\\}$ is the same
to the setup in Section 4.2.

We perform stimuli classification based on the labeled spike data
$\\{S^{(i)}(t),Y^{(i)}\\}_{i=1}^{n}$. Specifically, we follow the proposed
integrated latent alignment framework by first utilizing an autoencoder to
achieve dimension reduction of the spike data $\\{S^{(i)}(t)\\}$ for each
dataset in a sliding time interval with four time points. Then classifiers are
trained on the integrated data points, and are compared with the individual
learning method. We consider the proposed method integrating spike data from
different numbers of datasets and time points, and Table 3 illustrates the
average misclassification rates where Ind. denotes the individual learning
approach, and Prop. denotes the proposed scheme. Table 3 shows that the
proposed temporal data integration method can reduce the misclassification
errors for both classifiers QDA and LP compared with the independent training
scheme. For both classifiers, the improvement from the proposed method
increases as $T$ or $M$ increases, and the improvement maximizes at about
$25\%$ when $T=16$ and $M=6$. The simulation results demonstrate the
effectiveness of the temporal alignment in reducing the systematic
heterogeneity in subject-wise spike data.

Length of time Num of datasets Classifier LP QDA Ind. Prop. Ind. Prop. $T=4$
$M=2$ 0.403(0.02) 0.413(0.02) 0.484(0.01) 0.465(0.01) $M=4$ 0.411(0.01)
0.365(0.03) 0.463(0.01) 0.435(0.02) $M=6$ 0.459(0.01) 0.408(0.05) 0.466(0.01)
0.420(0.02) $T=10$ $M=2$ 0.383(0.02) 0.371(0.02) 0.401(0.01) 0.381(0.03) $M=4$
0.410(0.02) 0.307(0.01) 0.414(0.01) 0.325(0.03) $M=6$ 0.437(0.02) 0.319(0.03)
0.428(0.01) 0.321(0.03) $T=16$ $M=2$ 0.378(0.01) 0.359(0.01) 0.376(0.03)
0.346(0.03) $M=4$ 0.385(0.01) 0.299(0.02) 0.392(0.01) 0.283(0.03) $M=6$
0.411(0.01) 0.310(0.02) 0.406(0.01) 0.294(0.03)

Table 3: Misclassification rates from the individual learning approach (Ind.)
and the proposed temporal data integration method (Prop.) under different
number of time points $T$, where the proposed method integrates $M$ datasets.
Here QDA and LP stand for quadratic discriminative analysis and label
propagation algorithm, respectively. The spike data are simulated based on
Poisson generalized linear modeling and latent variable modeling.

##  5 Real Data Application

In this section, we apply the proposed ILA framework to the rodent
electrophysiological study introduced in Section 2. Here we use spike data
from five rodents recorded as the rodents smelled different odors in the
sequence memory task. The spike data from the $l$th rat is reformulated as a
three-dimensional data array $\mathcal{A}^{(l)}=\\{A_{ijk}\\},l=1,2,3,4,5$
where the $i,j,k$ represent the trials, neurons, and time points for each
rodent, respectively. The spike data $\mathcal{A}^{(l)}$ consists of binary
elements where $A_{ijk}=1$ indicates that the specific neuron fired, and
$A_{ijk}=0$ that it did not fire. In each trial, the subject samples one of
five odors labeled "A", "B", "C", "D", "E". We then extract 4 seconds of spike
data per trial (2 s before and after odor onset) in bins of 10 ms, resulting
in 400 time points for each trial. In the following, we investigate a subset
of the original spike data from trials corresponding to the four labels "A",
"B", "C", "D" (odor E was excluded because of its lower number of trials) and
time interval between the $250$ millisecond and $500$ millisecond after the
onset of odors. This time interval of the data was shown to exhibit the
strongest differentiation of odor information, as well as significant
differentiation of neuronal coding for different temporal order and trial
outcome information ([29]). In summary, the dimensions corresponding to
trials, neurons, and time points of spike data for the five rodents are
$168\times 46\times 25$, $152\times 49\times 25$, $207\times 104\times 25$,
$133\times 92\times 25$, and $174\times 79\times 25$, respectively.

In the following, we perform the first stage in the ILA framework to compress
the original temporal spike data into a low-dimensional latent space for each
rodent by using autoencoder, which is a nonlinear dimensionality reduction
method based on neural networks. The implementation details of the autoencoder
are in the Section 10 of the Supplemental Materials. Through the above
process, the spike data from each $100$ms sub-window are compressed into two-
dimensional latent features. We collect the preprocessed spike data from the
16 sub-windows for the $s$th rodent $(s=1,2,3,4,5)$ as
$(X^{(s)}(1),X^{(s)}(2),\cdots,X^{(s)}(16),Y^{(s)})$, where $X^{(s)}(t)\in
R^{n_{l}\times 2}$ denotes the latent neural representation of the $s$th rat
at the $t$th sub-windows $(t=1,\cdots,16)$ and $n_{l}$ is the number of trials
from the $l$th rat. In addition, $Y^{(l)}\in\\{1,2,3,4\\}^{n_{l}}$ denotes the
odor labels "A", "B", "C", "D" corresponding to each trial of the $l$th rat.
In the following, we first apply the proposed integrated supervised learning
method in Section 3.1 to the preprocessed spike data, and compare the results
with that of individual learning on the odor classification. Specifically, we
sample $20\%$ of the trials from each rat as training data, and the remaining
$80\%$ of trials serve as testing data. Then we train a classifier by our
integration method and the individual learning method separately, and
investigate the corresponding misclassification rates for the two methods on
each sub-window. We use the K-nearest neighbors algorithm (KNN) as our
classifier following the suggestion in [29] due to the nonlinear boundary
among the latent neural representations from different odors. The number of
neighbors for KNN is set as 25 labeled data points for both the individual
learning method and our method; i.e., the label of a specific data point is
determined by the labels of the 25 nearest data points, where the number of
neighbourhoods follows the setup in [29].

Figure 3 illustrates the classification accuracy of KNN trained at each
$100ms$ sub-windows by the individual learning method (blue line) and the
proposed integration learning method (red line), respectively. The odor
classification analysis is performed for each 100-ms time window based on the
output from the data compression stage, i.e., the autoencoder. Specifically,
the autoencoder compresses the original neural spike data within each 100-ms
time window into two-dimensional latent features, and the 100-ms time window
slides from 250 ms to 500 ms with a step of 10 ms. The results show that the
proposed method achieves higher odor classification accuracy for all five
rodents over the sequential sub-windows. The improvement from the proposed
data integration suggests the neuronal ensembles from the five rodents have in
common the fact that they exhibit a differentiation of odor information in
that time window. More specifically, it shows that our approach can detect the
common neuronal coding among the different subjects, and utilizes the shared
pattern to enhance the performance of odor classification. In addition, the
results show significant improvements from the proposed method for rat 1, rat
4, and rat 5 across all time intervals. For rat 2 and rat 3, the improvements
are more significant in the second half of the time intervals, while
improvements in the first half of the time intervals are modest. We also
investigate the improvement of the average classification accuracy among the
five rats, which is illustrated in the bottom right of Figure 3. The result
shows that the data integration method achieves about $10\%$ higher average
classification accuracy compared with the individual learning method. For our
specific application, this is a substantial improvement.

Figure 3: The four-odor classification accuracy based on the three rodents’
latent neural features extracted by autoencoder within time interval $250\sim
500$ milliseconds. The blue line denotes the classification accuracy from KNN
trained by the individual learning method at each time point. The red line
denotes KNN’s classification accuracy trained by the proposed data integration
method in Section 3.1.

In the following, we investigate the odor classification performance of the
temporal integration method introduced in Section 3.2, and compare the
performance of the individual learning approach based on the five rodents
spike samples within the time interval $250\sim 500$ milliseconds. For each
rodent, we sample the odor labels and the corresponding latent neuronal coding
from $20\%$ of the trials as training data, and the remaining $80\%$ of the
trials as testing data. Then we train the KNN classifier by our temporal
integration method and individual learning method separately, and examine the
misclassification rates of the two methods on each 16 sub-windows. In
addition, the number of template samples in the aligned space is set as
$n=200$ for the temporal alignment method.

We use the individual learning scheme as the baseline, and the classification
comparison is illustrated in Figure 4. The results show about $10\%$ accuracy
improvements from the proposed method for rodent 2 and rodent 4 across all
time intervals. For rodent 1, rodent 3, and rodent 5, the accuracy of the
proposed method is higher than the individual learning method at most of the
time points, and the improvements can be about $10\%$ at some time windows. In
addition, compared with the individual learning method, the accuracy of our
method changes more smoothly over the time because the proposed method
utilizes the time dependency among latent representations of spike samples.
This is another advantage of the proposed method in scientific interpretation
in addition to improving classification accuracy. Notice that the non-temporal
data alignment method introduced in Section 3.1 is designed for response
classification at each time point, while the temporal alignment method
introduced in Section 3.2 aims to capture the common dynamic patterns of
covariate-response association better and more smoothly. In other words, the
temporal alignment method in Section 3.2 does not necessarily achieve better
timepoint-wise response classification than the non-temporal method, and vice
versa.

Figure 4: The four-odor classification accuracy within time interval $250\sim
500$ milliseconds based on rodents’ electrophysiological data. The blue line
denotes the classification accuracy from KNN trained by the individual
learning method at each time point. The red line denotes KNN’s classification
accuracy trained by the proposed temporal integration method in Section 3.2.

Notice that Figure 3 and 4 correspond to the classification-oriented alignment
method presented in Section 3.1, and the covariate-response association
dynamic alignment method presented in Section 3.2. The non-temporal data
alignment method introduced in Section 3.1 is designed for response
classification at each time point, while the temporal alignment method
introduced in Section 3.2 aims to capture the common dynamic patterns of
covariate-response association better and more smoothly. In other words, the
temporal alignment method in Section 3.2 does not necessarily achieve better
timepoint-wise response classification than the non-temporal method and vice
versa. Therefore, instead of comparing alignment methods in Section 3.1 and
Section 3.2, Figure 3 and Figure 4 compare each of the proposed methods with
the individual classification scheme, and demonstrate that both of the
proposed methods can improve task-relevant response classification over the
individual method. In contrast, the advantage of temporal alignment method in
Section 3.2 is to allow one to investigate the underlying common covariate-
response association that smoothly changes over time, which is illustrated in
following Figure 5.

In summary, the proposed alignment approach can integrate information from
both different subjects and different time windows. The improvement from the
temporal data integration on the five subjects implies the existence of a
common structure in the neuronal coding associated with the processing of the
different odors. More importantly, the dynamic of the neuronal coding
clustering is also shared across subjects to some extent. This finding is
consistent with prior evidence that hippocampal neurons represent information
about task events critical for accurate performance [3, 16].

We also utilize the temporal data integration in Section 3.2 to investigate
the common dynamic pattern of neuronal coding among the five rodents.
Specifically, we first obtain the rodent-wise alignment mapping
$\\{\Omega^{(s)}\\}_{s=1}^{5}$ and temporal similarity matrices
$\\{D(t)\\}_{t=1}^{16}$ of 16 sub-windows from Stage 2 in the temporal
integration where we use the true odor labels to construct the label distance
matrices $\\{M^{(s)}(t)\\},s=1,\cdots,5,\;t=1,\cdots,16$. Then we reconstruct
the aligned latent space of neuronal coding for the five rodent at the
$250\sim 500$ms time windows by following Stage 3 based on $D(t)$ and
$\Omega^{(s)}$. Finally, we fix the estimated cross-subject latent space from
the time window at $250$ to $500$ms, and then project neuron spike samples
from the other time windows throughout $0$ to $1500$ms into this latent space.
Figure 5 shows neural activity trajectories for all 5 rodents throughout the
$0$ to $1500$ms period. Notice that we investigated the odor misclassification
rates by using common latent spaces at different dimensions $k=2,3,4$, and
found that different dimensions provide very close rates. For the purpose of
visualization, we adopt $k=2$ to illustrate the common dynamic patterns of the
latent neural coding.

Figure 5: Odor-specific neuronal coding and dynamics in the aligned latent
space. a. Differentiation of the neural activity across odor presentations on
the aligned latent dimensions (odors A-D; data aggregated across all 5
subjects). Each dot indicates the location of the neural activity for a 100-ms
bin from a given subject, with the color denoting the odor label of the
corresponding trial. b. Neural activity trajectory (mean across subjects) on
the latent dimensions during each odor trial type. c. Neural activity
trajectories for each subject showing comparable patterns across subjects.
Ellipses are consistent across panels and represent the spread of the latent
coding for each odor type, specifically the two largest square-rooted singular
values from the covariance of latent coding within each cluster.

The results show that the neural activity on odor A trials is very distinct
than that on odors B, C, and D trials. This can be explained by the fact that
the sequence always starts with odor A (so animals are not required to
remember it), and it is preceded by a running period (so that activity
dynamics associated with running extend into the odor period). Therefore, we
focus on comparing the coding properties and dynamics observed during odors
B-D, in which the behavioral and cognitive demands are matched.

Neural trajectories have been reported in the motor system in the past [23].
In fact, the motor system was an ideal system to initially demonstrate such
trajectories (or proof of concept) because of its simplicity due to a high
correspondence between motor system and direction of movement. However,
capturing neural trajectories for cognitive information is this paper, is
considerably more difficult – neural activity in the hippocampus is several
synapses (connections) away from the external world, and thus represents
information that is more abstract, multi-dimensional and variable across
subjects.

Notice that existing literature has utilized dimension-reduction method to
visualize trajectory representing the response over time to each odor
presentation [9]. However, our data integration approach allows, for the first
time, integration of such trajectories from different subjects into a common
(aligned) latent space (Fig 5b).

Specifically, the activity on odors B, C, and D trials starts at a similar
location in the latent space and heads toward the center of their respective
clusters (0-250 ms period), reaching maximum differentiation in the 250-500 ms
period. Toward the end of the odor presentations (750-1,250 ms period), the
activity for each odor then converges back toward the beginning of their
trajectory (faster for odor D, than C, than B), before reaching a distinct
region of the latent space in the post-trial activity (1,500ms; indicated by
arrowheads). Our data integration approach also allows the visualization of
each subject’s neural activity in the aligned latent space (Fig 5c), which
reveals comparable trajectories across the 5 subjects. Critically, these
results would have been difficult, if not impossible, to interpret without
this data integration approach, as the neural activity of each subject would
otherwise be in distinct latent spaces that could not be directly compared
across subjects.

##  6 Discussion

In this paper, we propose a novel and general data integration framework based
on optimal transport for supervised learning tasks, which is motivated by a
rodent electrophysiological experiment investigating non-spatial memory
mechanisms. Our method is effective when the data domains and dimensionalities
of each dataset are different, and the integration is effective even if the
number of datasets is small.

We apply the proposed framework to electrophysiological data collected from a
rodent non-spatial memory study. The proposed method identifies the existence
of common neuronal coding patterns across subjects during the processing of
different odors. This discovery is consistent with previous findings that
hippocampal neurons represent information about task events critical for
accurate performance [3, 16]. Furthermore, our results suggest that
heterogeneous neuron datasets can be systematically integrated to increase the
discriminating power of an electrophysiological experiment, which may have
broad implications in neuroscience research and beyond.

## References

  * Agueh and Carlier [2011] Agueh, M. and Carlier, G. (2011).  Barycenters in the Wasserstein space.  SIAM Journal on Mathematical Analysis, 43(2):904–924.
  * Akil et al. [2011] Akil, H., Martone, M. E., and Van Essen, D. C. (2011).  Challenges and opportunities in mining neuroscience data.  Science, 331(6018):708–712.
  * Allen and Fortin [2013] Allen, T. A. and Fortin, N. J. (2013).  The evolution of episodic memory.  Proceedings of the National Academy of Sciences, 110(2):10379–10386.
  * Allen et al. [2014] Allen, T. A., Morris, A. M., Mattfeld, A. T., Stark, C. E., and Fortin, N. J. (2014).  A sequence of events model of episodic memory shows parallels in rats and humans.  Hippocampus, 24(10):1178–1188.
  * Almotiri et al. [2017] Almotiri, J., Elleithy, K., and Elleithy, A. (2017).  Comparison of autoencoder and principal component analysis followed by neural network for e-learning using handwritten recognition.  In 2017 IEEE Long Island Systems, Applications and Technology Conference (LISAT), pages 1–5. IEEE.
  * Bengio et al. [2013] Bengio, Y., Courville, A., and Vincent, P. (2013).  Representation learning: A review and new perspectives.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(8):1798–1828.
  * Bengio et al. [2009] Bengio, Y. et al. (2009).  Learning deep architectures for AI.  Foundations and Trends in Machine Learning, 2(1):1–127.
  * Berry and Tkačik [2020] Berry, M. J. and Tkačik, G. (2020).  Clustering of neural activity: A design principle for population codes.  Frontiers in Computational Neuroscience, 14:20.
  * Brown et al. [2005] Brown, S. L., Joseph, J., and Stopfer, M. (2005).  Encoding a temporally structured stimulus with a temporally structured neural representation.  Nature neuroscience, 8(11):1568–1576.
  * Burbank [2015] Burbank, K. S. (2015).  Mirrored stdp implements autoencoder learning in a network of spiking neurons.  PLoS computational biology, 11(12):e1004566.
  * Carroll et al. [2009] Carroll, M. K., Cecchi, G. A., Rish, I., Garg, R., and Rao, A. R. (2009).  Prediction and interpretation of distributed neural activity with sparse models.  NeuroImage, 44(1):112–122.
  * Churchland and Shenoy [2007] Churchland, M. M. and Shenoy, K. V. (2007).  Temporal complexity and heterogeneity of single-neuron activity in premotor and motor cortex.  Journal of Neurophysiology, 97(6):4235–4257.
  * Courty et al. [2017] Courty, N., Flamary, R., Habrard, A., and Rakotomamonjy, A. (2017).  Joint distribution optimal transportation for domain adaptation.  Advances in Neural Information Processing Systems, 30.
  * Courty et al. [2014] Courty, N., Flamary, R., and Tuia, D. (2014).  Domain adaptation with regularized optimal transport.  In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pages 274–289. Springer.
  * Ebrahimi et al. [2020] Ebrahimi, M., Calarco, N., Campbell, K., Hawco, C., Voineskos, A., and Khisti, A. (2020).  Time-resolved fMRI shared response model using Gaussian process factor analysis.  arXiv preprint arXiv:2006.05572.
  * Eichenbaum [2014] Eichenbaum, H. (2014).  Time cells in the hippocampus: A new dimension for mapping memories.  Nature Reviews Neuroscience, 15(11):732–744.
  * Hackam and Redelmeier [2006] Hackam, D. G. and Redelmeier, D. A. (2006).  Translation of research evidence from animals to humans.  JAMA, 296(14):1727–1732.
  * Hinton and Salakhutdinov [2006] Hinton, G. E. and Salakhutdinov, R. R. (2006).  Reducing the dimensionality of data with neural networks.  Science, 313(5786):504–507.
  * Huang and Toyoizumi [2016] Huang, H. and Toyoizumi, T. (2016).  Clustering of neural code words revealed by a first-order phase transition.  Physical Review E, 93(6):062416.
  * Keeley et al. [2020] Keeley, S. L., Zoltowski, D. M., Aoi, M. C., and Pillow, J. W. (2020).  Modeling statistical dependencies in multi-region spike train data.  Current Opinion in Neurobiology, 65:194–202.
  * Kingma et al. [2019] Kingma, D. P., Welling, M., et al. (2019).  An introduction to variational autoencoders.  Foundations and Trends in Machine Learning, 12(4):307–392.
  * Kleindienst et al. [2011] Kleindienst, T., Winnubst, J., Roth-Alpermann, C., Bonhoeffer, T., and Lohmann, C. (2011).  Activity-dependent clustering of functional synaptic inputs on developing hippocampal dendrites.  Neuron, 72(6):1012–1024.
  * Lindén et al. [2022] Lindén, H., Petersen, P. C., Vestergaard, M., and Berg, R. W. (2022).  Movement is governed by rotational neural dynamics in spinal motor networks.  Nature, 610(7932):526–531.
  * Marinelli and McCutcheon [2014] Marinelli, M. and McCutcheon, J. E. (2014).  Heterogeneity of dopamine neuron activity across traits and states.  Neuroscience, 282:176–197.
  * Mémoli [2014] Mémoli, F. (2014).  The Gromov–Wasserstein distance: A brief overview.  Axioms, 3(3):335–341.
  * Perel et al. [2007] Perel, P., Roberts, I., Sena, E., Wheble, P., Briscoe, C., Sandercock, P., Macleod, M., Mignini, L. E., Jayaram, P., and Khan, K. S. (2007).  Comparison of treatment effects between animal experiments and clinical trials: Systematic review.  BMJ, 334(7586):197.
  * Peyré et al. [2019] Peyré, G., Cuturi, M., et al. (2019).  Computational optimal transport: With applications to data science.  Foundations and Trends in Machine Learning, 11(5-6):355–607.
  * Sejnowski et al. [2014] Sejnowski, T. J., Churchland, P. S., and Movshon, J. A. (2014).  Putting big data to good use in neuroscience.  Nature Neuroscience, 17(11):1440–1441.
  * Shahbaba et al. [2022] Shahbaba, B., Li, L., Agostinelli, F., Saraf, M., Cooper, K. W., Haghverdian, D., Elias, G. A., Baldi, P., and Fortin, N. J. (2022).  Hippocampal ensembles represent sequential relationships among an extended sequence of nonspatial events.  Nature Communications, 13(1):1–17.
  * Speiser et al. [2017] Speiser, A., Yan, J., Archer, E. W., Buesing, L., Turaga, S. C., and Macke, J. H. (2017).  Fast amortized inference of neural activity from calcium imaging data with variational autoencoders.  Advances in neural information processing systems, 30.
  * Titouan et al. [2019] Titouan, V., Courty, N., Tavenard, R., and Flamary, R. (2019).  Optimal transport for structured data with application on graphs.  In International Conference on Machine Learning, pages 6275–6284. PMLR, 2019.
  * Vayer et al. [2020] Vayer, T., Chapel, L., Flamary, R., Tavenard, R., and Courty, N. (2020).  Fused gromov-wasserstein distance for structured objects.  Algorithms, 13(9):212.
  * Villani [2009] Villani, C. (2009).  Optimal transport: Old and new, volume 338.  Springer.
  * Wang et al. [2016] Wang, Y., Yao, H., and Zhao, S. (2016).  Auto-encoder based dimensionality reduction.  Neurocomputing, 184:232–242.
  * Wu et al. [2018] Wu, T., Zhao, W., Keefer, E., and Yang, Z. (2018).  Deep compressive autoencoder for action potential compression in large-scale neural recording.  Journal of neural engineering, 15(6):066019.
  * Zhang et al. [2018] Zhang, J., Nguyen, T., Cogill, S., Bhatti, A., Luo, L., Yang, S., and Nahavandi, S. (2018).  A review on cluster estimation methods and their application to neural spike data.  Journal of Neural Engineering, 15(3):031003.
  * Zhang et al. [2017] Zhang, Q., Borst, J. P., Kass, R. E., and Anderson, J. R. (2017).  Inter-subject alignment of MEG datasets in a common representational space.  Technical report, Wiley Online Library.