* [83] Victor-Bogdan Popescu, Jose Angel Sanchez-Martin, Daniela Schacherer, Sadra Safadoust, Negin Majidi, Andrei Andronescu, Alexandru Nedea, Diana Ion, Eduard Mititelu, Eugen Czeizler, and Ion Petre.  Netcontrol4biomed: a platform for biomedical network synthesis and analysis.  under review, 2021.
  * [84] Jun Ren, Jianxin Wang, Min Li, and Fangxiang Wu.  Discovering essential proteins based on PPI network and protein complex.  International journal of data mining and bioinformatics, 12 (1):24–43, 2015.
  * [85] Yannick Rochat.  Closeness centrality extended to unconnected graphs: The harmonic centrality index.  In Applications of Social Network Analysis, ASNA2009, 2009.
  * [86] Murielle Roussel, Philippe Moreau, Benjamin Hebraud, Kamel Laribi, Arnaud Jaccard, Mamoun Dib, Borhane Slama, Véronique Dorvaux, Bruno Royer, Laurent Frenzel, Sonja Zweegman, Saskia K Klein, Annemiek Broijl, Kon-Siong Jie, Jianping Wang, Veronique Vanquickelberghe, Carla de Boer, Tobias Kampfenkel, Katharine S Gries, John Fastenau, and Pieter Sonneveld.  Bortezomib, thalidomide, and dexamethasone with or without daratumumab for transplantation-eligible patients with newly diagnosed multiple myeloma (cassiopeia): health-related quality of life outcomes of a randomised, open-label, phase 3 trial.  The Lancet Haematology, 7(12):e874–e883, 2020.
  * [87] Britta Ruhnau.  Eigenvector-centrality — a node-centrality?  Social Networks, 22(4):357–365, 2000.
  * [88] Gert Sabidussi.  The centrality index of a graph.  Psychometrika, 31(4):581–603, 1966.
  * [89] R. Sahoo, T.S. Rani, and S.D. Bhavani.  Chapter 17 - Differentiating Cancer From Normal Protein-Protein Interactions Through Network Analysis.  In Quoc Nam Tran and Hamid R. Arabnia, editors, Emerging Trends in Applications and Infrastructures for Computational Biology, Bioinformatics, and Systems Biology, Emerging Trends in Computer Science and Applied Computing, pages 253–269. Morgan Kaufmann, Boston, 2016.
  * [90] Rafael Santos Coelho, Phablo Moura, and Yoshiko Wakabayashi.  The k-hop connected dominating set problem: approximation and hardness.  Journal of Combinatorial Optimization, 34, 11 2017.
  * [91] Mansoor Saqi, Johann Pellet, Irina Roznovat, Alexander Mazein, Stéphane Ballereau, Bertrand De Meulder, and Charles Auffray.  Systems medicine: The future of medical genomics, healthcare, and wellness.  Methods Mol Biol, 1386:43–60, 2016.
  * [92] P. Shannon, A. Markiel, O. Ozier, N. S. Baliga, J. T. Wang, D. Ramage, N. Amin, B. Schwikowski, and T. Ideker.  Cytoscape: a software environment for integrated models of biomolecular interaction networks.  Genome research, 13:2498–2504, 2003.
  * [93] P. Shannon, A. Markiel, O. Ozier, N. S. Baliga, J. T. Wang, D. Ramage, N. Amin, B. Schwikowski, and T. Ideker.  Cytoscape: An open source platform for complex network analysis and visualization.  https://cytoscape.org/, 2021.  [Online; accessed 4-February-2021].
  * [94] Marvin E. Shaw.  Group structure and the behavior of individuals in small groups.  The Journal of Psychology, 38(1):139–149, 1954.
  * [95] R. W. Shields and J. B. Pearson.  Structural controllability of multi-input linear systems.  IEEE Transactions on Automatic Control, 21(2):203–212, 1976.
  * [96] Q Tian, N D Price, and L Hood.  Systems cancer medicine: towards realization of predictive, preventive, personalized and participatory (p4) medicine.  J Intern Med, 271(2):111–121, Feb 2012.
  * [97] Rodger E. Tiedemann, Yuan Xao Zhu, Jessica Schmidt, Chang Xin Shi, and et al.  Identification of molecular vulnerabilities in human multiple myeloma cells by rna interference lethality screening of the druggable genome.  Cancer Research, 72, 2012.
  * [98] Feng Wang, Hongwei Du, Erika Camacho, Kuai Xu, Wonjun Lee, Yan Shi, and Shan Shan.  On positive influence dominating sets in social networks.  Theoretical Computer Science, 412(3):265–269, 2011.  Combinatorial Optimization and Applications.
  * [99] R. Wester, B. van der Holt, E. Asselbergs, and et al.  Phase ii study of carfilzomib, thalidomide, and low-dose dexamethasone as induction and consolidation in newly diagnosed, transplant eligible patients with multiple myeloma; the carthadex trial.  Haematologica, 104(11):2265‐2273, 2019.
  * [100] David S Wishart, Yannick D Feunang, An C Guo, Ana Marcu, and et al.  Drugbank 5.0: a major update to the drugbank database for 2018.  Nucleic Acids Research, 46:1074–1082, 2017.
  * [101] Jie Wu, Mihaela Cardei, Fei Dai, and Shuhui Yang.  Extended dominating set and its applications in ad hoc networks using cooperative communication.  IEEE Trans. Parallel Distrib. Syst., 17(8):851–864, August 2006\.
  * [102] Stefan Wuchty.  Controllability in protein interaction networks.  Proceedings of the National Academy of Sciences, 111(19):7156–7160, 2014.
  * [103] Stefan Wuchty and Peter F. Stadler.  Centers of complex networks.  Journal of Theoretical Biology, 223(1):45–53, 2003.
  * [104] Rong Yang and Leyla Zhuhadar.  Extensions of closeness centrality?  In Proceedings of the 49th Annual Southeast Regional Conference, ACM-SE ’11, pages 304–305, New York, NY, USA, 2011. Association for Computing Machinery.
  * [105] Haiyuan Yu, Philip M Kim, Emmett Sprecher, Valery Trifonov, and Mark Gerstein.  The importance of bottlenecks in protein networks: correlation with gene essentiality and expression dynamics.  PLoS computational biology, 3(4):e59, April 2007.
  * [106] Jiguo Yu, Nannan Wang, Guanghui Wang, and Dongxiao Yu.  Connected dominating sets in wireless ad hoc and sensor networks – a comprehensive survey.  Computer Communications, 36(2):121–134, 2013.
  * [107] Peng Zhang, Lin Tao, Xian Zeng, Chu Qin, Shangying Chen, Feng Zhu, Zerong Li, Yuyang Jiang, Weiping Chen, and Yu-Zong Chen.  A protein network descriptor server and its use in studying protein, disease, metabolic and drug targeted networks.  Briefings in Bioinformatics, 18(6):1057–1070, 08 2016.
  * [108] Xiao-Fei Zhang, Le Ou-Yang, Yuan Zhu, Meng-Yun Wu, and Dao-Qing Dai.  Determining minimum set of driver nodes in protein-protein interaction networks.  BMC Bioinformatics, 16(1):146, May 2015.


# Latent Variable Double Gaussian Process Model for Decoding Complex Neural
Data ††thanks: This project is sponsored by the Defense Advanced Research
Projects Agency (DARPA) under cooperative agreement No. N660012324016. The
content of the information does not necessarily reflect the position or the
policy of the Government, and no official endorsement should be inferred.

Navid Ziaei  Department of Electrical
& Computer Eng.
Isfahan University of Technology
Department of Computer Science
Worcester Polytechnique Institute
Worcester, MA, USA
<EMAIL_ADDRESS>Joshua J. Stim  Department of Psychiatry
& Behavioral Sciences
University of Minnesota
Minneapolis, MN, USA
<EMAIL_ADDRESS>Melanie D. Goodman-Keiser  Department of Psychiatry
& Behavioral Sciences
University of Minnesota
Minneapolis, MN, USA
<EMAIL_ADDRESS>Scott Sponheim  Department of Psychiatry
& Behavioral Sciences
University of Minnesota
Minneapolis Veterans Affairs
Health Care System
Minneapolis, MN, USA
<EMAIL_ADDRESS>Alik S. Widge  Department of Psychiatry & Behavioral
Sciences
University of Minnesota
Minneapolis, MN, USA
<EMAIL_ADDRESS>Sasoun Krikorian  Department of Computer Science
Worcester Polytechnique Institute
Worcester, MA, USA
<EMAIL_ADDRESS>Ali Yousefi  Department of Computer Science
Worcester Polytechnique Institute
Worcester, MA, USA
<EMAIL_ADDRESS>
###### Abstract

Non-parametric models, such as Gaussian Processes (GP), show promising results
in the analysis of complex data. Their applications in neuroscience data have
recently gained traction. In this research, we introduce a novel neural
decoder model built upon GP models. The core idea is that two GPs generate
neural data and their associated labels using a set of low- dimensional latent
variables. Under this modeling assumption, the latent variables represent the
underlying manifold or essential features present in the neural data. When GPs
are trained, the latent variable can be inferred from neural data to decode
the labels with a high accuracy. We demonstrate an application of this decoder
model in a verbal memory experiment dataset and show that the decoder accuracy
in predicting stimulus significantly surpasses the state-of-the-art decoder
models. The preceding performance of this model highlights the importance of
utilizing non-parametric models in the analysis of neuroscience data.

###### Index Terms:

Latent Variable Model, Gaussian Process Model, Dimensionality Reduction,
Manifold Representation, Neural Decoder

##  I Introduction

With advancements in neural recording technology, we observe consistent growth
in the size and modality of neural data being recorded. In parallel, machine
learning (ML) tools and techniques have shown unprecedented growth and
advancement in the last decade. This parallel growth holds the promise for
further discoveries of brain function. Specifically, we anticipate the
emergence of new neural decoder models that achieve a higher level of
prediction accuracy in inferring underlying processes or factors influencing
neural activity. For example, we expect a more precise estimation of a rat’s
movement trajectory by observing the spiking activity of an ensemble of place
cells, or decoding arm movement observing motor cortex neural activity.

GPs are non-parametric models often used for both regression and
classification tasks [1]. These models exhibit a high expressive power;
furthermore, they draw a more precise picture of their estimations as they
provide a measure of uncertainty in their predictions. Extensions of GPs, such
as the latent variable GP, have been recently introduced for dimensionality
reduction, showing promising results in inferring low-dimensional manifolds
representing neural activity [2]. While these methods open the door to better
characterizing neural data, they also have their limitations. For instance,
most of these approaches draw a point estimate in their inference step, making
them prone to overfitting. Additionally, stimuli information is often
discarded in the estimation of the underlying manifold, limiting their
decoding accuracy.

In this research, we develop a novel latent variable GP model that
characterizes both neural activity and associated stimuli information [3]. Two
GPs, derived through a common latent variable, characterize neural activity
and associated stimuli (or labels). For the model, we derive a Bayesian
inference solution allowing us to have a robust inference of the underlying
manifold. We show that the inferred manifold captures essential features
present in neural recordings. We further show these features can be used in
predicting the data label or stimulus with high accuracy. The double GP and
its inference solution resemble an auto-encoder framework, where the input to
the encoder is neural activity, and the decoder generates both the label and
neural data. However, the fundamental difference lies in the framework’s
ability to run on datasets with a limited sample size. With our model, we can
find the proper dimension of the manifold, whilst this is not the case for the
auto-encoder models. More importantly, our proposed framework provides
probabilistic inference and prediction, leading to a more robust overall
performance.

##  II Latent Varaible Double GP

In this section, we first define the model. We then discuss model training,
inference, and finally the decoding processes.

###  II-A Model Definition

Let’s assume we observe neural data as continuous measurement defined by
$\mathbf{Y}^{n}\in\mathbb{R}^{N\times D}$, where $N$ is the number of data
points and $D$ represents the number of neural features. We also have
$\mathbf{Y}^{s}\in\mathbb{R}^{N\times K}$ , with $K$ as the number of possible
categories or stimuli associated with each data point. We further assume there
is a latent variable $\mathbf{X}\in\mathbb{R}^{N\times Q}$, where $Q\ll D$,
that captures underlying essential features deriving both $\mathbf{Y}^{n}$ and
$\mathbf{Y}^{s}$. Each row of $\mathbf{X}$ is assumed to be generated from a
multi-variate Gaussian distribution with mean $\mathbf{0}$ and covariance
$\mathbf{I}_{Q}$ and these latent variables are I.I.D; thus, their prior
distribution is defined by:

$p\left(\mathbf{X}\right)=\prod^{N}_{i=1}\mathcal{N}\left(\mathbf{x_{i}};\mathbf{0},\mathbf{I}_{Q}\right),$
(1)

where $\mathbf{x}_{i}$ is the vector of the latent variable for the $i$th
sample. We then assume that there are two GP prior on function values
$\mathbf{F}^{n}\in\mathbb{R}^{N\times D}$ and
$\mathbf{F}^{s}\in\mathbb{R}^{N\times K}$ defined by

$\displaystyle
p\left(\mathbf{F}^{n}\mid\mathbf{X}\right)=\prod^{D}_{d=1}\mathcal{N}\left(f^{n}_{:,d};\mathbf{0},\mathbf{K}^{n}_{NN}\right)$
(2) $\displaystyle
p\left(\mathbf{F}^{s}\mid\mathbf{X}\right)=\prod^{K}_{k=1}\mathcal{N}\left(f^{s}_{:,k};\mathbf{0},\mathbf{K}^{s}_{NN}\right),$
(3)

where $\mathbf{K}^{n}_{NN}$ and $\mathbf{K}^{s}_{NN}$ are the $N\times N$
covariance matrices defining functions behind the continuous and discrete
observations.
$\mathbf{f}^{n}_{:,d}\triangleq\left[f^{n}_{d}\left(\mathbf{x}_{1}\right),\ldots,f^{n}_{d}\left(\mathbf{x}_{N}\right)\right]$
and
$\mathbf{f}^{s}_{:,d}\triangleq\left[f^{s}_{d}\left(\mathbf{x}_{1}\right),\ldots,f^{s}_{d}\left(\mathbf{x}_{N}\right)\right]$
are function values representing the $d$th and $k$th columns of the functions.
Elements of covariance matrices are defined by a kernel function, with
$k^{n}_{i,j}=k_{\phi}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)$ for different
pairs of $i$ and $j$ indices. Similarly, we have
$k^{s}_{i,j}=k_{\theta}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)$ for
$\mathbf{K}^{s}_{NN}$ covariance. Here, $\phi$ and $\theta$ are the free
parameters of the kernel function. We finally assume:

$y_{i,d}=f^{n}_{d}\left(\mathbf{x_{i}}\right)+\epsilon_{d},$ (4)

where $\epsilon_{d}\sim\mathcal{N}\left(0,\sigma^{2}_{d}\right)$ and
$\mathbf{F}^{n}=\left[\mathbf{f}^{n}_{:,1},\ldots,\mathbf{f}^{n}_{:,D}\right]$
is the GP defined by $\mathbf{K}^{n}_{NN}$. For the discrete process, we have

$y_{i,k}=g\left(f^{s}_{k}\left(\mathbf{x}_{i}\right)\right),$ (5)

in which $g\left(\cdot\right)$ is representing conditional probability of the
$y_{i,k}$ defined by a squeezing function and
$\mathbf{F}^{s}=\left[\mathbf{f}^{s}_{:,1},\ldots,\mathbf{f}^{s}_{:,D}\right]$
is the GP defined by $\mathbf{K}^{s}_{NN}$.

$\mathbf{X}$$\mathbf{F}^{n}$$\mathbf{Y}^{n}$$\mathbf{F}^{s}$$\mathbf{Y}^{s}$$\phi$$\theta$$\sigma$
Figure 1:  Graphical Model Depicting Latent Variable Double GP: The latent
variable $\mathbf{X}$ with a Gaussian prior derives both continuous and
discrete covariance structures $\mathbf{K}^{n}_{NN}$ and $\mathbf{K}^{s}_{NN}$
for $\mathbf{F}^{n}$ and $\mathbf{F}^{s}$ value functions. These functions
then define the generative model creating the discrete and continuous
measurements.

Eq. (1) to (5) define the generative model describing how samples of
continuous and discrete observations are created as a function of the latent
variables. Model parameters are tuned to enhance the likelihood of observed
data, which will be elaborated in the following subsection on the training
algorithm.

###  II-B Model Training

For the training, we require to find the maximum likelihood estimate of kernel
parameters $\phi$, $\theta$, and measurement noise $\sigma_{d}$. In our model,
$\mathbf{X}$ is unobserved; if $\mathbf{X}$ were known, the conditional
distribution of observation could be defined by:

$\displaystyle\textstyle
p\left(\mathbf{Y}^{n},\mathbf{Y}^{s}\mid\mathbf{F}^{n},\mathbf{F}^{s},\mathbf{X}\right)=p\left(\mathbf{Y}^{n}\mid\mathbf{F}^{n},\mathbf{X}\right)p\left(\mathbf{Y}^{s}\mid\mathbf{F}^{s},\mathbf{X}\right)$
$\displaystyle
p\left(\mathbf{Y}^{n}\mid\mathbf{F}^{n},\mathbf{X}\right)=\textstyle\prod_{i=1}^{N}\textstyle\prod_{d=1}^{D}p\left(y_{i,d}\mid
f^{n}_{d}\left(\mathbf{x}_{i}\right),\mathbf{X}\right)$ (6) $\displaystyle
p\left(\mathbf{Y}^{s}\mid\mathbf{F}^{s},\mathbf{X}\right)=\textstyle\prod_{i=1}^{N}\textstyle\prod_{k=1}^{K}p\left(y_{i,k}\mid
f^{s}_{k}\left(\mathbf{x}_{i}\right),\mathbf{X}\right)$ (7)

The evidence for $\mathbf{Y}^{n}$, $\mathbf{Y}^{s}$ is defined by:

$p\left(\mathbf{Y}^{n},\mathbf{Y}^{s}\right)=\int
p\left(\mathbf{Y}^{n},\mathbf{Y}^{s}\mid\mathbf{F}^{n},\mathbf{F}^{s},\mathbf{X}\right)p\left(\mathbf{X}\right)\,d\mathbf{F}^{n}\,d\mathbf{F}^{s}\,d\mathbf{X},$
(8)

which needs to be maximized by adjusting $\phi$ and $\theta$. Solving this
integral is an interactable optimization problem; instead, we use doubly
stochastic variational inference method introduced in [4] to estimate the
parameters. In this approach, we introduce $\mathbf{Z}^{n}$ and
$\mathbf{Z}^{s}$ as sets of inducing points along with corresponding inducing
variables $\mathbf{U}^{n}$ and $\mathbf{U}^{s}$, each assumed to have a
Gaussian Process (GP) prior. The goal is to draw an approximate posterior
distribution over latent features
$p\left(\mathbf{X}\mid\mathbf{Y}^{n},\mathbf{Y}^{s}\right)$ using a
variational distribution $q_{\psi}\left(\mathbf{X}\right)$. To find this
approximate posterior, we approximate
$p\left(\mathbf{F}^{n},\mathbf{U}^{n}\mid\mathbf{Y}^{n},\mathbf{X}\right)$
with
$q\left(\mathbf{F}^{n},\mathbf{U}^{n}\right)=p\left(\mathbf{F}^{n}\mid\mathbf{U}^{n},\mathbf{X}\right)q_{\lambda}\left(\mathbf{U}^{n}\right)$,
and $p\left(\mathbf{F}^{s},\mathbf{U}^{s}\mid\mathbf{Y}^{s},\mathbf{X}\right)$
with
$q\left(\mathbf{F}^{s},\mathbf{U}^{s}\right)=p\left(\mathbf{F}^{s}\mid\mathbf{U}^{s},\mathbf{X}\right)q_{\gamma}\left(\mathbf{U}^{s}\right)$.
In this framework, $q_{\psi}\left(\mathbf{X}\right)$,
$q_{\lambda}\left(\mathbf{U}^{n}\right)$, and
$q_{\lambda}\left(\mathbf{U}^{s}\right)$ represent the variational
distribution, facilitating the decomposition of the model’s complexity into
manageable components. These distributions are parameterized as follows:

$\displaystyle q_{\psi}\left(\mathbf{X}\right)$
$\displaystyle=\textstyle\prod^{N}_{i=1}\mathcal{N}\left(\mathbf{x}_{i}\mid\bm{\mu}_{i},\mathbf{s}_{i}I_{Q}\right)$
(9) $\displaystyle q_{\lambda}\left(\mathbf{U}^{s}\right)$
$\displaystyle=\textstyle\prod^{K}_{k=1}\mathcal{N}\left(\mathbf{u}^{s}_{k}\mid\mathbf{m}^{s}_{k},\mathbf{S}^{s}_{k}\right)$
(10) $\displaystyle q_{\gamma}\left(\mathbf{U}^{n}\right)$
$\displaystyle=\textstyle\prod^{K}_{k=1}\mathcal{N}\left(\mathbf{u}^{n}_{k}\mid\mathbf{m}^{n}_{k},\mathbf{S}^{n}_{k}\right),$
(11)

where $\psi\triangleq\left\\{\bm{\mu}_{i},\mathbf{s}_{i}\right\\}^{N}_{i=1}$,
$\lambda\triangleq\left\\{\mathbf{m}^{s}_{k},\mathbf{S}^{s}_{k}\right\\}^{K}_{k=1}$,
and
$\gamma\triangleq\left\\{\mathbf{m}^{n}_{k},\mathbf{S}^{n}_{k}\right\\}^{K}_{k=1}$
are variational parameters that need to be optimized. Utilizing these
variational distributions and parameters, we can derive a lower bound for the
likelihood function, known as the Evidence Lower Bound (ELBO), defined by:

$\displaystyle\mathrm{ELBO}=$
$\displaystyle\mathbb{E}_{q_{\psi}\left(\mathbf{X}\right)}\left[\mathbb{E}_{p\left(\mathbf{F}^{s}\mid\mathbf{U}^{s},\mathbf{X}\right)q_{\lambda}\left(\mathbf{U}^{s}\right)}\left[\log
p\left(\mathbf{Y}^{s}\mid\mathbf{F}^{s}\right)\right]\right]$ $\displaystyle+$
$\displaystyle\mathbb{E}_{q_{\psi}\left(\mathbf{X}\right)}\left[\mathbb{E}_{p\left(\mathbf{F}^{s}\mid\mathbf{U}^{s},\mathbf{X}\right)q_{\lambda}\left(\mathbf{U}^{s}\right)}\left[\log
p\left(\mathbf{Y}^{s}\mid\mathbf{F}^{s}\right)\right]\right]$ $\displaystyle-$
$\displaystyle\mathrm{KL}\left(q_{\lambda}\left(\mathbf{U}^{s}\right)\|p\left(\mathbf{U}^{s}\right)\right)-\mathrm{KL}\left(q_{\lambda}\left(\mathbf{U}^{n}\right)\|p\left(\mathbf{U}^{n}\right)\right)$
$\displaystyle-$
$\displaystyle\mathrm{KL}\left(q_{\psi}\left(\mathbf{X}\right)\|p\left(\mathbf{X}\right)\right)$
$\displaystyle\triangleq$
$\displaystyle\mathrm{ELL}^{s}+\mathrm{ELL}^{n}-\mathrm{KL}_{u^{s}}-\mathrm{KL}_{u^{n}}-\mathrm{KL}_{X}.$
(12)

The Expected Log Likelihood (ELL) terms in the lower bound, defined in
equation (12), can be calculated using sampling method. The details of these
calculations are presented in [3]. Using a re-parameterization technique, we
use gradient-based approaches such as Adam to optimize the model’s free
parameters along with variational parameters. In summary, we build a recursive
solution to estimate variational parameters ($\psi$, $\lambda$, and $\gamma$)
and other model’s free parameters.

###  II-C Inference Process

In the inference step, our objective is to represent one or several test
points in the latent space. These test points may be neural data $Y^{n}_{*}$,
stimulus labels $Y^{s}_{*}$, or both. For new points, we define a latent
representation $X_{*}$ and initialize its variational parameters randomly.
Given a new test point, the posterior distribution
$p\left(\mathbf{X}_{*}\mid\mathbf{Y}^{n}_{*},\mathbf{Y}^{s}_{*},\mathbf{Y}^{n},\mathbf{Y}^{s}\right)$
is approximated by
$q\left(\mathbf{X}_{*}\right)=\mathcal{N}\left(\mathbf{X}_{*}\mid\bm{\mu}_{*},\mathbf{s}_{*}\mathbf{I}_{Q}\right)$.
The optimization process minimizes the negative lower bound with respect to
the variational parameters $\bm{\mu}_{*}$ and $\mathbf{s}_{*}$, while the
inducing points and other model variables are kept constant. This optimization
is expressed by:

$\operatorname*{arg\,max}_{\bm{\mu}_{*},\mathbf{s}_{*}}\left\\{\sum_{d}\mathrm{ELL}^{n}_{*,d}+\sum_{k}\mathrm{ELL}^{s}_{*,k}-\mathrm{KL}\left(q_{\psi}\left(\mathbf{X}_{*}\right)\|p\left(\mathbf{X}_{*}\right)\right)\right\\}.$
(13)

By solving this optimization problem, we obtain the estimate of the posterior
over the test points’ latent variables, i.e. mean and variance.

###  II-D Decoding Process

In the decoding step, we possess partial information about the test points,
the neural data $\mathbf{Y}^{n}_{*}$, and our goal is to determine the
probability distribution over the stimulus $\mathbf{Y}^{s}_{*}$. Similar to
the inference step, we define a corresponding latent variable for the test
point. However, this time we optimize the following expression:

$\operatorname*{arg\,max}_{\bm{\mu}_{*},\mathbf{s}_{*}}\left\\{\sum_{d}\mathrm{ELL}^{n}_{*,d}-\mathrm{KL}\left(q_{\psi}\left(\mathbf{X}_{*}\right)\|p\left(\mathbf{X}_{*}\right)\right)\right\\}.$
(14)

given we lack the label data. By optimizing the function, we determine the
optimal values for the variational parameters of $X_{*}$. Through the
generative path for the label, we can find
$p\left(\mathbf{Y}^{c}_{*}|\mathbf{X}_{*}\right)$ and effectively decode
labels.

##  III Neural Dataset

In this study, we examine the verbal memory (VerbMem) dataset recorded by
Scott Sponheim et al. [5]. The verbal memory task involves three key stages,
known as encoding, recall, and lexical decision blocks. During encoding,
participants judge the size of objects shown as words. In the recall stage,
they try to remember and recount as many words as possible. Finally, in the
lexical decision part, they distinguish between previously seen words and new
words, while their EEG data is recorded to analyze verbal memory and
recognition processes. The EEG data was collected using a 128-channel device,
and participants’ responses included both correct and incorrect choices, along
with their response times.

Our focus here is on the lexical decision blocks within the experiment. The
block consists of both new and old (seen in the encoding phase) words.
Participants must distinguish between these two categories. The goal is to
determine whether a participant is responding to new or old words based on
neural activity recordings. On average, participants’ decisions are incorrect
$17.6\pm 9.6\%$ of the time. Incorrect responses are excluded as we are
interested in instances where participants are confident in their choices. We
construct a decoder model for each individual patient, who may have a
diagnosis of schizophrenia (SZ) or are non-psychiatric control participants.
In our decoding, we remain agnostic to this information. The neural features
include time-domain features like the mean amplitudes of the N200, P300, and
post-P300 components, coherence features such as total coherency across
frequencies, coherence vectors in different brain regions, and frequency-
domain features like the average power in delta, theta, alpha, beta, and gamma
bands, totaling 1705 features per trial [6]. We selected the 25 most
correlated features based on the point- biserial correlation coefficient and
used them as input to our decoder model. In the next section, we’ll explain
how our neural decoder model works.

##  IV Decoding Results in VerbMem Dataset

In our model—and for classification purposes—we assume the latent variable
dimension is 8. For the kernels, we use Gaussian Automatic Relevance
Determination (ARD) kernels. The $\phi$ and $\theta$ parameters include length
scale and variance, which are estimated along with the approximate posterior
distributions of $\mathbf{X}$. We aggregate three blocks of data, yielding
approximately 120 data points per participant. Following the methodology
outlined in Section III, trials with incorrect decisions are excluded,
resulting in 90 to 105 trials per participant. We employ 5-fold cross-
validation to evaluate decoding accuracy. After optimizing the model
parameters on the train set, the ARD kernel parameter, as illustrated in Fig.
2.A, reveals three dominant dimensions that significantly impact the model’s
predictions. This finding indicates that the model automatically identifies
the necessary dimensions for the latent space. Fig. 2.B illustrates the
inferred $\mathbf{X}$ for training data, influenced by both neural features
(continuous path) and labels (discrete path). The test set $\mathbf{X}_{*}$ is
derived solely from neural features (see Fig. 2.C). Then it goes through the
discrete path of the generative model to determine the probability of trials
being old or new. This discrimination yields a test set accuracy of 92%.

Figure 2: Neural Features and Inferred $\mathbf{X}$: (A) ARD coefficients, indicating the relevance of each dimension in the input data to the model. (B-C) Scatter plots of latent variables for training (B) and testing (C) phases. Blue dots represent trials with new words, and red dots indicate trials with old words. TABLE I: Classification Accuracy in VERBMEM Dataset | XGBoost | Proposed model
---|---|---
| Accuracy | F-measure | Accuracy | F-measure
1 | $0.71\pm 0.11$ | $0.72\pm 0.12$ | $\mathbf{0.85\pm 0.05}$ | $\mathbf{0.86\pm 0.06}$
2 | $0.77\pm 0.15$ | $0.77\pm 0.14$ | $\mathbf{0.82\pm 0.05}$ | $\mathbf{0.82\pm 0.06}$
3 | $\mathbf{0.74\pm 0.14}$ | $\mathbf{0.74\pm 0.15}$ | $0.72\pm 0.11$ | $0.71\pm 0.11$
4 | $0.73\pm 0.08$ | $0.72\pm 0.08$ | $\mathbf{0.75\pm 0.07}$ | $\mathbf{0.75\pm 0.07}$
5 | $0.67\pm 0.04$ | $0.68\pm 0.04$ | $\mathbf{0.68\pm 0.06}$ | $\mathbf{0.69\pm 0.06}$

The training and decoding procedures were applied to data from five
participants; the mean and variance of their classification accuracy from
5-fold cross-validation are detailed in Table I The model shows consistent
performance across participants, with an average accuracy of 76.4%. To compare
the decoder performance with state-of-the-art approaches, we examined the
XGBoost classifier [7]. As shown in Table I, the model performance is
significantly better than that of the XGBoost, and this holds for 4 out of 5
participants’ data. The decoder’s performance achieves an average F-score of
76.6% across five participants, surpassing the XGBoost’s F-score of 72.6% by
4%.

##  V Discussion

In this research, we introduced a novel non-parametric neural decoder model.
The model combines GP with latent variables, letting us simultaneously
characterize neural activity and associated stimuli – or labels – as a
function of the latent variables. The framework is well suited in analysis of
datasets with limited sample sizes such as those datasets that are recorded in
neuroscience experiments. We examined the model performance in decoding
stimulus labels in the VerbMem dataset, reaching an average accuracy of 76.6%.
The performance is not the sole benefit of the model; with the model, we can
also find the optimal dimension of latent variables. Furthermore, we can
reduce the number of training samples to reduce the computational complexity
of the decoding process without a loss of accuracy. We applied the XGBoost
classifier to the same dataset; the XGBoost average prediction accuracy is
72.6% which is about 4% below our proposed approach. The preceding performance
and other model attributes highlight its potential for high-dimensional neural
data analysis. Whilst we listed the pros of the model, the model is not
without its limitations. The model solution including inference, training, and
decoding is not as straightforward as XGBoost and other ML techniques; this
might hinder its wide usage in the neuroscience community. To address this, we
have created an open-source package with synthetic and benchmark dataset
examples, which we think will facilitate the usage of the framework in other
datasets. The source code and additional resources are publicly available at
our GitHub repository 111https://github.com/Navid-Ziaei/LDGD. In this
research, we mainly focused on decoding accuracy; we think it is crucial to
study whether inferred latent variables have ethological interpretation and
how they can be studied from a biomarker discovery point of view. A challenge
with the GP model is its computational complexity when the size of the dataset
grows; whilst this was not the case in the VerbMem dataset, it is important to
examine how the model performance and computational cost scale with the sample
size and dimension of the input. These aspects of the model present new
avenues for future research, notably in exploring the model’s effectiveness in
clinical settings and its comparison with a broader range of ML techniques.
Addressing these questions will be crucial in understanding the full potential
and limitations of this innovative approach in the field of neuroscience. The
preceding performance of our latent variable Gaussian Process (GP) model in
neural data analysis opens an intriguing pathway for future research in the
realm of neural data analysis and neural decoder models. One of the crucial
directions in our research is extending the model to handle time series data.
This extension could significantly enhance the model’s applicability, as we
can detour feature extraction and pre-processing steps in the analysis of
neural data and decoding step. In parallel, we are interested in addressing
modeling challenges of latent GP models such as those that are listed in the
discussion section.

## References

  * [1] C. K. Williams and C. E. Rasmussen, Gaussian processes for machine learning, vol. 2.  MIT press Cambridge, MA, 2006.
  * [2] M. Titsias and N. D. Lawrence, “Bayesian gaussian process latent variable model,” in Proceedings of the thirteenth international conference on artificial intelligence and statistics, pp. 844–851, JMLR Workshop and Conference Proceedings, 2010.
  * [3] N. Ziaei, B. Nazari, U. T. Eden, A. Widge, and A. Yousefi, “A bayesian gaussian process-based latent discriminative generative decoder (ldgd) model for high-dimensional data,” 2024.
  * [4] H. Salimbeni and M. Deisenroth, “Doubly stochastic variational inference for deep gaussian processes,” 2017.
  * [5] C. A. Marquardt, V. J. Pokorny, S. G. Disner, N. W. Nelson, K. A. McGuire, and S. R. Sponheim, “Inefficient attentional control explains verbal-memory deficits among military veterans with posttraumatic reexperiencing symptoms,” Clinical Psychological Science, vol. 10, no. 3, pp. 499–513, 2022.
  * [6] L. Hu and Z. Zhang, “Eeg signal processing and feature extraction,” 2019.
  * [7] T. Chen and C. Guestrin, “Xgboost: A scalable tree boosting system,” in Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, pp. 785–794, 2016.


# MOD-Finder: Identify multi-omics data sets related to defined chemical
exposure

Sebastian Canzler,1 Jörg Hackermüller,1,2 Jana Schor1,∗

1Young Investigators Group Bioinformatics and Transcriptomics,
Department of Molecular Systems Biology, Helmholtz-Centre for Environmental
Research - UFZ,
Permoser Straße 15, Leipzig, 04317, Germany
2Bioinformatics Department, Universität Leipzig, Härtelstraße 16-18,
04107 Leipzig, Germany.

∗To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS>
###### Abstract

Summary: Integration of multi-omics data on chemical exposure of cells or
organisms promises a more complete representation of the responding pathways
than single omics data. Data of different omics layers, like transcriptome or
proteome is deposited in different repositories. Additionally, precisely
specifying a chemical of interest that was used in the exposure experiments
suffers from different nomenclatures and non-uniquely mapping of chemical
identifiers. The manual search for corresponding omics data sets of different
layers for exposure with a chemical of interest is thus a tedious task. We
have developed MOD-Finder (Multi-Omics Data set Finder) to efficiently search
for chemical-related omics data sets in several publicly available databases
in an automated manner. A plain and simple presentation of the returned omics
data sets is augmented with effect information that are assumed to be
triggered by the chemical of interest.
Availability and Implementation: MOD-Finder is implemented in R using the
Shiny package. The web service is available at
https://webapp.ufz.de/mod\\_finder and the source code under the GNU GPL v3
license at https://github.com/yigbt/MOD-Finder.
Supplementary information: Supplementary data are available at
https://www.ufz.de/index.php?en=44919

##  1 Introduction

Exposure of cells or organisms to hazardous chemicals triggers a series of
molecular events that may eventually lead to an adverse outcome. Understanding
the mechanisms that propagate the initial effect at the molecular level to
subsequent molecular and cellular effects and finally to the adverse outcome,
i.e. defining adverse outcome pathways, is a major field of research in
toxicology. Omics techniques have been instrumental for these analyses and as
cellular pathways combine biomolecules of different types, integration of
multi-omics data is increasingly considered in toxicology, see e.g. Buesen et
al. (2017).

Omics data are deposited in repositories that are specific for one omics
layer, e.g. transcriptomics, or proteomics and several repositories may be
used by the community for one layer. Assembling corresponding multi-omics data
thus is a tedious task that highly depends on the quality of meta data. The
Omics Discovery Index (OmicsDI) framework facilitates the aggregation of
multiple omics data sets across a variety of sources (Perez-Riverol et al.,
2017).

Identifying corresponding multi-omics data sets for exposure with a defined
chemical of interest (CoI) remains a challenge also when using OmicsDI:
Chemicals exhibit a considerable degree of naming variability, which may
include different, in some cases non-unique trivial and trademark names,
acronyms, or identifiers following standard nomenclature guidelines as
provided by the IUPAC (Krallinger et al., 2015). The reliable _digital_
identification of a CoI is therefore a prerequisite for assembling multi-omics
data. Annotated effect data of a chemical can assist in selecting the right
digital identifier for a CoI in case of naming or identifier mapping
collisions. To the best of our knowledge, a tool for assembling multi-omics
data for exposure with a CoI is currently not available.

Here, we present MOD-Finder (Multi-Omics Data set Finder). An easy-to-use web
application, to search for multi-omics data sets related to a user defined
CoI. It supports the reliable identification of the correct chemical, provides
integrated meta-information and functional annotation of the queried chemical,
and lists detected multi-omics data sets from various sources for
transcriptomic, proteomic, and metabolomic data. MOD-Finder visualizes
described molecular effects on gene and pathway level as well as related
diseases.

MOD-Finder provides known functional annotation of chemicals of interest
together with multi-omics data sets intended to be used in multi-omics data
analysis approaches in toxicology. Examples are benchmarking multi-omics data
integration software, evaluation of the benefit relative to single-omics data
analysis, hypotheses generation, and elucidation of potential mode of actions.

##  2 MOD-Finder \- Methods and workflow

The workflow is split into three parts that are visualized in Figure 1.

###  _Digital_ identification of the chemicals of interest.

CoIs are referenced using several independent digital identifiers, which often
cannot be mapped to one another straightforwardly. Many-to-many and many-to-
none relations are common. The first step in MOD-Finder considers the correct
identification of the CoI provided by the user. We implemented a search
strategy that queries all provided identifiers, or text-based entries in the
CompTox Chemistry Dashboard (Williams et al., 2017), a source of high-quality,
structure-curated, open chemical data, for a number of matching CoIs. For CoIs
not found in CompTox, MOD-Finder applies a similar search strategy in the
Pubchem111https://pubchem.ncbi.nlm.nih.gov/.

Figure 1: The MOD-Finder workflow: top) _digital_ identification and selection
of a CoI, middle) collection, integration, and visualization of chemical
related information and effect annotation, bottom) search of related omics
data sets.

### Functional annotation of chemicals of interest.

Prior knowledge on biological effects of the CoI may facilitate selection of
relevant data sets. MOD-Finder therefore queries the Comparative
Toxicogenomics Database (CTDbase) (Davis et al., 2019) for known relations
between the CoI and genes, pathways, and diseases and visualizes this data.

### Query omics data bases.

MOD-Finder queries NCBI GEO (Barrett et al., 2013) and ArrayExpress
(Kolesnikov et al., 2015), for gene expression data sets (transcriptome). The
PRIDE database (Vizcaíno et al., 2016) is queried as one of the most prominent
data repositories of mass spectrometry based proteomics data and the
MetabolomeXchange 222http://www.metabolomexchange.org international data
aggregation and notification service for metabolomics. Different data base
layouts, degrees, and types of access to data, e.g. RESTful APIs or
downloadable SQL data base files, are handled appropriately by MOD-Finder.

### Availability and maintenance.

MOD-Finder is provided as an R Shiny application hosted at
https://webapp.ufz.de/mod_finder. The source code is available at
https://github.com/yigbt/MOD-Finder. MOD-Finder is updated periodically and
upgraded to include further public omics databases or additional omics layers.
MOD-Finder releases are tracked and separately available at github.

###

MOD-Finder presents a comprehensive, clearly structured overview of the
identified data sets using visuals and customizable tables. Links to
subsequently download the data sets are provided within the app in separate
output panels.

##  3 Conclusion and discussion

Unraveling the pathway response to chemical exposure is a major aim of current
toxicological research. Integration of multi-omics data promises a much more
complete representation of pathways and therefore supports not only a deeper
mechanistic understanding of toxicological processes, but also facilitates
predictive approaches to chemical risk assessment (Buesen et al., 2017). MOD-
Finder simplifies the entry into such data-intensive approaches by enhancing
the systematic omics data set retrieval and providing prior knowledge of known
effects of the chemicals in a clear visualization. This combination is unique
and makes MOD-Finder a useful tool for multi-omics data approaches in
toxicology. To the best of our knowledge, the only other publicly available
multi-omics data retrieval tool is OmicsDI. It integrates information on
different biological entities, including genes, proteins, metabolites, and the
corresponding publications from PubMed but lacks integration of chemical
knowledge. MOD-Finder has a clear focus on retrieving multi-omics data for
chemical exposure experiments. This includes the correct computational
identification of a chemical of interest, knowledge retrieval of described
effects combined with the selection of related omics data sets from
transcriptomics, proteomics, and metabolomics measurements.

## Funding

This work was funded in part by the Cefic Long-Range Research Initiative
Program (Project C5-XomeTox).

## References

  * Barrett et al. (2013) Barrett, T. et al. (2013).  NCBI GEO: archive for functional genomics data sets–update.  Nucleic Acids Res, 41(Database issue), D991–5.
  * Buesen et al. (2017) Buesen, R. et al. (2017).  Applying ’omics technologies in chemicals risk assessment: Report of an ecetoc workshop.  Regul Toxicol Pharmacol, 91 Suppl 1, S3–S13.
  * Davis et al. (2019) Davis, A.P. et al. (2019).  The comparative toxicogenomics database: update 2019.  Nucleic Acids Res, 47(D1), D948–D954.
  * Kolesnikov et al. (2015) Kolesnikov, N. et al. (2015).  ArrayExpress update–simplifying data submissions.  Nucleic Acids Res, 43(Database issue), D1113–6.
  * Krallinger et al. (2015) Krallinger, M. et al. (2015).  Chemdner: The drugs and chemical names extraction challenge.  Journal of cheminformatics, 7, S1.
  * Perez-Riverol et al. (2017) Perez-Riverol, Y. et al. (2017).  Discovering and linking public omics data sets using the omics discovery index.  Nat Biotechnol, 35(5), 406–409.
  * Vizcaíno et al. (2016) Vizcaíno, J.A. et al. (2016).  2016 update of the pride database and its related tools.  Nucleic Acids Res, 44(D1), D447–56.
  * Williams et al. (2017) Williams, A.J. et al. (2017).  The CompTox Chemistry Dashboard: a community data resource for environmental chemistry.  J Cheminform, 9(1), 61.