# CustOmics: A versatile deep-learning based strategy for multi-omics
integration

Hakim Benkirane
Laboratory of Mathematics and Computer Science
CentraleSupélec, Université Paris-Saclay
GIf-sur-Yvettes, 91170
<EMAIL_ADDRESS>
Yoann Pradat
Laboratory of Mathematics and Computer Science
CentraleSupélec, Université Paris-Saclay
Gif-sur-Yvettes, 91170
<EMAIL_ADDRESS>
Stefan Michiels
CESP-Oncostat Team
GustaveRoussy
Villejuif, 91170
<EMAIL_ADDRESS>
Paul-Henry Cournède
Laboratory of Mathematics and Computer Science
CentraleSupélec, Université Paris-Saclay
Gif-sur-Yvettes, 91170
<EMAIL_ADDRESS>
Corresponding authorCorresponding author

###### Abstract

Recent advances in high-throughput sequencing technologies have enabled the
extraction of multiple features that depict patient samples at diverse and
complementary molecular levels. The generation of such data has led to new
challenges in computational biology regarding the integration of high-
dimensional and heterogeneous datasets that capture the interrelationships
between multiple genes and their functions. Thanks to their versatility and
ability to learn synthetic latent representations of complex data, deep
learning methods offer promising perspectives for integrating multi-omics
data. These methods have led to the conception of many original architectures
that are primarily based on autoencoder models. However, due to the difficulty
of the task, the integration strategy is fundamental to take full advantage of
the sources’ particularities without losing the global trends. This paper
presents a novel strategy to build a customizable autoencoder model that
adapts to the dataset used in the case of high-dimensional multi-source
integration. We will assess the impact of integration strategies on the latent
representation and combine the best strategies to propose a new method,
CustOmics (https://github.com/HakimBenkirane/CustOmics). We focus here on the
integration of data from multiple omics sources and demonstrate the
performance of the proposed method on test cases for several tasks such as
classification and survival analysis.

_Keywords_ Multi-omics $\cdot$ Variational Autoencoders $\cdot$ Representation
Learning

##  1 Introduction

With the advent of high-throughput technologies, multiple omics data are more
and more available to characterize the molecular portraits of patients,
notably in oncology. Analyzing this mix of data sources and leveraging their
information to improve our understanding of biological phenomena remain
challenging.

Several studies and projects have made available cohorts of data that are
characterized by several molecular sources. The Cancer Genome Atlas (TCGA,
https://portal.gdc.cancer.gov/) has profiled thousands of tumor samples for
multiple molecular assays and made available several types of data such as
genome sequencing, RNA sequencing, DNA methylation, proteomics, etc.
Integrating such diversified data is crucial to reduce the uncertainty over a
biological process due to various experimental setups and to unravel
interactions that a single source cannot extract.

This high diversity of data comes with two significant challenges. The first
one is linked to the high-dimensionality of the data. Due to the genetic
complexity of the human molecular profile, omics data generally suffers from
the ’curse of dimensionality’ resulting from the high number of features
compared to the smaller number of samples Misra et al. (2018). This high-
dimensional space often contains correlated features that result in high
redundancy, reducing the prediction performance of algorithms James et al.
(2013). The second issue lies in data heterogeneity: having its origins from
different sources and expressing different phenomena in humans’ biological
system, omics data is very diverse Bersanelli et al. (2016). For example,
transcriptomics and proteomics are normalized differently from other omics
data and use different scaling before their analysis, leading to different
ranges and data distributions. Finally, omics data like metabolomics can also
generate sparsity, as some variables can be below the detection limit and
hence be assigned null values.

To overcome those challenges and alleviate overfitting in prediction tasks,
standard methods either manually select a small subset of molecular features
based on domain knowledge Bush et al. (2009) or use dimensionality reduction
techniques before downstream analysis. However, due to the high heterogeneity
of the sources, those methods can overlook some genome-wide hidden patterns
that justify the need for more flexible methods.

In the past few years, multi-omics integration has been a very active research
subject in health science and precision medicine (Zitnik et al., 2019).
Multiple statistical learning methods have been proposed from the desire to
investigate various complex molecular systems behind cancer. One of the most
well-known methods engraved in the statistical field is the Principal
Component Analysis (PCA) (extensively presented, for example, in Jolliffe
(2002)). From this standard model, several variations have been explored,
which include Multiple Factor Analysis (Sánchez et al., ), Consensus PCA, and
multi-block PCA (Hassani et al., 2013). All those methods try to extend the
Principal Component approach to a multi-source framework by considering each
source as a distinct block. Another variation of the PCA is the Non-negative
Matrix Factorization (NMF) first introduced in Zhang et al. (2011) which uses
the same philosophy as PCA but considers a non-negative constraint instead of
an orthogonality one. Additionally, we can consider the joint Dimensionality
Reduction methods (jDR) as an extension of factorial methods studied in a
comprehensive benchmark in Cantini et al. (2021). Such methods can be divided
into several categories including Bayesian methods (Lock and Dunson, 2013)
that use assumptions on data distribution and dependencies to build a
statistical model, as well as network-based methods that rely on a network
representation and can identify modules of the disease-associated mechanisms.
The genes are represented by nodes, and edges represent the association
between those genes (Vaske et al., 2010; Bonnet et al., 2014). Other more
specific methods have been explored, mainly for clustering. For example Newman
and Cooper (2010) that presented AutoSOME which is a framework using the Self-
organizing Map approach and couples it with a density equalizer; this
equalizer reduces dimension and preserves the local topology of the gene
expressions. Another well-known multi-omics clustering method is iCluster
(Shen et al., 2012), which uses the same idea as NMF but without the non-
negative constraint and allows more diversity in data distributions.

With the rise of deep learning as the new state-of-the-art approach in many
medical applications, several studies have explored the interest of these
methods for multi-omics data integration. The primary use of these methods has
been with autoencoders Rumelhart et al. (1986), notably by Chaudhary et al.
Chaudhary et al. (2017) in the context of survival prediction. The framework
OmiVAE (Zhang et al., 2019) introduced new improvements for standard
representation learning and OmiEmbed (Zhang et al., 2021) for multitask
learning. (Simidjievski et al., 2019) also proposed multiple deep learning
architectures based on variational frameworks for multiple tasks. Huang et al.
(Huang et al., 2019) developed a neural network architecture to integrate
multi-omics data for survival by feeding the eigengene matrix instead of the
raw data. Based on these methods, several application studies on specific
cancer cases were presented, like Hira et al. (2021) that used the OmiVAE
framework to conduct an in-depth study of the Ovarian Cancer’s mechanisms.
Tong et al. (2020) also used standard autoencoders to explore multi-view
learning in survival analysis for Breast Cancer. To our knowledge, no
benchmarking study has explored and compared the different deep learning
approaches and strategies for multi-omics data integration in the context of
multitask learning.

In this work, we discuss different strategies for the integration of high-
dimensional multi-source data to learn low-dimensional latent representation
from multi-omics datasets. Moreover, we introduce a new architecture called
CustOmics that combines the advantages of the different strategies and
alleviates some of their limitations. We evaluate the impact of this new
method on different test cases: we first apply it to a pan-cancer dataset and
then study how it handles smaller datasets by evaluating it on specific
cohorts. We also provide a new package to help bioinormaticians and
computational biologists build deep learning architectures that allow an easy
switch between the multiple strategies to adapt to specific use cases.

## Materials and methods

### Representation Learning for Multi-Omics Integration

Representation Learning is a field of statistical learning that aims to
automatically discover relevant representations of the input data (Bengio et
al., 2013). In the field of multiple omics data integration, it can help
synthesize the heterogeneous distributions into a common space, revealing the
interaction between the different sources. From a biological point of view,
each omics source represents a specific view of the processes involved.
Therefore, the main goal behind integration is to uncover a broader view of
the phenomenon of interest and to better understand the underlying processes.

Each patient is characterized by $K$ omics vectors $(\mathbf{x}_{k})_{1\leq
kK}\in\mathbb{R}^{M}_{k}$ with $M_{k}$ the number of features of the $k$-th
source. The models presented in this paper aim at mapping the set of omics
vectors to a vector $\mathbf{z}\in\mathbb{R}^{m}$ with
$m\ll\sum_{k=1}^{K}M_{k}$, the latent representation. The latent features are
the components of the vector $\mathbf{z}$. We will use factorial methods as a
basis of comparison as they have proved to be the reference method for mapping
input data into a latent space of lower dimension, notably for multi-omics
integration, as illustrated by the recent review of Cantini et al. Cantini et
al. (2021).

For factorial analysis, we consider approaches based on the study of principal
components using linear operations through Multiple Factor Analysis (MFA).
Deep learning methods use autoencoder architectures to build the latent
representation by jointly training encoder and decoder functions. The types of
architecture can be divided into three main categories:

  * •

Early Integration (EI): Refers to the set of methods that aim at merging the
different sources before the dimensionality reduction. Instead of giving as
model inputs a set of separate vectors for each source, we first concatenate
the vectors of all the sources. It is characterized by its simplicity, but
finds its limitations when some sources bear more significant signals than the
others, which makes it difficult for the model to learn interactions between
sources.

  * •

Joint Integration (JI): In this setting, sub-representations are created
inside the same model for each source before learning the output. This is the
most promising strategy theoretically and is the most popular in the
literature. However, it is challenged by the sources’ heterogeneity, as they
do not necessarily follow the same learning dynamics and may need different
approaches with specific losses.

  * •

Late Integration (LI): Consists in learning an output for each source
separately, with its own representation model. It can adapt well to the
specificities of each source, but does not retrieve any cross-modality
interaction.

Figure 1: a. Early Integration: Sources are concatenated before being fed to a
single model. b. Joint Integration: Sub-representations of each source are
learned jointly before inputting into the model. c. Late Integration: Each of
the sources outputs its prediction using its independent model, the
predictions are then aggregated. d. Mixed Integration: Representation of the
mixed integration approach. In phase 1, a specific model is trained for each
source independently and embeds a sub-representation adapted to the source’s
specificities. In phase 2, the specific models are trained jointly to create
the final output.

To alleviate the issues brought by those different integration methods, we
propose a hybrid strategy, named mixed-integration, that comes halfway between
the joint and late integration. It consists of two learning phases: the first
phase independently trains the network of each source with an adapted loss to
create sub-representations. Those specific models will then be jointly trained
in a second phase to build a global output representation. That way, we can
adapt our training to each source that may have different learning dynamics in
the first phase, while still learing cross-modality interactions in the second
phase.

### Variational Autoencoders

A Variational Autoencoder (VAE) is a deep generative model which can learn
meaningful data representations from high-dimensional input data. It is an
extension of standard autoencoders in which the encoder encodes the input as a
distribution over a latent space instead of a single point.

In this case, the encoding function $q$ represents a variational distribution
(known as encoding distribution) $q_{\phi}(\mathbf{z}|\mathbf{x})$ in which
$\phi$ is the parameter to estimate and the decoding function $p$ represents
the posterior $p_{\theta}(\mathbf{x}|\mathbf{z})$.

The particularity of a VAE is its ability to encode a distribution. After the
encoding phase, there is a sampling phase in which we sample points from the
distribution $q_{\phi}(\mathbf{z}|\mathbf{x})$.

Traditionally, the distributions in the VAE architecture are supposed
Gaussian: the encoder function will learn the two parameter vectors
$\mathbf{\mu}$ and $\mathbf{\sigma}$ (respectively the means and standard
deviations) that are used to generate samples in
$q_{\phi}(\mathbf{z}|\mathbf{x})$ using the reparametrization trick
$\mathbf{z}=\mathbf{\mu}+\mathbf{\sigma}\odot\mathbf{\epsilon}$ where
$\mathbf{\epsilon}\sim\mathcal{N}(0,I)$ and $\odot$ denotes the element-wise
product.

The loss function for this architecture can be written as the sum of two
distinct losses. First, a reconstruction loss that focuses on the
autoencoder’s ability to reconstruct the data:

$\mathcal{L}_{recon}=\mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x})}[p_{\theta}(\mathbf{x}|\mathbf{z})]$
(1)

It can be interpreted as the conditional entropy of $x$ over $z$ which
quantifies the amount of uncertainty that one has over the joint distribution
$(x,z)$, knowing $z$. In a more practical way, it is related to the quantity
of information of $x$ that is retained by $z$ (thus qualifying a
reconstruction potential).

Second, a regularization loss aims at getting the encoding distribution as
close as possible to the theoretical distribution of the latent vector.
Traditionally, a Kullback-Leibler divergence is used:
$\displaystyle\mathcal{L}_{reg}=D_{KL}(q_{\phi}(\mathbf{z}|\mathbf{x})||p_{\theta}(\mathbf{z}))$.

The total loss will benefit from the framework introduced in Higgins et al.
(2017), $\displaystyle\mathcal{L}=\mathcal{L}_{recon}+\beta\mathcal{L}_{reg}$,
adding a weight on the regularization term to balance between the two parts of
the loss function, with a hyperparameter $\beta$ to optimize.

Figure 2: a. Early Integration VAE: Variational Autoencoder architecture with
early integration strategy. b. Joint Integration VAE: Variational Autoencoder
architecture with joint integration strategy. c. Late Integration VAE:
Variational Autoencoder architecture with late integration strategy.d. Mixed-
Integration/CustOmics: This is a hierarchical architecture composed of
specific per-source autoencoders that converges into a central variational
autoencoder.

### Early Integration

An early integration autoencoder is an autoencoder applied to the concatenated
data from various sources. This type of architecture was reported in Viaud et
al. (2021) as the SDAE (Standard Deep Autoencoder). The drawback of such an
integration is that the resulting concatenation of every dataset results in a
noisier and more complex high-dimensional matrix that makes the training more
difficult. Moreover, this type of integration does not consider the specific
data distribution of each source separately and thus leads the model to learn
irrelevant patterns that can, in the worst case, mask the essential signals.
For those reasons, even though this strategy has been introduced in the
benchmark, it is not recommended in a real setting for complex multi-source
datasets such as omics data. For comparison purposes, we adapt the SDAE
architecture to a Variational Autoencoder instead of a standard autoencoder as
described in Figure 2

### Joint Integration

This is a standard way to approach multi-source integration using Variational
Autoencoders. This model has been used in multiple designs in literature
(Zhang et al. (2019); Zhang et al. (2021); Simidjievski et al. (2019)) and is
probably the most popular as it is relatively straightforward. The idea is to
reduce through a series of the hidden layer the dimension of each source
separately, then concatenate them in a central layer before sampling. In this
architecture, we would learn the mean and standard deviation of the
concatenated data such that:

  * •

$\mathbf{z}=\mathcal{N}(\mu,\sigma^{2}I)$

  * •

$\mu=q_{\mu}(T(q^{(1)}(\mathbf{x}_{1};\theta_{e,1}),...,q^{(K)}(\mathbf{x}^{(K)};\theta_{e,K})))$

  * •

$\sigma=q_{\sigma}(T(q^{(1)}(\mathbf{x}^{(1)};\theta_{e,1}),...,q^{(K)}(\mathbf{x}^{(K)};\theta_{e,K})))$

### Late Integration

A late integration approach for variational autoencoders is a variation of the
architecture proposed in Simidjievski et al. (2019) with the Hierarchical VAE
architecture. The idea is to use variational inference on each layer
separately and input the resulting representations in an aggregation network
learned independently of the other VAEs, which is the main difference from the
joint integration approach. This architecture is described in Figure 2.

### Mixed Integration / CustOmics

This integration strategy will be the foundation for building the CustOmics
architecture. The proposed method is a hierarchical mixed-integration that
consists of an autoencoder for each source that creates a sub-representation
that will then be fed to a central variational autoencoder. This new
integration strategy benefits from two training phases. The first phase will
act as a normalization process: each source will train separately to learn a
more compact representation that synthesizes its information with less noise.
This will help the integration as we will lose all imbalance issues between
the sources and avoid loss of focus when a source has an inferior
dimensionality or weaker signal than the others. The second phase will
constitute a simple joint integration between the learned sub-representations,
while still training all the encoders to fine-tune those representations as
some signals are enhanced in the presence of other sources. Regarding the
regularization loss for the central layer, the KL divergence can be an
obstacle to generalization. As stated in Zhao et al. (2017), the KL divergence
suffers from various problems. The first is that the model can fail to learn a
meaningful representation of the input. Indeed, the KL divergence can
sometimes be too restrictive and naturally tends to make the latent code a
random sample from $p_{\theta}(Z)$. The second is that the KL divergence can
make the model overfit and learn a latent code with variance tending to
infinity.

So instead, we will use the Maximum Mean Discrepancy (MMD) to assess the
distance between the distributions. This distance stands on the foundation
that two distributions are identical if and only if their moments are the
same. Let $p,q$ be two distributions, the MMD distance is given as follows:

$MMD(p(x)||q(x))=\mathbb{E}_{p(x),p(x^{\prime})}(\kappa(x,x^{\prime}))+\mathbb{E}_{q(x),q(x^{\prime})}(\kappa(x,x^{\prime}))-2\mathbb{E}_{p(x),q(x^{\prime})}(\kappa(x,x^{\prime}))$
(2)

where $\kappa$ is a kernel function. We will choose a Gaussian kernel
$\kappa(x,x^{\prime})=e^{\frac{-||x-x^{\prime}||^{2}}{2\sigma^{2}}}$.

Whereas all the models mentioned previously build the latent representation in
an unsupervised way, we also create latent features adapted to specific tasks
like classification or survival. This idea has been used multiple times in the
literature, for example in Zhang et al. (2021).

The solution relies on adding a task-related loss to the autoencoder objective
function. Therefore, we denote by $\mathcal{L}_{task}$ the loss such that
$\displaystyle\mathcal{L}_{tot}=\mathcal{L}_{AE}+\alpha\mathcal{L}_{task}$.

For the classification task, we use a categorical cross-entropy loss defined
by $\displaystyle\mathcal{L}_{class}=\sum_{i}y_{i}log(\hat{y}_{i})$, where
$y_{i}$ is the ground truth for the $i^{th}$ sample, and $\hat{y}_{i}$ its
estimation with the downstream model.

For the survival task, we use the DeepSurv loss function. This nonlinear
proportional hazard model is introduced in Katzman et al. (2018). The model is
built by using the negative partial log-likelihood formula that translates in
our case into:

$L(\theta)=-\sum_{i:E_{i}=1}(\hat{\mu}(t,x_{i};\theta)-log\sum_{j\in\mathcal{R}(T_{i})}e^{\hat{\mu}(x_{i};\theta)})$
(3)

where $\hat{\mu}(x;\theta)$ is the risk function estimated by the output layer
of the network, $\mathcal{R}(t)$ is the risk set, that is the set of patients
still at risk of failure after time $t$.

Additionally, this loss will be assigned to each omic-specific network in the
first training phase to create adequate sub-representations before the joint
integration phase.

### Interpretability

To build a more interpretable architecture, we adapt the method introduced in
Withnell et al. (2021) to compute SHAP values (Lundberg and Lee, 2017) for
deep variational autoencoders. Whereas this method has only been conceived for
single source inputs, we expand it to the multi-source setting of CustOmics by
adapting it to any deep autoencoder and applying it to each source autoencoder
thanks to the dual-phase approach characterizing the mixed-integration
strategy. After training the CustOmics network, we pass each intermediate
autoencoder to a DeepSHAP explainer that was modified to make sure that it
could take either each component of the latent vector or the prediction output
for the contribution analysis. For each prediction, SHAP values corresponding
to each gene or latent dimension are computed to estimate the overall
contribution by averaging them over a group of samples with similar features.
When combined with the mixed-integration approach through CustOmics, we can
compute for each source two types of features importance. During the first
training phase, we compute the features importances of each source separately
to assess the importance of genes when only considering this source. During
the second training phase, the features’ importance can change as it will then
consider other sources, showing how some genes can have increasing or
decreasing importance when introducing complementary views of the biological
system.

### Test Cases and Datasets

In this study, we use datasets extracted from the Genomic Data Commons (GDC)
pan-cancer multi-omics study Grossman et al. (2016). It is one of the most
comprehensive datasets for multi-omics analysis, with high-dimensional omics
data and corresponding phenotype data from The Cancer Genome Atlas (TCGA). In
our experiments, we use three types of omics data, Copy Number Variations
(CNV), RNA-Seq gene expressions, and DNA methylation. The CNV dataset
comprises Gistic2 measurements on a total of 19,729 genes. The RNA-Seq
expression dataset profile comprises around 60,484 identifiers referring to
corresponding exons and measuring log2 transformed Fragments Per Kilobase of
transcript per Million mapped reads (FPKM). Finally, the DNA methylation
dataset was produced using the Infinium HumanMethylation450 BeadChip (450K)
arrays with 485,578 probes in which beta values of probes indicate the
methylation ratio of corresponding CpG sites. Moreover, we also evaluate our
method on six smaller cohorts from TCGA: Bladder Urothelial Carcinoma (BLCA,
n=437), Breast Invasive Carcinoma (BRCA, n=1022), Lung Adenocarcinoma (LUAD,
n=498), Glioblastoma & Lower Grade Glioma (GBMLGG, n=515) and Uterine Corpus
Endometrial Carcinoma (UCEC, n=538). We will perform in this study 4 different
evaluations on 4 test cases for clasification and survival. The first task in
our set of experiments is classifying the different tumor types in the pan-
cancer study. The classification performance was measured using five metrics:
Accuracy, Macro-averaged F1 score, precision, recall, and ROC-AUC Yang et al.
(2021). We also perform a second classification task to validate our findings
on a smaller dataset and test the robustness of our method. This task aims at
performing a tumor subtype classification based on the PAM50 classification
(IlluminlA, ILluminalB, Basal, and HER2), we use the same setup as the pan-
cancer case. The third test case will be a survival study of the Pancancer
dataset. Finally, the fourth test case will be about evaluating the survival
performances on the 5 datasets presented earlier. For all those test cases, We
compare the CustOmics model to several reference methods for multi-omics
integration: first with a combination of dimensionality reduction methods,
Multiple Factor Analysis (MFA), Uniform Manifold Approximation and Projection
(UMAP) and non-negative matrix factorization (NMF), and also with different
deep learning methods corresponding to the various strategies described in the
Methods section.

### Data Preprocessing

The data required some preprocessing before analysis.

  * •

For RNA gene expression profiles, 594 exons located on the Y chromosome were
removed, along with 1,904 ones with zero expression and 248 with missing
values.

  * •

For DNA methylation data, the same strategy as with gene expression profiles
was used, in addition to removing probes that cannot be mapped to the human
reference genome. It leaves us with 438,831 CpG sites.

Afterward, we intersected each combination of omics data in order to retrieve
the maximum number of samples for each test case. We then identified and
removed features with missing/consistently zero/NA values for other omic
files. Finally, for non-normalized datasets such as CNVs and RNA-Seq data, we
applied a min-max normalization to ensure that each omic source was scaled
identically and thus would have the same importance during integration.

### Implementation Details

The CustOmics framework is based on the Pytorch deep-learning library. It can
be applied to any combination of high-dimensional datasets and incorporates
different integration strategies depending on the type of data and task to
perform. As done in Zhang et al. Zhang et al. (2021), DNA methylation data can
be divided into 23 separate blocks, each feeding a hidden layer corresponding
to a chromosome to avoid overfitting and save GPU memory.

The whole architecture is built using fully-connected blocks. We use a batch
normalization technique in each layer composing the neural network to address
the internal covariate shift problemIoffe and Szegedy (2015). Also, to avoid
overfitting problems, we use dropout Srivastava et al. (2014); its rate is
considered a hyperparameter.

The input dataset was randomly split into training, validation, and testing
set (60-20-20%) by using a stratified 5-fold cross-validation so that the
proportion of samples in each tumor type between the different sets is
preserved in all the folds. We perform Bayesian optimization Snoek et al.
(2012) using the validation set to find our model’s best possible combination
of hyperparameters.

All models were trained using an Nvidia Tesla V100S with 32 GB memory.

## Results

### Classification Results

We first perform the classification task on the pan-cancer dataset. Each
architecture is coupled with an artificial neural network classifier composed
of two hidden layers with 256 and 128 neurons.

Figure 3: a. pan-cancer and PAM50 classification results: Overall classification results for the pan-cancer tumor classification test case and the PAM50 subtype classification for breast cancer. b. T-SNE vizualization for each omic source separately, along with the latent representation constructed by CustOmics. We see that the constructed layer representation succeeds at separating the data into 4 distinct clusters that we couldn’t distinguish with each omic source alone. c. PAM50 gene importance: Computed SHAP values on the RNA-Seq data of the most relevant genes responsible for the discrimination between each subtype against the others using CustOmics for both integration phases. Table 1: The classification performance for the pan-cancer dataset is evaluated with 5 standard metrics for UMAP, NMF, MFA, and deep-learning methods. We evaluate the performances on the final predicted output of the downstream classifier. Model | Accuracy | F1-score | Precision | Recall | ROC-AUC
---|---|---|---|---|---
UMAP | $0.7598\pm 0.0036$ | $0.7149\pm 0.0029$ | $0.7200\pm 0.0031$ | $0.7598\pm 0.0032$ | $0.8740\pm 0.0012$
NMF | $0.8599\pm 0.0017$ | $0.8406\pm 0.0013$ | $0.8460\pm 0.0018$ | $0.8599\pm 0.0021$ | $0.9266\pm 0.0019$
MFA | $0.9167\pm 0.0012$ | $0.9025\pm 0.0014$ | $0.8945\pm 0.0008$ | $0.9167\pm 0.0013$ | $0.9565\pm 0.0003$
Early Int. VAE | $0.9337\pm 0.0079$ | $0.9314\pm 0.0086$ | $0.9367\pm 0.0067$ | $0.9337\pm 0.0079$ | $0.9655\pm 0.0041$
Joint Int. VAE | $0.9610\pm 0.0032$ | $0.9600\pm 0.0032$ | $0.9631\pm 0.0043$ | $0.9610\pm 0.0032$ | $0.9898\pm 0.0005$
Late Int. VAE | $0.9492\pm 0.0115$ | $0.9464\pm 0.0111$ | $0.9498\pm 0.0079$ | $0.9492\pm 0.0115$ | $0.9737\pm 0.0060$
CustOmics | $\mathbf{0.9788\pm 0.0025}$ | $\mathbf{0.9705\pm 0.0033}$ | $\mathbf{0.9728\pm 0.0041}$ | $\mathbf{0.9685\pm 0.0034}$ | $\mathbf{0.9918\pm 0.0001}$

Figure 6 and table 1 show the overall classification results (More details can
be found in Table S1). Among the factorial methods, the MFA achieved the best
results, so we coupled this method with the same ANN classifier used with the
deep-learning representation methods as a basis for comparison. However, it
does not perform as well as most deep-learning methods. It is because MFA
cannot uncover nonlinear relationships between different sources, unlike deep-
learning architectures. Moreover, as the MFA is an early integration, it
suffers from the problems mentioned above when we introduced this strategy.

As hinted earlier, we can see that for the deep learning strategies, early
integration is behind the others in terms of performance. It can be explained
by the fact that RNA-Seq data hold more signals when determining tumor types
or subtypes. Thus, concatenating the sources before feeding them to the VAE
overshadows the other sources, and the learned representation depends mostly
on RNA-Seq data without leveraging the other modalities. It is illustrated in
Table S2, showing that the classification results using RNA-Seq data only are
very close to those obtained with early integration, indicating that the model
may overlook the interactions between sources. Late integration is not optimal
either, since interactions between sources are not properly learned. Joint
integration performs well in most cases, but we see that the best results
displayed in Table S2 are achieved by the combination of only two sources,
RNA-Seq and Methylation data, as it seems that CNV data only adds noise to the
latent representation, meaning that its information is not handled well with
this strategy.

Table 2: Classification performance for PAM50 classification on the TCGA-BRCA dataset, with 5 standard metrics. We compare machine learning methods like UMAP, NMF, and MFA with deep-learning methods. We evaluate the performances on the final predicted output of the downstream classifier. Model | Accuracy | F1-score | Precision | Recall | ROC-AUC
---|---|---|---|---|---
UMAP | $0.6815\pm 0.0152$ | $0.6612\pm 0.0140$ | $0.6637\pm 0.0157$ | $0.6815\pm 0.0151$ | $0.8482\pm 0.0055$
NMF | $0.7025\pm 0.0132$ | $0.6955\pm 0.0135$ | $0.7032\pm 0.0110$ | $0.7025\pm 0.0131$ | $0.8576\pm 0.0119$
MFA | $0.7532\pm 0.0164$ | $0.7460\pm 0.0160$ | $0.7562\pm 0.0162$ | $0.7531\pm 0.0165$ | $0.8884\pm 0.0037$
Early Int. VAE | $0.8063\pm 0.0152$ | $0.7840\pm 0.0167$ | $0.8228\pm 0.0167$ | $0.8063\pm 0.0150$ | $0.9552\pm 0.0077$
Joint Int. VAE | $0.8518\pm 0.0184$ | $0.8488\pm 0.0189$ | $0.8589\pm 0.0161$ | $0.8518\pm 0.0151$ | $0.9734\pm 0.0035$
Late Int. VAE | $0.8312\pm 0.0174$ | $0.8124\pm 0.0201$ | $0.8429\pm 0.0215$ | $0.8312\pm 0.0176$ | $0.9689\pm 0.0066$
CustOmics | $\mathbf{0.8758\pm 0.0162}$ | $\mathbf{0.8728\pm 0.0141}$ | $\mathbf{0.9012\pm 0.0137}$ | $\mathbf{0.8758\pm 0.0130}$ | $\mathbf{0.9828\pm 0.0022}$

These results confirm the interest of the CustOmics architecture, as it gives
the best performances for all the test cases. Not only does it perform better
than the others, but it also takes advantage of the complementarity and
interactions between sources: As shown in Table S1, all sources bring an
additional amount of information. Figure 6 gives a visualization of the
different sources: Even though the initial sources are quite entangled, the
CustOmics latent representation succeeds in separating the clusters using the
mutual information between modalities, resulting in better classification
performances. We also use the interpretability property of CustOmics
introduced in the Methods section to highlight the most relevant features for
the discrimination between PAM50 subtypes by computing their respective SHAP
values for each source. We do it for both phases: in phase 1, we retrieve the
relevant genes considered when using a single omic source, whereas, in phase
2, we investigate how the addition of other sources’ signals changes the
genes’ importance. Figure 6 references the results of such explanations on
RNA-Seq data, Figure S1 and S2 show the results for CNV and methylation data.
We witness some well-referenced biomarkers for breast cancer like TFF1 Yi et
al. (2020), suggesting that our method can retrieve relevant biological
information.

### Survival Analysis

The second task in this study is survival analysis. The goal is to predict the
risk score associated with each patient from the corresponding high-
dimensional omics data. Two standard metrics evaluate the performance of this
downstream task: The c-index, which is a generalization of the AUC metric for
censored data (Harrell et al., 1982), and the Integrated Brier Score (Graf et
al., 1999).

Figure 4 and table 3 show the results for the different methods (more details
can be found in Table S3) for the survival task. The same observations as for
the classification task can be done regarding the differences between
integration strategies. Here again, we also evaluated the performance of the
CustOmics method for each combination of omics sources (Table S3).

Figure 4: a. Survival Analysis Performances: We evaluate the performances of the survival model for the pan-cancer dataset using both C-index and the Integrated Brier Score (IBS). Here again, our model outperforms the other integration strategies on both metrics. b. Log-rank test: We compute the p-value associated to the log-rank test between high and low risk groups for every integration strategy on a validation set for the pan-cancer survival test case, and we compare it to mono-omic survival predictions. c. Kaplan Meier Curves : We draw the Kaplan Meier curves and display the p-value associated to the log-rank test as computed previously for each dataset using the predicted hazard from the CustOmics model and stratify the population into high and low risk on the test set for the predicted hazard ratio. This figure shows that our method succeeds in stratifying the patients into risk subgroups. Table 3: The survival analysis performance for the pan-cancer dataset is evaluated with 2 standard metrics, C-index and IBS. We compare classical methods like UMAP, NMF, and MFA with deep-learning methods and evaluate the performances on the final predicted output of the downstream survival network. Model | C-index | IBS
---|---|---
UMAP | $0.5948\pm 0.0231$ | $0.2486\pm 0.0327$
NMF | $0.6012\pm 0.0204$ | $0.2207\pm 0.0264$
MFA | $0.6127\pm 0.0164$ | $0.2192\pm 0.0203$
Early Int. VAE | $0.6578\pm 0.0103$ | $0.2106\pm 0.0117$
Joint Int. VAE | $0.6709\pm 0.0041$ | $0.1802\pm 0.0072$
Late Int. VAE | $0.6629\pm 0.0086$ | $0.2112\pm 0.0088$
CustOmics | $\mathbf{0.6841\pm 0.0033}$ | $\mathbf{0.1745\pm 0.0052}$

The last task consists in evaluating the model performances for survival
analysis for several specific cancer types of the TCGA datasets described in
the dataset section. The objective is to evaluate the robustness of the models
when dealing with smaller datasets.

Table 4: Performances of state-of-the-art integration methods for survival
analysis, using concordance index on 5 TCGA cohorts: Bladder Urothelial
Carcinoma (BLCA), Breast Invasive Carcinoma (BRCA), Glioblastoma & Lower Grade
Glioma (GBMLGG), Lung Adenocarcinoma (LUAD) and Ulterine Corpus Endometrial
Carcinoma (UCEC).

Model | BLCA | BRCA | GBMLGG | LUAD | UCEC | Overall
---|---|---|---|---|---|---
UMAP | $0.527\pm 0.048$ | $0.524\pm 0.039$ | $0.557\pm 0.028$ | $0.530\pm 0.017$ | $0.543\pm 0.025$ | $0.536$
NMF | $0.553\pm 0.060$ | $0.560\pm 0.036$ | $0.584\pm 0.067$ | $0.589\pm 0.052$ | $0.570\pm 0.068$ | $0.571$
MFA | $0.591\pm 0.052$ | $0.599\pm 0.043$ | $0.605\pm 0.036$ | $0.597\pm 0.058$ | $0.586\pm 0.020$ | $0.596$
Early Int. VAE | $0.603\pm 0.054$ | $0.618\pm 0.039$ | $0.628\pm 0.056$ | $0.612\pm 0.041$ | $0.609\pm 0.032$ | $0.614$
Joint Int. VAE | $0.616\pm 0.072$ | $0.627\pm 0.030$ | $0.635\pm 0.020$ | $0.608\pm 0.038$ | $0.630\pm 0.021$ | $0.624$
Late Int. VAE | $0.610\pm 0.055$ | $0.620\pm 0.057$ | $0.621\pm 0.067$ | $0.595\pm 0.012$ | $0.627\pm 0.022$ | $0.615$
CustOmics | $\mathbf{0.637\pm 0.050}$ | $\mathbf{0.633\pm 0.018}$ | $\mathbf{0.642\pm 0.028}$ | $\mathbf{0.625\pm 0.037}$ | $\mathbf{0.667\pm 0.022}$ | $\mathbf{0.640}$

Finally, we perform a more thorough analysis on the survival results that we
display in Figure 4. We leave out 20% of our datasets for validation purpose,
and we perform 5-fold cross-validation on the remaining 80% to compute the
p-values associated to the log-rank test for different combinations of the Pan
cancer test case. We show the ability of CustOmics to stratify the patients
into distinct risk groups using the predicted hazard ratio. This
stratification ability was also measured quantitatively using the p-value
associated to the log-rank test between the different categories. Even though
the comparison between joint, late and mixed integration is not significant in
this case, it is interesting to note that the addition of multiple omics
sources has greatly affected the p-value as it was nearly multiplied by a
factor $10^{-5}$ for CNV data and $10^{-2}$ for RNA-Seq and methylation data.
Those results also show that the early integration is the strategy with
minimal enrichment from other sources, which corroborates our previous
intuition and is coherent with the results found in the different experiments.

## Discussion

In this work, we presented a range of integration strategies for multi-source
data that can handle both the high dimensionality and the heterogeneity of the
data. To take the best of all those strategies, we presented the mixed-
integration and the CustOmics framework to alleviate the limitations of the
existing methods. This new framework can achieve better latent representations
and lead to a more robust and generalizable architecture, as shown with the
systematic better results than alternative strategies. Importantly, our method
is capable of adapting to each omic source by handling the training
independently in a first phase, which solves the issue of unbalanced signals
between the sources by standardizing the representations before starting to
learn cross-modality interactions. Our fusion model demonstrated better
performances for both classification and survival outcome prediction across
all test cases. It not only achieved great results on pan-cancer data, but it
also did on smaller dataset for particular cohorts, showing the robustness of
our method in situations with fewer samples. Interestingly, the performances
of CustOmics in comparison to other integration strategies for different
combinations of omic sources show that the mixed-integration handles better
the fusion between CNV, RNA-Seq and methylation data. This aspect could
probably be further enhanced by adding prior knowledge of those sources in the
intermediate autoencoders, for example introducing a negative binomial prior
in the RNA-Seq autoencoder. Furthermore, by adapting the SHAP method to our
architecture, we were able to highlight genes of importance for specific
tasks. However, there is room for improvement from a computational point of
view, since the interpretability method is very expensive.

As clinical applications often suffer from a lack of samples to build
efficient deep learning models and as all the modalities are not always
extracted for all patients, it would be beneficial to study the benefits of
per-source transfer learning during phase 1 to pick-up weaker signals and
eliminate the noise. In fact, one of the advantages of phase 1 learning the
sub-representations independently is that it is not necessary for all patients
to have all omic sources available at all times, which helps us take the best
advantage of the data at hand.

## Conclusion

In conclusion, our generic and interpretable multi-source deep learning
framework improves on state-of-the-art integration strategies by proposing a
hybrid approach that fits well with multi-omics data. The framework is
available to use on Github: https://github.com/HakimBenkirane/CustOmics.
Possible steps to refine CustOmics are a more thorough study on the
interpretability part of the framework and adding prior knowledge to
individual autoencoders to enhance the sub-representations of omic sources
during phase 1 of mixed-integration.

## Supplementary Materials

### Classification Results

Table 5: Classification performances for multiple combinations of omics data using CustOMICS on the Pancancer dataset. We can see that the best performances are brought by Transcriptomics data, but the addition of other omics data increases the performances, showing that the integration is relevant. Omics | Accuracy | F1-score | Precision | Recall | ROC-AUC
---|---|---|---|---|---
CNV | $0.47\pm 0.03$ | $0.47\pm 0.03$ | $0.47\pm 0.03$ | $0.48\pm 0.01$ | $0.75\pm 0.02$
RNAseq | $0.92\pm 0.01$ | $0.93\pm 0.01$ | $0.93\pm 0.01$ | $0.92\pm 0.01$ | $0.96\pm 0.00$
methyl | $0.68\pm 0.02$ | $0.68\pm 0.02$ | $0.68\pm 0.02$ | $0.68\pm 0.02$ | $0.82\pm 0.01$
CNV + RNAseq | $0.93\pm 0.01$ | $0.92\pm 0.01$ | $0.92\pm 0.01$ | $0.92\pm 0.01$ | $0.96\pm 0.00$
CNV + methyl | $0.71\pm 0.02$ | $0.69\pm 0.02$ | $0.70\pm 0.02$ | $0.69\pm 0.02$ | $0.85\pm 0.01$
RNAseq + methyl | $0.94\pm 0.01$ | $0.94\pm 0.01$ | $0.94\pm 0.01$ | $0.94\pm 0.0132$ | $0.97\pm 0.00$
CNV + RNAseq + methyl | $0.98\pm 0.01$ | $0.97\pm 0.01$ | $0.97\pm 0.01$ | $0.97\pm 0.01$ | $0.99\pm 0.00$
Table 6: Classification performances for multiple combinations of omics data using joint integration on the Pancancer dataset. Omics | Accuracy | F1-score | Precision | Recall | ROC-AUC
---|---|---|---|---|---
CNV | $0.47\pm 0.03$ | $0.47\pm 0.03$ | $0.47\pm 0.03$ | $0.48\pm 0.01$ | $0.75\pm 0.02$
RNAseq | $0.92\pm 0.01$ | $0.93\pm 0.01$ | $0.93\pm 0.01$ | $0.92\pm 0.01$ | $0.96\pm 0.00$
methyl | $0.68\pm 0.02$ | $0.68\pm 0.02$ | $0.68\pm 0.02$ | $0.68\pm 0.02$ | $0.82\pm 0.01$
CNV + RNAseq | $0.90\pm 0.02$ | $0.89\pm 0.02$ | $0.90\pm 0.02$ | $0.88\pm 0.03$ | $0.93\pm 0.00$
CNV + methyl | $0.70\pm 0.02$ | $0.69\pm 0.02$ | $0.70\pm 0.02$ | $0.70\pm 0.02$ | $0.85\pm 0.01$
RNAseq + methyl | $0.96\pm 0.01$ | $0.96\pm 0.01$ | $0.96\pm 0.01$ | $0.96\pm 0.0132$ | $0.99\pm 0.00$
CNV + RNAseq + methyl | $0.94\pm 0.01$ | $0.94\pm 0.01$ | $0.95\pm 0.01$ | $0.94\pm 0.01$ | $0.97\pm 0.00$
Figure 5: PAM50 gene importance: Computed SHAP values on the CNV data of the
most relevant genes responsible for the discrimination between each subtype
against the others using CustOmics for both integration phases. Figure 6:
PAM50 gene importance: Computed SHAP values on the methylation data of the
most relevant genes responsible for the discrimination between each subtype
against the others using CustOmics for both integration phases.

### Survival Results

Table 7: Survival performances for multiple combinations of omics data using CustOMICS on the Pan cancer dataset. We can see that the best performances are brought by Transcriptomics data, but the addition of other omics data increases the performances, showing that the integration is relevant. Omics | Accuracy | F1-score
---|---|---
CNV | $0.54\pm 0.05$ | $0.25\pm 0.06$
RNASeq | $0.63\pm 0.02$ | $0.20\pm 0.02$
methyl | $0.59\pm 0.02$ | $0.23\pm 0.03$
CNV + RNAseq | $0.64\pm 0.02$ | $0.19\pm 0.02$
CNV + methyl | $0.61\pm 0.03$ | $0.21\pm 0.03$
RNAseq + methyl | $0.64\pm 0.02$ | $0.19\pm 0.02$
CNV + RNAseq + methyl | $0.68\pm 0.01$ | $0.17\pm 0.01$

## References

  * Misra et al. [2018] Biswapriya Misra, Carl Langefeld, Michael Olivier, and Laura Cox.  Integrated omics: Tools, advances, and future approaches.  _Journal of Molecular Endocrinology_ , 61(1):R21–45, 2018.
  * James et al. [2013] Gareth James, Daniela Witten, Trevor Hastie, and Robert Tibshirani.  _An Introduction to Statistical Learning: with Applications in R_.  Springer, 2013.
  * Bersanelli et al. [2016] Matteo Bersanelli, Ettore Mosca, Daniel Remondini, Enrico Giampieri, Claudia Sala, Gastone Castellani, and Luciano Milanesi.  Methods for the integration of multi-omics data: Mathematical aspects.  _BMC Bioinformatics_ , 17(2):S15, 2016.
  * Bush et al. [2009] William S Bush, Scott M Dudek, and Marylyn D Ritchie.  Biofilter: a knowledge-integration system for the multi-locus analysis of genome-wide association studies.  In _Biocomputing 2009_. World Scientific, 2009.
  * Zitnik et al. [2019] Marinka Zitnik, Francis Nguyen, Bo Wang, Jure Leskovec, Anna Goldenberg, and Michael M Hoffman.  Machine learning for integrating data in biology and medicine: Principles, practice, and opportunities.  _Information Fusion_ , 50:71–91, 2019.
  * Jolliffe [2002] I.T. Jolliffe.  _Principal Component Analysis_.  Springer, 2002.
  * [7] Alex Sánchez, José Fernández-Real, Esteban Vegas, Francesc Carmona, Jacques Amar, Remy Burcelin, Matteo Serino, Francisco Tinahones, M Carmen Ruíz de Villa, Antonio Minãrro, et al.  Multivariate methods for the integration and visualization of omics data.  In _Spanish Symposium on Bioinformatics_ , page 382. Springer.
  * Hassani et al. [2013] Sahar Hassani, Mohamed Hanafi, El Mostapha Qannari, and Achim Kohler.  Deflation strategies for multi-block principal component analysis revisited.  _Chemometrics and Intelligent Laboratory Systems_ , 120:154–68, 2013.
  * Zhang et al. [2011] Shihua Zhang, Qingjiao Li, Juan Liu, and Xianghong Zhou.  A novel computational framework for simultaneous integration of multiple types of genomic data to identify microrna-gene regulatory modules.  _Bioinformatics_ , 27(13):i401–09, 2011.
  * Cantini et al. [2021] Laura Cantini, Pooya Zakeri, Celine Hernandez, Aurelien Naldi, Denis Thieffry, Elisabeth Remy, and Anaïs Baudot.  Benchmarking joint multi-omics dimensionality reduction approaches for the study of cancer.  _Nature Communications_ , 12(1):124, 2021.
  * Lock and Dunson [2013] Eric Lock and David Dunson.  Bayesian consensus clustering.  _Bioinformatics_ , 29(20):2610–16, 2013.
  * Vaske et al. [2010] Charles Vaske, Stephen Benz, J Sanborn, Dent Earl, Christopher Szeto, Jingchun Zhu, David Haussler, and Joshua Stuart.  Inference of patient-specific pathway activities from multi-dimensional cancer genomics data using paradigm.  _Bioinformatics_ , 26(12):i237–45, 2010.
  * Bonnet et al. [2014] Eric Bonnet, Laurence Calzone, and Tom Michoel.  Integrative multi-omics module network inference with lemon-tree.  _PLoS computational biology_ , 11(2):e1003983, 2014.
  * Newman and Cooper [2010] Aaron Newman and James Cooper.  Autosome: A clustering method for identifying gene expression modules without prior knowledge of cluster number.  _BMC bioinformatics_ , 11(1):117, 2010.
  * Shen et al. [2012] Ronglai Shen, Quincy Mo, Nikolaus Schultz, Venkatraman Seshan, Adam Olshen, Jason Huse, Marc Ladanyi, and Chris Sander.  Integrative subtype discovery in glioblastoma using icluster.  _PloS one_ , 7(4):e35236, 2012.
  * Rumelhart et al. [1986] D. E. Rumelhart, G. E. Hinton, and R. J. Williams.  _Learning Internal Representations by Error Propagation_.  MIT Press, 1986.
  * Chaudhary et al. [2017] Kumardeep Chaudhary, Olivier Bertrand Poirion, Liangqun Lu, and Lana Garmire.  Deep learning-based multi-omics integration robustly predicts survival in liver cancer.  _Clinical Cancer Research_ , 24(6):1248–59, 2017\.
  * Zhang et al. [2019] Xiaoyu Zhang, Jingqing Zhang, Kai Sun, Xian Yang, Chengliang Dai, and Yike Guo.  Integrated multi-omics analysis using variational autoencoders: Application to pan-cancer classification.  In _2019 IEEE International Conference on Bioinformatics and Biomedicine (BIBM)_ , page 765–69, 2019.
  * Zhang et al. [2021] Xiaoyu Zhang, Yuting Xing, Kai Sun, and Yike Guo.  Omiembed: A unified multi-task deep learning framework for multi-omics data.  _Cancers_ , 13(12):3047, 2021.
  * Simidjievski et al. [2019] Nikola Simidjievski, Cristian Bodnar, Ifrah Tariq, Paul Scherer, Helena Andres Terre, Zohreh Shams, Mateja Jamnik, and Pietro Liò.  Variational autoencoders for cancer data integration: Design principles and computational practice.  _Frontiers in Genetics_ , 10:1205, 2019.
  * Huang et al. [2019] Zhi Huang, Xiaohui Zhan, Shunian Xiang, Travis S. Johnson, Bryan Helm, Christina Y. Yu, Jie Zhang, Paul Salama, Maher Rizkalla, Zhi Han, and Kun Huang.  Salmon: Survival analysis learning with multi-omics neural networks on breast cancer.  _Frontiers in Genetics_ , 10:166, 2019.
  * Hira et al. [2021] Muta Hira, Mohammad Abdur Razzaque, Claudio Angione, James Scrivens, Saladin Sawan, and Mosharraf Sarkar.  Integrated multi-omics analysis of ovarian cancer using variational autoencoders.  _Scientific Reports_ , 11(1):6265, 2021.
  * Tong et al. [2020] Li Tong, Jonathan Mitchel, Kevin Chatlin, and May Wang.  Deep learning based feature-level integration of multi-omics data for breast cancer patients survival analysis.  _BMC medical informatics and decision making_ , 20(1):225, 2020.
  * Bengio et al. [2013] Yoshua Bengio, Aaron Courville, and Pascal Vincent.  Representation learning: A review and new perspectives.  _IEEE transactions on pattern analysis and machine intelligence_ , 35(8):1798–828, 2013.
  * Higgins et al. [2017] Irina Higgins, Loic Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Matthew Botvinick, Shakir Mohamed, and Alexander Lerchner.  beta-VAE: Learning basic visual concepts with a constrained variational framework.  In _International Conference on Learning Representations_ , 2017.
  * Viaud et al. [2021] G. Viaud, P. Mayilvahanan, and P. Cournede.  Representation learning for the clustering of multi-omics data.  _IEEE/ACM Transactions on Computational Biology and Bioinformatics_ , 19:135–45, 2021.
  * Zhao et al. [2017] Shengjia Zhao, Jiaming Song, and Stefano Ermon.  Infovae: Information maximizing variational autoencoders.  _arXiv preprint arXiv:1706.02262_ , 2017.
  * Katzman et al. [2018] Jared Katzman, Uri Shaham, Alexander Cloninger, Jonathan Bates, Tingting Jiang, and Yuval Kluger.  Deepsurv: Personalized treatment recommender system using a cox proportional hazards deep neural network.  _BMC Medical Research Methodology_ , 18(1):24, 2018.
  * Withnell et al. [2021] Eloise Withnell, Xiaoyu Zhang, Kai Sun, and Yike Guo.  Xomivae: an interpretable deep learning model for cancer classification using high-dimensional omics data.  _Briefings in Bioinformatics_ , 22:315, 2021.
  * Lundberg and Lee [2017] Scott M Lundberg and Su-In Lee.  A unified approach to interpreting model predictions.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems 30_ , pages 4765–4774. Curran Associates, Inc., 2017.
  * Grossman et al. [2016] Robert L. Grossman, Allison P. Heath, Vincent Ferretti, Harold E. Varmus, Douglas R. Lowy, Warren A. Kibbe, and Louis M. Staudt.  Toward a shared vision for cancer genomic data.  _New England Journal of Medicine_ , 375(12):1109–12, 2016.
  * Yang et al. [2021] Zhiyong Yang, Qianqian Xu, Shilong Bao, and Xiaochun Cao.  Learning with multiclass auc: Theory and algorithms.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , PP:1–1, 2021.
  * Ioffe and Szegedy [2015] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In _International conference on machine learning_ , page 448–456. PMLR, 2015.
  * Srivastava et al. [2014] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  _Journal of Machine Learning Research_ , 15(1):1929–1958, 2014.
  * Snoek et al. [2012] Jasper Snoek, Hugo Larochelle, and Ryan P. Adams.  Practical bayesian optimization of machine learning algorithms.  In _Proceedings of the 25th International Conference on Neural Information Processing Systems - Volume 2_ , 2012.
  * Yi et al. [2020] Jie Yi, Liwen Ren, Dandan Li, Jie Wu, Wan Li, Guanhua Du, and Jinhua Wang.  Trefoil factor 1 (tff1) is a potential prognostic biomarker with functional significance in breast cancers.  _Biomedicine & Pharmacotherapy_, 124:109827, 2020.
  * Harrell et al. [1982] Jr Harrell, Frank E., Robert M. Califf, David B. Pryor, Kerry L. Lee, and Robert A. Rosati.  Evaluating the Yield of Medical Tests.  _JAMA_ , 247:2543–2546, 1982.
  * Graf et al. [1999] Erika Graf, Claudia Schmoor, Willi Sauerbrei, and Martin Schumacher.  Assessment and comparison of prognostic classification schemes for survival data.  _Statistics in medicine_ , 18 17-18:2529–45, 1999.


# Sub-0.1 degree phase locking of a single-photon interferometer

Vojtěch Švarc  1,* Martina Nováková  1 Michal Dudka  1 and Miroslav Ježek1
1Department of Optics, Faculty of Science, Palacký University, 17. listopadu
12, 77146 Olomouc, Czech Republic<EMAIL_ADDRESS>
###### Abstract

We report a single-photon Mach-Zehnder interferometer stabilized to a phase
precision of 0.05 degrees over 15 hours. To lock the phase, we employ an
auxiliary reference light at a different wavelength than the quantum signal.
The developed phase locking operates continuously, with negligible crosstalk,
and for an arbitrary phase of the quantum signal. Moreover, its performance is
independent of intensity fluctuations of the reference. Since the presented
method can be used in a vast majority of quantum interferometric networks it
can significantly improve phase-sensitive applications in quantum
communication and quantum metrology.

††journal: opticajournal

##  1 Introduction

Interferometers locked to a specific phase are widely used in many fields of
photonics research and technology, such as quantum information processing [1,
2], quantum metrology [3, 4, 5], quantum communication [6, 7, 8, 9, 10], and
tests of fundamental physics [11, 12, 13]. In these applications, precise
phase control of single-photon interferometers and photonic routers is crucial
[8, 12, 14]. For instance, the phase stability determines the amount of
quadrature squeezing required for protocols in quantum metrology [15, 16]. In
quantum cryptography and quantum communication, phase instabilities increase
the bit error rate, which deteriorates the speed and security [17]. Generally,
sub-degree stability is considered to be sufficient for most applications
except for ultra-sensitive ones such as loop-based protocols where the phase
error cumulates each roundtrip [18, 19, 14].

Although in classical interferometry, the sub-degree phase stabilization is
feasible, at the single-photon level it is extremely challenging due to
Poissonian photodetection noise. As a consequence, the typical precision of
single-photon phase locking is a few degrees, and the response is slower than
0.1 s [17, 20]. An efficient strategy to overcome the photodetection noise is
to lock the phase via a bright reference light co-propagating with the signal
photons. To address the signal and the reference individually, they must
differ in a selected parameter. The universal approach, suitable for both
free-space and fiber setups, is to differentiate the signal and the reference
by wavelength [6, 8, 9]. An alternative approach, suitable for free-space
setups only, is to add a small transversal displacement between the signal and
the reference beams [1, 11].

Optimal phase locking should possess the following features: excellent phase
stability, long-term performance, continuous operation, and minimal crosstalk
between the signal and the reference. Most of the works in the field of
single-photon interferometer phase locking did not satisfy more than two of
these requirements [17, 20, 21, 15, 22, 6]. Specifically, most of them are
limited to short-term operation and few-degree stability. Two works
demonstrated continuous operation, reasonably low crosstalk, and 0.2-degree
stability in the short term [8] or a 4-degree stability in the long term [23].
However, achieving all of the desired features simultaneously has not been
reported so far.

Here we demonstrate continuous sub-0.1 degree phase locking of a single-photon
Mach-Zehnder interferometer (MZI) in the dual-wavelength configuration. To
lock the phase, we use spectrally stable 1.5 nW reference light adaptively
corrected for intensity fluctuations, and as a result, we reach exceptional
phase stability of 0.05 degrees for 15 hours demonstrated for the signal.
Moreover, using integrated electro-optic modulators embedded in MZI, we
achieve sub-ns phase switching. The single-photon signal is virtually free of
crosstalk from the reference and exhibits visibility of 99.6$\%$ for 3 nm
spectral width, enabling short-pulse operation. These properties make our
approach directly applicable in many quantum protocols requiring fast photonic
routing or low-latency feedforward operation [14, 13]. Furthermore, the
presented phase stabilization can facilitate demanding loop-based protocols
[24, 25, 18, 14, 19] and enable ultra-stable operation of time-bin quantum
communication protocols [12, 6, 7, 8].

##  2 Methods

The phase lock is commonly carried out by an active feedback loop consisting
of a detector, a proportional-integral-derivative (PID) controller, and a
phase modulator. The principle is following: A phase fluctuation is imprinted
to the detected light as an intensity fluctuation. Subsequently, the produced
electronic signal is evaluated by the PID controller, and an appropriate
voltage is fed into the phase modulator in order to compensate for the initial
phase fluctuation. The precision and speed of this approach depend on noise,
frequency response, and dynamic range of the feedback loop. However, there are
additional challenges one has to count on. In this section, we will focus on
intensity and wavelength-related issues and ways how to suppress them.
Although these issues are not commonly addressed, they can significantly limit
the performance of the active feedback loop, especially while aiming at long-
term sub-degree stability.

Conventionally, the PID controller evaluates the actual phase error by
comparing the intensity detected at a single interferometer output with a
constant setpoint. However, as depicted in Fig. 1(a), this solution cannot
distinguish between intensity and phase fluctuations. Therefore the setpoint
does not correspond to the desired phase if the intensity of the reference
fluctuates. This issue can be overcome by evaluating the setpoint using the
actual intensity of the reference. As shown in Fig. 1(b), this is possible by
monitoring both interferometer outputs, summing their outcomes, and using a
fraction of the sum as the setpoint. Apart from suppressing the influence of
intensity fluctuations, this adaptive setpoint improves the signal-to-noise
ratio (SNR) by a factor of 2 since two independent noisy detectors are
averaged.

Figure 1: Phase locking with (a) constant setpoint and (b) adaptive setpoint.
Phase versus intensity (interference fringe) is plotted for time $t_{1}$ and
$t_{2}$. The phase is set to 90 degrees in $t_{1}$, and then in $t_{2}$, the
total intensity decreases. (a) If the phase is locked using one detector and a
constant setpoint, the phase lock cannot distinguish between intensity and
phase drift. Consequently, decreased intensity in $t_{2}$ leads to a phase
locking error. (b) By detecting both outputs, the setpoint can be adaptively
corrected for the intensity fluctuations. Therefore, the phase is kept
constant despite decreased intensity in $t_{2}$. (c) An extended model of MZI
including LSD parameters $\mu_{\mathrm{O}}$ and $\mu_{\mathrm{A}}$.

Although the adaptive setpoint eliminates efficiently the total intensity
changes, it does not sufficiently respond to local intensity fluctuations
caused by an asymmetric misalignment of the setup (e.g. detuning in a single
arm). This effect, which we call local setup detuning (LSD), is typically much
weaker than the total intensity drift, but it is harder to determine and
correct for. LSD spoils the phase estimation, which consequently leads to a
phase error since the setpoint no longer corresponds to the desired phase. We
divide LSD into two types described by parameters $\mu_{\mathrm{O}}$ and
$\mu_{\mathrm{A}}$ as illustrated in Fig. 1(c).

The first type of LSD is relative detuning between interferometer outputs
which is equivalent to a time-dependent loss $\mu_{\mathrm{O}}$ at a single
output port. The phase estimation $\varphi_{\mathrm{EST}}$ using output
intensities $I_{1}$ and $I_{2}$ is described as

$\varphi_{\mathrm{EST}}=\arccos\left[\dfrac{\dfrac{I_{1}}{V_{2}}-\dfrac{I_{2}}{V_{1}}\left(1-\mu_{\mathrm{O}}\right)}{I_{1}+I_{2}\left(1-\mu_{\mathrm{O}}\right)}\right],$
(1)

where $V_{1}$ and $V_{2}$ are visibilities at the corresponding output ports.
If $\mu_{\mathrm{O}}=0$, $\varphi_{\mathrm{EST}}$ is equivalent to the real
phase $\varphi$. However, if $\mu_{\mathrm{O}}\neq 0$, the phase is estimated
with an error given as $\varphi_{\mathrm{EST}}-\varphi$. For example, let us
assume that $\mu_{\mathrm{O}}=1\%$ is introduced in the reference. Then the
incorrect phase estimation causes a phase error of 0.3 degrees. Here we assume
a high-visibility regime and phase locked close to $\frac{\pi}{2}$ (parameters
relevant to our setup). This example illustrates that even small detuning
would cause non-negligible phase error.

The second type of LSD is described by detuning of a single arm in the
interferometer, equivalent to time-dependent loss $\mu_{\mathrm{A}}$, which
impacts the visibility and intensity distribution between the outputs. As a
consequence, the phase estimation is biased since we have no access to the
actual visibility during the stabilization process. To describe this effect,
it is necessary to express the intensities and visibilities in Eq. (1) using
transmittances $T_{1}$ and $T_{2}$ of beam splitters employed. For simplicity,
let us assume perfect coherence and $\mu_{\mathrm{O}}=0$. Then the
corresponding equation for the estimated phase in MZI reads

$\varphi_{\mathrm{EST}}=\arccos\left[\frac{T_{2}(1-T_{2})(1-2T_{1})\mu_{\mathrm{A}}+2\sqrt{T_{1}T_{2}(1-T_{1})(1-T_{2})(1-\mu_{\mathrm{A}})}\cos\varphi}{2\sqrt{T_{1}T_{2}(1-T_{1})(1-T_{2})}(1-T_{1}\mu_{\mathrm{A}})}\right].$
(2)

From this analysis, the phase estimation error is negligible if
$\mu_{\mathrm{A}}$ is small (less then $10\%$) and the first beam splitter is
balanced perfectly ($T_{1}=0.5$). Unfortunately, fiber couplers show a strong
dependence on wavelength, so the splitting ratio for either the signal or the
reference typically differs from 50:50. We will demonstrate this effect using
parameters relevant to our setup. Let us assume that the reference phase is
locked at $\frac{\pi}{2}$ and the interferometer arm losses are relatively
detuned by $1\%$, which corresponds to $\mu_{\mathrm{A}}=1\%$. This detuning
can potentially impact both the signal and the reference. But in our case, the
fiber couplers are designed for the signal with the wavelength 810 nm, thus
$T_{1,2}\approx 0.5$ and the signal phase is estimated correctly. However, for
the reference with the wavelength 840 nm, the fiber couplers are unbalanced
($T_{1,2}=0.35$). According to Eq. (2), the reference phase is estimated with
a 0.1-degree error. Consequently, the phase is locked with a 0.1-degree error.

So far, we have discussed phase errors caused by intensity fluctuations. Now
let us briefly discuss the impact of wavelength fluctuations. As the phase
$\varphi$ and the wavelength $\lambda$ are fundamentally connected by equation
$\varphi=2\pi\frac{\Lambda}{\lambda}$, $\Lambda$ being optical path
difference, wavelength instabilities in general affect the phase. In the dual-
wavelength configuration, the phase error is present if the phase between the
signal and the reference drifts relatively. The relative phase drift
$\varphi_{\mathrm{s}}-\varphi_{\mathrm{r}}$ for small wavelength drift
$\Delta\lambda$ is described as

$\varphi_{\mathrm{s}}-\varphi_{\mathrm{r}}=2\pi\left(\frac{\Lambda_{\mathrm{s}}}{\lambda_{\mathrm{s}}^{2}}\Delta\lambda_{\mathrm{s}}-\frac{\Lambda_{\mathrm{r}}}{\lambda_{\mathrm{r}}^{2}}\Delta\lambda_{\mathrm{r}}\right),$
(3)

where indices s and r denote the signal and the reference, respectively. To
avoid the relative phase drift, the wavelengths have to be locked, or the
interferometer needs to be balanced perfectly. The latter also requires
precise dispersion compensation since the zero optical path difference is
necessary for both drifting wavelengths.

##  3 Experiment

The experimental setup is depicted in Fig. 2. Single-mode coupled light-
emitting diode (LED) provides an intensity and spectrally stable source of the
reference. The spectral stability is further improved by adding interference
filters with a central wavelength of 840 nm and full width at half maximum
(FWHM) of 21 nm. In the signal ports, arbitrary light within the 3 nm FWHM
around 810 nm can be injected. For active phase locking characterization, we
injected the signal in a single input port. The signal is generated by a
single-mode coupled and spectrally filtered LED, which perfectly simulates
spectral properties of typical quantum photonic sources such as a spontaneous
parametric down-conversion pumped by ultra-short pulses. Before coupling into
the MZI, both the signal and the reference are horizontally polarized.

Figure 2: Experimental setup for an interferometer with active phase locking.
The signal and the reference are injected into the interferometer and co-
propagate through 8 m of fibers and 0.5 m of free space in each arm. The
interferometer is assembled to maximize interference visibility for a broad
spectrum and simultaneously provide fast phase switching. At the
interferometer outputs, the signal and the reference are separated and
detected. The detected reference is processed by a PID controller, and
consequently, the phase is locked via a fiber stretcher. Legend: light-
emitting diode (LED), interference filter (IF), polarizing beam splitter
(PBS), fiber coupler (FC), phase dispersion compensator (PDC), fiber stretcher
(FS), electro-optic phase modulator (EOM), quarter-wave plate (QWP),
photodiode (PD), sensitive photodiode or single-photon detector (SPD),
proportional-integral-derivative controller (PID), analog-to-digital convertor
(ADC).

The signal and the reference are merged at a balanced fiber coupler (FC1) and
co-propagate through the MZI. To achieve simultaneously high visibility and
fast switching capabilities, the MZI combines fiber and free-space elements.
Single-mode polarization-maintaining fiber part of a total length of 8 m
provides inherent spatial and polarization overlap together with the phase
control. Path balancing, polarization filtration, and dispersion compensation
are achieved in 0.5 m of free space. Balanced fiber couplers optimized for 810
nm are used to split and merge the interfering light. Their splitting ratios
are estimated to be 51:49 for the signal and 35:65 for the reference. The
splitting ratio of FC1 can be effectively modified to perfect 50:50 operation
by adjusting losses in one MZI arm. The arm lengths are balanced in the free
space by micrometric translation of one coupling lens. Thanks to the 21 nm
broad spectrum of the reference, we can balance MZI arms within approximately
two wavelengths, i.e., with phase precision of $\pm 2\pi$. To improve
polarization degree of light, two pairs of polarizing beam splitters (PBSs)
with extinction greater than 30 dB are used. One pair of PBSs is placed in the
free-space part, and the other pair is included in the wavelength separation
part. The MZI allows for fast phase switching via 10 GHz integrated electro-
optic modulators (EOMs). This capability is crucial in applications of the
MZI, however, it significantly increases chromatic dispersion in the MZI. To
cancel out the dispersion, we compose MZI arms symmetrically. Namely, the
fiber lengths are balanced with mm precision, and similar EOMs are placed in
both arms. Additional dispersion compensation and manipulation are performed
in free space with a custom component called phase dispersion compensator
consisting of two elements, PDC1 and PDC2. As shown in Fig. 3(a), the PDC1
consists of two 1 cm thick high-dispersion SF10 glass plates mounted on
contra-directional tilting stages. The PDC1 allows for manipulating the
relative optical path between the wavelengths without any misalignment of
their fiber coupling. Thus we can set an arbitrary relative phase between the
signal and the reference. This enables the reference to be locked at the most
sensitive point of the interference fringe (phase $\frac{\pi}{2}$) regardless
of the signal phase. PDC2 consists of a fixed SF10 glass plate with 3 cm
thickness to compensate for the residual dispersion. As a result, we reach the
visibility of 99.6$\%$ for both MZI outputs for the 3 nm broad signal and
98.6$\%$ for one output of the 21 nm broad reference.

Figure 3: (a) Phase dispersion compensator PDC1 used for tuning the relative
phase between the signal and the reference. (b) Fiber stretcher used for the
phase locking. By applying voltage to a piezoelectric crystal, the fiber
coiled around two aluminium segments is stretched and the phase is changed.

At each MZI output, the wavelengths are demultiplexed via a dichroic optical
isolator assembled from a PBS, a quarter-wave plate (QWP), and an interference
filter. We reach negligible crosstalk of 30 photons/s from the reference to
the signal, equivalent to extinction better than 80 dB. The separation occurs
at the narrow-band interference filter centered at 810 nm, where the signal
passes through while the reference is reflected back. Subsequently, the
reference is reflected at the PBS. This is possible since the reference
polarization is transformed from horizontal to vertical by bi-directional
passing through a 45-degree oriented QWP and $\pi$ phase shift caused by
reflection at the interference filter. Since the QWP undesirably changes the
signal polarization, we use a compensatory QWP after the wavelength separation
part. Before the detection, both the signal and the reference are coupled into
a single-mode fiber. For the signal, polarization-maintaining fibers are used,
which facilitates its transfer to following experimental stages. During active
phase locking tests, the signal with the power of about 1 pW ($4\cdot 10^{6}$
photons/s) impinges directly the ultra-sensitive silicon photodiodes with
noise equivalent power (NEP) of 1.4 fW/$\mathrm{\sqrt{Hz}}$ and 30 Hz
bandwidth. For single-photon measurements, we use single-photon avalanche
diodes with a 50 Hz dark count rate plugged into a time tagger.

The active phase lock is carried out by a feedback loop composed of the
photodiodes, PID controller, and a fiber stretcher. The reference with the
power of 1.5 nW impinges silicon photodiodes with NEP=9
fW/$\mathrm{\sqrt{Hz}}$ and 2 kHz bandwidth. The resulting electronic signals
are processed by an in-house developed analog PID controller based on low-
noise operational amplifiers. The circuit is provided with the adaptive
setpoint described in Methods. Therefore, we lock the phase regardless of the
intensity fluctuations of the reference. The PID controller is fine-tunable in
each parameter via 10-turn precision potentiometers. As a result, the target
phase is set with high accuracy. To provide correct phase estimation for the
adaptive setpoint, we include additional tunability to cancel out detector
offsets and to balance detection efficiencies. The PID controller is designed
to drive a custom-made fiber stretcher depicted in Fig. 3(b). It has 7.2
$\mu$F capacitance, 1 kHz bandwidth, 0.11 V half-wave voltage, and dynamic
range of $\pm 275$ $\mathrm{\mu m}$ [26]. Since the extent of the phase drift
is much smaller in our case, we use only $\pm 110$ $\mathrm{\mu m}$
corresponding to the output voltage range of the PID controller.

To reduce air flux and temperature instabilities, we cover the whole setup
with styrofoam plates during all measurements. The passive protection reduces
the temperature drift from $2^{\circ}$C peak-to-peak per 2 hours to a peak-to-
peak drift of $0.3^{\circ}$C per 15 hours. Our setup is placed on a 60x60x6 cm
breadboard with no additional damping of mechanical vibrations. For the data
acquisition, we use an 8-channel 18-bit analog-to-digital converter set to a
sampling frequency of 16 Hz. We increase SNR by averaging the data samples to
1 Hz except for the detailed spectral characterization (Fig. 4(b)). We collect
data at each detector, and additionally, we monitor the voltage applied to the
fiber stretcher. The knowledge of the driving voltage, together with the
information on half-wave voltage, allows for the evaluation of the compensated
phase drift. The collected data are processed by an algorithm based on Eq.
(1). Apart from the output intensities, it requires additional parameters such
as the visibility of each MZI output and a ratio of losses in the MZI outputs.
For measurements with a non-stabilized phase, these parameters are
periodically recovered as interference fringes are scanned spontaneously due
to the phase drift. While the stabilization process is running, interference
fringe scanning is not possible, thus the additional parameters are
independently measured before the stabilization starts.

##  4 Results and discussion

We compare the performance of the MZI with and without the active
stabilization and present the phase stability in three ways: in the time
domain, in the frequency domain, and via the Allan deviation. The results in
the time domain are depicted in Fig. 4(a). The typical phase drift of non-
stabilized MZI is shown on the left side. On the right side, we show phase
noise for stabilized MZI (orange and red lines) and the phase compensated
during the stabilization process (black line). For further analysis, we use
only black, orange, and red datasets since they were achieved within the same
environmental conditions. While the phase lock is active, a phase drift of
4000 degrees per 15 hours is virtually eliminated. In terms of standard
deviation accumulated over a typical single-photon sampling frequency of 1 Hz,
we reach 0.05 degrees for the signal and 0.002 degrees for the reference
within a 15-hour interval.

Further, we analyze the spectral power density of the phase noise shown in
Fig. 4(b). During the stabilization process, the phase noise is suppressed by
9 orders of magnitude for frequencies in the range $10^{-3}$ to $10^{-5}$ Hz.
For frequencies above $10^{-3}$ Hz, the phase noise is primarily given by
detection noise. Furthermore, we evaluate the Allan deviation of the phase
noise depicted in Fig. 4(c). Allan deviation corresponds to root mean square
phase drift for a variable time interval $\tau$ [27]. For the signal, we reach
Allan deviation below $2\times 10^{-2}$ degrees for time intervals from 10 s
to more than 5 hours. The minimum value $2\times 10^{-3}$ degrees is achieved
for $10^{3}$ s. Allan deviation of the reference is evaluated to be smaller
than $10^{-3}$ degrees.

Figure 4: Comparison of stabilized and non-stabilized MZI. In stabilized MZI,
the reference and the signal are represented by the red line and the orange
line, respectively. (a) Stability in the time domain. The blue line represents
the phase in non-stabilized MZI and the black line corresponds to the phase
drift compensated during the phase locking (see the right vertical axis). (b)
Spectral power density $G$ of the phase noise. (c) Allan deviation of the
phase noise. The gray line illustrates the noise level of the signal
acquisition.

The results show that the reference exhibits higher phase precision than the
signal. This discrepancy has two reasons. Firstly, the signal is affected by a
higher amount of detection noise than the reference. Specifically, the
detection noise superimposes the signal up to $10^{3}$ s as illustrated by the
gray line in Fig. 4(c). Secondly, the signal is slightly affected by the
residual LSD leading to increased phase error in the long term. Considering
the aforementioned reasons, the Allan deviation of the signal should be
interpreted as an upper bound of the real phase error in the MZI. On the other
hand, the Allan deviation of the reference gives the lower bound of the real
phase error since it does not reflect LSD. We estimate that in the short term,
where LSD is negligible, the real Allan deviation follows the reference curve,
whereas, in the long term, it approaches the signal curve. The measured phase
precision of the reference illustrates the ultimate performance of the phase
lock. Therefore, the direct demonstration of sub-$10^{-3}$ degree Allan
deviation for time intervals up to $10^{4}$ s is feasible provided that we
completely cancel out LSD and improve SNR of the signal detection.

Figure 5: Complete phase tunability via the phase dispersion compensator
demonstrated for a classical signal (a) and single photons (b). Although the
reference is locked at the most sensitive phase $\frac{\pi}{2}$, we are able
to tune arbitrary splitting ratio $T$:$R$ for the signal.

Although, in general, the performance of phase locking is dependent on the
reference phase (the best performance is achieved in $\pi/2$), we are not
limited to a specific phase value of the signal. In Fig. 5, we show the phase
locking for an arbitrary signal phase while being locked with the reference
close to $\pi/2$. This is possible due to full phase tunability between the
signal and the reference achieved by the phase dispersion compensator element
(see Experiment for details). Fig. 5(a) shows the full tunability for a
classical signal, while in Fig. 5(b) we demonstrate the tunability for the
single-photon signal (count rate approx. $5\times 10^{5}$ counts/s,
integration time 1 s). Additional results at the single-photon level were
achieved previously with preliminary versions of the setup. Firstly, we
demonstrated balanced time-bin multiplexing for photon-number resolving
detection using a loop configuration [14]. Secondly, we demonstrated the Fock
state conversion exploiting feedforward control of the interferometer [13].
These experiments were performed before a significant improvement of the phase
lock, hence they do not reflect the unprecedented phase stability presented
here.

Let us compare our results with other state-of-the-art experiments across
various interferometer configurations and target applications. Niwa et al.
reported 20 pm phase noise between 0.3 mHz and 1 Hz in 5x5 cm free-space MZI
[28]. The results are achieved for strong light and classical detectors. They
use advanced passive and active methods of protection against environment
changes: the setup is placed in a vacuum chamber with 1 mK thermal stability,
and the MZI elements are composed of ultralow expansion glass. In comparison,
we reach phase noise around 100 pm within a much larger frequency interval
(0.02 mHz to 1 Hz) for an 8 m long fiber interferometer placed in a standard
lab environment. In the field of single-photon experiments, the phase
stability is typically in order of a few degrees [17, 20, 21, 15, 22, 6, 29,
23]. Five of the aforementioned works are limited to short-term operation (up
to few minutes)[20, 21, 15, 22, 6], whereas long-term stability (up to several
hours) is achieved only in three works [17, 29, 23]. To the best of our
knowledge, the most precise stabilization of a single-photon interferometer
has been reported by Toliver et al. [8] and Roztocki et al. [30]. Toliver et
al. achieved a 200 s phase lock with a standard deviation of $0.2$ degrees. In
contrast, we demonstrate a standard deviation of 0.05 degrees for more than 15
hours. Roztocki et al. reached Allan deviation below 0.2 degrees from
$10^{-1}$ s to $10^{4}$ s, but these results were demonstrated only for a
bright reference light (the phase was not reconstructed using an independent
signal as in our case). In contrast, we reach Allan deviation below $5\times
10^{-3}$ degrees for the reference within the same interval.

Our results demonstrate the unprecedented performance of dual-wavelength phase
locking across single-photon interferometers. Our approach will work for all
types of interferometers and complex interferometric networks, including fiber
and free-space propagation. Regarding free-space interferometers, phase
locking might be even simpler, as a small spatial displacement between the
signal and the reference can substitute the wavelength separation. However,
this scheme is not free from the relative phase drift between the signal and
the reference since the index of refraction in air locally changes for each of
the displaced paths [29]. Further, our phase-locking method can be utilized
for unbalanced interferometers, provided that a sufficiently coherent source
is employed as the reference. However, highly coherent sources can bring
additional challenges due to random interference on optical surfaces and
elements, leading to intensity fluctuation throughout the setup.

##  5 Conclusion

We have implemented a novel technique of continuous phase locking for single-
photon circuits and networks. Our approach extends dual-wavelength
stabilization methods by the adaptive setpoint, enabling precise phase locking
despite the fluctuating intensity. We show that the adaptive setpoint,
together with thorough optimization in all possible degrees of freedom, leads
to unprecedented stability even with faint optical signals. We have
experimentally demonstrated sub-0.1 degree phase locking of 8.5 m long fiber-
based MZI for more than 15 hours. Specifically, we have reached a phase
standard deviation of 0.05 degrees, which is equivalent to the path
displacement of about 100 pm. To the best of our knowledge, our results
represent the most precise phase locking of single-photon interferometers.
Moreover, we have reached complete phase tunability of the device and
negligible crosstalk from the reference to the single-photon signal. Since our
method can be implemented in a vast majority of interferometric networks, it
can significantly improve the performance of many phase-sensitive applications
in information processing, quantum metrology, and quantum communication.

## Funding

Czech Science Foundation (project 19-19189S); Ministry of Education, Youth,
and Sports of the Czech Republic (projects
CZ.02.2.69/0.0/0.0/19$\\_$073/0016713, DSGC-2021-0118).

## Acknowledgments

We thank Petr Novák for his help with the 3D printing of the phase dispersion
compensator mount. Student grant DSGC-2021-0118 is funded under the OPIE
project „Improvement of Doctoral Student Grant Competition Schemes and their
Pilot Implementation“, reg. no. CZ.02.2.69/0.0/0.0/19$\\_$073/0016713 in the
grant results.

## Disclosures

The authors declare that there are no conflicts of interest related to this
article.

## Data Availability

Data underlying the results presented in this paper are available in Ref.
[31].

## References

  * [1] X. song Ma, S. Zotter, N. Tetik, A. Qarry, T. Jennewein, and A. Zeilinger, “A high-speed tunable beam splitter for feed-forward photonic quantum information processing,” Opt. Express 19, 22723–22730 (2011).
  * [2] M. Miková, H. Fikerová, I. Straka, M. Mičuda, J. Fiurášek, M. Ježek, and M. Dušek, “Increasing efficiency of a linear-optical quantum gate using electronic feed-forward,” Phys. Rev. A 85, 012305 (2012).
  * [3] V. Giovannetti, S. Lloyd, and L. Maccone, “Advances in quantum metrology,” Nature Photonics 5, 222–229 (2011).
  * [4] J. Aasi, J. Abadie, B. Abbott, R. Abbott, T. Abbott, M. Abernathy, C. Adams, T. Adams, P. Addesso, R. Adhikari _et al._ , “Enhanced sensitivity of the LIGO gravitational wave detector by using squeezed states of light,” Nature Photonics 7, 613–619 (2013).
  * [5] G. Thekkadath, M. Mycroft, B. Bell, C. Wade, A. Eckstein, D. Phillips, R. Patel, A. Buraczewski, A. Lita, T. Gerrits _et al._ , “Quantum-enhanced interferometry with large heralded photon-number states,” NPJ Quantum Information 6, 1–6 (2020).
  * [6] G. B. Xavier and J. P. von der Weid, “Stable single-photon interference in a 1 km fiber-optic Mach–Zehnder interferometer with continuous phase adjustment,” Opt. Lett. 36, 1764–1766 (2011).
  * [7] G. Carvacho, J. Cariñe, G. Saavedra, Á. Cuevas, J. Fuenzalida, F. Toledo, M. Figueroa, A. Cabello, J.-Å. Larsson, P. Mataloni, G. Lima, and G. B. Xavier, “Postselection-loophole-free Bell test over an installed optical fiber network,” Phys. Rev. Lett. 115, 030503 (2015).
  * [8] P. Toliver, J. M. Dailey, A. Agarwal, and N. A. Peters, “Continuously active interferometer stabilization and control for time-bin entanglement distribution,” Opt. Express 23, 4135–4143 (2015).
  * [9] S.-B. Cho and H. Kim, “Active stabilization of a fiber-optic two-photon interferometer using continuous optical length control,” Opt. Express 24, 10980–10986 (2016).
  * [10] Y. Xu, J. Lin, Y. Li, H. Dai, S. Liao, and C. Peng, “Active Phase stabilization for the interferometer with 128 actively selectable paths,” IEEE Trans. Nucl. Sci. 66, 1076–1080 (2019).
  * [11] X.-s. Ma, S. Zotter, J. Kofler, R. Ursin, T. Jennewein, Č. Brukner, and A. Zeilinger, “Experimental delayed-choice entanglement swapping,” Nature Physics 8, 479–484 (2012).
  * [12] F. Vedovato, C. Agnesi, M. Tomasin, M. Avesani, J.-A. Larsson, G. Vallone, and P. Villoresi, “Post-selection-loophole-free Bell violation with genuine time-bin entanglement,” Physical Review Letters 121, 190401 (2018).
  * [13] V. Švarc, J. Hloušek, M. Nováková, J. Fiurášek, and M. Ježek, “Feedforward-enhanced Fock state conversion with linear optics,” Opt. Express 28, 11634–11644 (2020).
  * [14] V. Švarc, M. Nováková, G. Mazin, and M. Ježek, “Fully tunable and switchable coupler for photonic routing in quantum detection and modulation,” Opt. Lett. 44, 5844–5847 (2019).
  * [15] Y. Takeno, M. Yukawa, H. Yonezawa, and A. Furusawa, “Observation of -9 dB quadrature squeezing with improvement of phase stability in homodyne measurement,” Opt. Express 15, 4321–4327 (2007).
  * [16] T. Eberle, V. Händchen, and R. Schnabel, “Stable control of 10 dB two-mode squeezed vacuum states of light,” Opt. Express 21, 11546–11553 (2013).
  * [17] V. Makarov, A. Brylevski, and D. R. Hjelme, “Real-time phase tracking in single-photon interferometers,” Appl. Opt. 43, 4385–4392 (2004).
  * [18] S. Takeda and A. Furusawa, “Universal quantum computing with measurement-induced continuous-variable gate sequence in a loop-based architecture,” Physical Review Letters 119, 120504 (2017).
  * [19] S. Takeda, K. Takase, and A. Furusawa, “On-demand photonic entanglement synthesizer,” Sci. Adv. 5, eaaw4530 (2019).
  * [20] D. Pulford, C. Robillard, and E. Huntington, “Single photon locking of an all-fiber interferometer,” Rev. Sci. Instrum. 76, 063114 (2005).
  * [21] S. Suzuki, H. Yonezawa, F. Kannari, M. Sasaki, and A. Furusawa, “7 dB quadrature squeezing at 860 nm with periodically poled K Ti OPO4,” Appl. Phys. Lett. 89, 061116 (2006).
  * [22] S.-B. Cho and T.-G. Noh, “Stabilization of a long-armed fiber-optic single-photon interferometer,” Opt. Express 17, 19027–19032 (2009).
  * [23] S. Yanikgonul, R. Guo, A. Xomalis, A. N. Vetlugin, G. Adamo, C. Soci, and N. I. Zheludev, “Phase stabilization of a coherent fiber network by single-photon counting,” Opt. Lett. 45, 2740–2743 (2020).
  * [24] K. R. Motes, A. Gilchrist, J. P. Dowling, and P. P. Rohde, “Scalable boson sampling with time-bin encoding using a loop-based architecture,” Phys. Rev. Lett. 113, 120501 (2014).
  * [25] Y. He et al., “Time-bin-encoded boson sampling with a single-photon device,” Phys. Rev. Lett. 118, 190501 (2017).
  * [26] M. Nováková, L. Podhora, V. Švarc, and M. Ježek, “Polarization-maintaining 3D printed fiber stretcher,” in preparation (2022).
  * [27] D. W. Allan, “Statistics of atomic frequency standards,” Proceedings of the IEEE 54, 221–230 (1966).
  * [28] Y. Niwa, K. Arai, A. Ueda, M. Sakagami, N. Gouda, Y. Kobayashi, Y. Yamada, and T. Yano, “Long-term stabilization of a heterodyne metrology interferometer down to a noise level of 20 pm over an hour,” Appl. Opt. 48, 6105–6110 (2009).
  * [29] M. Mičuda, E. Doláková, I. Straka, M. Miková, M. Dušek, J. Fiurášek, and M. Ježek, “Highly stable polarization independent Mach-Zehnder interferometer,” Rev. Sci. Instrum. 85, 083103 (2014).
  * [30] P. Roztocki, B. MacLellan, M. Islam, C. Reimer, B. Fischer, S. Sciara, R. Helsten, Y. Jestin, A. Cino, S. T. Chu _et al._ , “Arbitrary phase access for stable fiber interferometers,” Laser & Photonics Reviews 15, 2000524 (2021).
  * [31] V.Švarc, “GitHub repository SvarcOpticsOlomouc/stabilization,” https://github.com/SvarcOpticsOlomouc/stabilization.