# Spatial Transcriptomics Dimensionality Reduction
using Wavelet Bases

Zhuoyan Xu
Department of Statistics
University of Wisconsin - Madison
<EMAIL_ADDRESS>
&Kris Sankaran
Department of Statistics
University of Wisconsin - Madison
<EMAIL_ADDRESS>

###### Abstract

Spatially resolved transcriptomics (ST) measures gene expression along with
the spatial coordinates of the measurements. The analysis of ST data involves
significant computation complexity. In this work, we propose gene expression
dimensionality reduction algorithm that retains spatial structure. We combine
the wavelet transformation with matrix factorization to select spatially-
varying genes. We extract a low-dimensional representation of these genes. We
consider Empirical Bayes setting, imposing regularization through the prior
distribution of factor genes. Additionally, We provide visualization of
extracted representation genes capturing the global spatial pattern. We
illustrate the performance of our methods by spatial structure recovery and
gene expression reconstruction in simulation. In real data experiments, our
method identifies spatial structure of gene factors and outperforms regular
decomposition regarding reconstruction error. We found the connection between
the fluctuation of gene patterns and wavelet technique, providing smoother
visualization. We develop the package and share the workflow generating
reproducible quantitative results and gene visualization. The package is
available at https://github.com/OliverXUZY/waveST.

_K_ eywords Spatial Transcriptomics $\cdot$ Wavelet Transformation $\cdot$
Empirical Bayes Matrix Factorization $\cdot$ Factor Gene

##  1 Introduction

Spatial resolved transcriptomics (ST) is a technology measuring spatial
variation in gene expression. Different technological platforms support
different genome coverage and spatial resolution (tissue-level measurement to
subcellular measurement). Examples of ST platforms include Spatial
Transcriptomics [Ståhl et al., 2016; Xia et al., 2019], 10x Genomics Visium,
Slide-seq [Rodriques et al., 2019], sci-Space [Srivatsan et al., 2021], and
MERFISH [Chen et al., 2015]. Spatial transcriptomics allows visualization and
quantitative analysis of the transcriptome with spatial resolution in
individual tissue sections. A series of studies combining gene expression and
spatial information has been brought to generate new insight in biological
analysis [Kuppe et al., 2020; Shah et al., 2016; Berglund et al., 2018].

Quantification of gene expression has wide applications in transcriptomics.
Understanding the spatial distribution of gene expression has helped to answer
fundamental questions in developmental biology [Asp et al., 2019; Rödelsperger
et al., 2021], cancer [Thrane et al., 2018; Moncada et al., 2020], and
neuroscience [Moffitt et al., 2016; Close et al., 2021]. Two widely-used
methods for gene expression quantification are fluorescent in situ
hybridization (FISH) and next-generation sequencing.

While these approaches have been made to measure gene expression while
preserving spatial information, there are statistical challenges in analyses
combining the gene and spatial information. Specifically, gene dimensionality
reduction guided by spatial information is still an active area of study. In
this work, we conduct a dimensionality reduction on spatially varying gene
data. After transformation, gene expression over locations is a square matrix,
which allows us to view it as an image. Therefore, tools from image processing
can be adapted – our method incorporates spatial information by wavelet
transformation, a multi-scale analysis decompose sequence into orthonormal
series. We apply wavelet transformation for each gene expression over
locations. This is a common technique in denoising the image. We thus apply
techniques analogous to image analysis. We evaluate performance using
reconstruction error and comparing resulting visualizations. Our analysis
pipeline build from off-the-shelf tools, and code is available for reproducing
our results. We also attach visualizations in our pipeline, giving example
interpretations of intermediate result.

A natural idea in spatial transcriptomic analysis is to identify gene types
across spatial locations [Abu-Jamous and Kelly, 2018]. For example, in
clustering the representations in each cluster achieve dimensionality
reduction under the assumption that genes within one cluster have the same
type. A parallel technique is clustering based on cell type instead. General
clustering methods can be combined into more sophisticated pipelines tailored
toward spatial single-cell analysis. SC3 [Kiselev et al., 2017] is an ensemble
clustering method. It calculates distance matrices across cell locations using
the Euclidean distance, then applies spectral clustering, and assigns
membership [Kiselev et al., 2018]. To identify the cell type in each cluster,
one can perform differential expression analysis between all pairs of
clusters. scGeneFit uses a label-aware compression method to find marker genes
[Dumitrascu et al., 2021]. Given the cell-by-gene expression matrix and cell
clustering membership, scGeneFit maps cells to lower-dimensional space where
cells within the same cluster are closer. For gene dimensionality reduction,
[Zhu and Sabatti, 2020] constructs a neighborhood graph from the spatial
coordinates, then applies a graph-based feature selection procedure to
determine spatially varying genes. They also provide the option to infer a
latent graph embedding for cells based on selected genes, applying spline
models to fit the gene’s expression on the latent embedding. Then they
leverage the fitted coefficients to reduce the dimensionality of each gene.
[Svensson et al., 2018] proposed pipelines using mixed-effect models
incorporating spatial information. The model contains two random effect terms:
a spatial variance term that parametrizes gene expression covariance by
pairwise distance between samples and a noise term that models nonspatial
variability. The model leverages efficient inference methods previously
developed for linear mixed models, and it is computationally efficient.

Our setting is closely aligned with the following recent works: [Shang and
Zhou, 2022] developed SpatialPCA, applying probabilistic PCA on ST data for
dimensionality reduction. They assume data are given as a location-by-genes
matrix and construct a regression model similar to factor analysis, where the
prior covariance matrix of factor genes is a distance matrix constructed with
a Gaussian kernel. [Velten et al., 2022] proposed MEFISTO, combining factor
analysis with the non-parametric framework of Gaussian processes to model
spatio-temporal dependencies in the latent space. [Townes and Engelhardt,
2021] developed nonnegative spatial factorization (NSF), combining a Gaussian
process prior over spatial locations and a Poisson or negative binomial
likelihood for count data, identifying generalizable spatial patterns of gene
expression. All these works impose spatial structure on the prior of the
factor genes. While these methods offer new dimensionality reduction
techniques to cluster the genes, the complex model structure and a large
number of hyperparameters introduce uncertainty and noise. Instead of imposing
structural assumptions on the prior of the factor genes, we impose structure
on the factor gene itself. Our contributions are the following:

  * •

We propose an approach, based on techniques from matrix decomposition and
image signal processing to perform gene dimensionality reduction that retains
inferred spatial structure.

  * •

We run simulations showing that wavelet-guided dimensionality reduction
performs better estimation than the singular value decomposition (SVD) under
low signal-to-noise (SNR) regime.

  * •

We perform real data experiments to identify the connection between wavelet
techniques and fluctuation of the gene expression, which would be useful in
selecting spatially related genes based on reconstruction error.

  * •

We provide a gene extraction pipeline capturing the global information of
spatially related genes. We provide smoother visualization of gene factors via
wavelet methods. We develop an R package waveST and share the workflow
generating reproducible quantitative results and gene visualization.

The diagram for workflow can be seen in Figure 1. The paper is organized in
following: In Section 2 we introduce background on required techniques,
including the wavelet transformation and matrix decompositions. In Section 3,
we formally define our problem under this setting, and section 4 introduces
our algorithms and analysis pipeline. In section 5, we implement simulations
showing the effect of wavelet transformation in reducing error. In section 6
we conduct our method on data from [Weber, 2021], showing the reconstruction
error in dimension reduction and visualization of lower-dimensional
representations.

Figure 1: A summary of the proposed workflow.

##  2 Background

In this section, we will cover background of methods we used in data pre-
processing and analysis.

###  2.1 Wavelet Transformation

Density estimation and function approximation is a fundamental problem in
statistics and machine learning. Non-parametric methods such as spline
regression [Perperoglou et al., 2019], Fourier transformation [Cochran et al.,
1967] and Wavelet transformation [Nason, 2008] has been used in such
scenarios.

Consider a model with one predictor: $y=f(x)+\epsilon$, where
$\mathbb{E}\epsilon^{2}=\sigma^{2}$. We want to estimate $f$, the _trend_ of
the function. We assume the predictors are ordered $x_{1}<x_{2}<\ldots<x_{n}$.
We have a signal or frequency to estimate.

Consider the _Haar mother_ wavelet:

$\psi(x)=\left\\{\begin{array}[]{cc}1&x\in\left[0,\frac{1}{2}\right)\\\
{}\hfil\\\ -1&x\in\left[\frac{1}{2},1\right)\\\ {}\hfil\\\ 0&\text{ otherwise
}\end{array}\right.$

satisfying $\int_{-\infty}^{\infty}\psi(x)dx=0$. Unlike the Fourier basis,
wavelets oscillate and decay fast, only contributing to a certain local area
and zero elsewhere. We can generate wavelets from the Haar mother wavelet
$\psi_{j,k}(x)=2^{j/2}\psi\left(2^{j}x-k\right)$ for integer $j,k$. These
wavelets form orthonormal sets. We can decompose the trend as
$f(x)=\sum_{j=-\infty}^{\infty}\sum_{k=-\infty}^{\infty}d_{j,k}\psi_{j,k}(x)$,
where $d_{j,k}=\int_{-\infty}^{\infty}f(x)\psi_{j,k}(x)dx=<f,\psi_{j,k}>$ are
the wavelet coefficients.

Wavelet based methods have several advantages among other non-parametric
methods, especially dealing with sparse data. One property wavelet has is
localization. If a sequence has a discontinuity, this will only influence the
wavelet basis around it. In contrast, for a Fourier basis consisting of sine
and cosine functions at different frequencies, every basis element will
interact with this discontinuity, hence influencing every Fourier coefficient.

The simplest discrete wavelet transformation calculates the difference and
sums between each adjacent pair. Suppose we have one vector with length $n$,
where $n$ is dyadic ($n=2^{J}$). We computed $d_{J-1,k}=y_{2k}-y_{2k-1}$ as
finest-level detail and $c_{J-1,k}=y_{2k}+y_{2k-1}$ as the finest-level
averages. We have $\left\\{d_{k}\right\\}$ containing $n/2=2^{J-1}$ as our
finest level coefficients. To obtain the next coarsest coeffients we set:

$\displaystyle d_{J-2,\ell}$ $\displaystyle=c_{J-1,2\ell}-c_{J-1,2\ell-1}$ (1)
$\displaystyle c_{J-2,\ell}$ $\displaystyle=c_{J-1,2\ell}+c_{J-1,2\ell-1}$ (2)

We continue this procedure until $j$ reaches one. The set of details
$\left\\{d_{k}\right\\}$ across levels are our wavelet coefficients.

We use $J$ to denote the scale of wavelets, with larger $J$ we have finer
scale and approximation. However, finer scale sometimes introduces more
parameters to capture minor details of the sequence (overfitting). The trade-
off in choosing $J$ is crucial in the wavelet transformation.

Regularization and smoothing is can be used to prevent overfitting.
Regularization is usually conducted by shrinking wavelet coefficients. The
concept of wavelet shrinkage was first proposed by [Donoho and Johnstone,
1994]. The motivation behind wavelet shrinkage is straightforward. Consider
empirical wavelet coefficients, the large coefficients usually contain true
signal and noise, whereas the small coefficients only contain noise. The
shrinkage is often used when the wavelet coefficients are assumed to be a
_sparse vector_.

Wavelet shrinkage is often conducted by setting a _threshold_ (denoted by
$\delta$) and only keeping the coefficients above the threshold. To choose a
threshold, a natural metric is the squared error between estimated function
and the truth: $\hat{M}=\sum_{i=1}^{n}(f(x_{i})-\hat{f}(x_{i}))^{2}$ and
$M=\mathbb{E}\hat{M}$. Donoho et al. [1994] proposed a _universal threshold_
$\delta=\sigma\sqrt{2\log n}$, which induced $M\leq\mathcal{O}(\log
n\sigma^{2})$. Donoho and Johnstone [1995] also proposed _Sure threshold_
based on Stein’s (1981) unbiased risk estimator. The optimal SURE threshold
can be obtained in $\mathcal{O}(n\log n)$ operations. Donoho and Johnstone
[1995] also noted that SURE sometimes failed when the true signal coefficients
are highly sparse, they proposed a hybrid scheme called _SureShrink_ ,
combining the SURE and universal thresholds, using them depending on certain
situations.

The extension of wavelet methods to 2D regularly spaced data (images) and such
data in higher dimensions was proposed by [Mallat, 1989]. We only consider 2-D
wavelet transformation since we decompose 2-D spatial gene expression data.
Suppose we have $n\times n$ matrix $A$ where $n=2^{J}$ is dyadic. A simple
discrete wavelet transformation on $A$ first applies procedure (1) and (2) to
the rows of the matrix. We then have two matrices of size
$n\times\frac{n}{2}$, called $H$ and $G$. Then we apply the same procedures to
both the columns of $H$ and $G$, resulting in four matrices $HH,GH,HG$, and
$GG$ each of size $\frac{n}{2}\times\frac{n}{2}$. These are our finest level
coefficients. $HH$ is the local average of the original matrix used for the
next level procedure.

###  2.2 Factor Gene

Clustering and dimensionality reduction are widely used in genomics. In a
gene-by-sample matrix, genes are often grouped into profiles, where genes from
the same profile have a similar function. Statistically, we can treat them as
correlated variables. We use the term factor genes as the principal component
in our gene-by-sample data. Factor genes compose the linear combinations of
genes. We consider each gene as a variable, and we aim to find variables
$(\boldsymbol{f_{1}},\ldots,\boldsymbol{f_{K}})$ such that each gene $g_{i}$
has form
$\boldsymbol{g_{i}}=a_{1}\boldsymbol{f_{1}}+a_{2}\boldsymbol{f_{2}}+\cdots
a_{K}\boldsymbol{f_{K}}$, where $a_{j}$ are coefficients.

Suppose we have gene($P$)-by-sample($N$) matrix $A\in\mathbb{R}^{N\times P}$
with sample covariance matrix $S=\frac{1}{n}\tilde{A}^{T}\tilde{A}$, where
$\tilde{A}$ is the column centered $A$. Consider the SVD on
$\tilde{A}=U\Lambda V^{T}$. Then we have $SV=V\Lambda^{2}$. The columns of $V$
are called _eigenarrays_ [Alter et al., 2000]. The first $K$ columns of $Z=AV$
are the _factor genes_. The factor genes capture the mutual underlying
information of genes.

###  2.3 Empirical Bayes Matrix Factorization

Matrix Factorization is often used in capture factor genes. We have a
formulation similar to [Wang and Stephens, 2021], consider the factorization
model on observed samples by gene data $Y\in\mathbb{R}^{N\times P}$

$Y=FL^{T}+E$

where $F\in\mathbb{R}^{N\times K}$ denotes the factors,
$L\in\mathbb{R}^{P\times K}$ denotes the loadings for each factor, and
$E\in\mathbb{R}^{N\times P}$ denotes Gaussian noise with zero mean. Typical
formulations would treat factors and loadings as fixed effects and use Maximum
Likelihood Estimation (MLE). In our high-dimensional setting, penalty-based
regularizers are often considered. Considering a prior for $L$ and $F$ under
Bayesian setting has a similar effect and several advantages over adding a
regularizer alone, such as simplifying hyperparameter search and selection of
the number of factors $K$.

One feature of empirical Bayes approaches in matrix factorization proposed by
[Bishop, 1999] is that the methods automatically select the number of factors
$K$. In [Wang and Stephens, 2021], the author add factors with prior one at a
time, estimating priors at each step until convergence. If the computed prior
of the newly added factor is almost point mass on 0, the algorithm eliminates
this factor and returns.

Following [Wang and Stephens, 2021], we consider _Empirical Bayes_ in our
setting, where we set a prior with unknown parameters of $L$ and $F$.
Empirical Bayes is not strictly Bayes, since the prior parameters are directly
estimated by the MLE of the data. One example would be normal distribution
$\mathcal{D}_{F}$ of $F$ and $\mathcal{D}_{L}$ of $L$ where coordinates are
independent, which is conjugate prior. We estimate parameters in
$\mathcal{D}_{F}$ and $\mathcal{D}_{L}$ by the maximize the marginal
likelihood calculated by integral out $L$ and $F$. Then we computed the
posterior distribution of $L$ and $F$.

##  3 Problem Setup

Consider $Y\in\mathbb{R}^{N\times P}$, where each row represents a sample and
each column represents a gene. Further, let $S\in\mathbb{R}^{N\times 2}$ store
spatial coordinates of each sample. We assume there are spatially-related
genes among all genes. We assume the genes fall into several profiles. The
genes in the same profile have the same expression over the sampled spatial
context. For each profile, we can summarize the gene expression for such group
by one representative factor. Consider there are $K$ profiles, we can then
form our model as:

$\displaystyle Y$
$\displaystyle=\sum_{k=1}^{K}\boldsymbol{f}_{k}\boldsymbol{l}_{k}^{T}+E$
$\displaystyle=FL^{T}+E$

where $\boldsymbol{f}_{k}\in\mathbb{R}^{n}$ is the factor gene in profile $k$
capturing the gene expression pattern in that group,
$\boldsymbol{l}_{k}\in\mathbb{R}^{p}$ indicates the loading coefficients of
each factor. We assume there is random noise $E\in\mathbb{R}^{n\times p}$ in
observations and $E_{ij}\overset{\text{iid}}{\sim}N(0,1/\tau)$. The goal is to
select spatially-related genes based on spatial information in $S$. We also
extract the gene factors which capture information across all genes based on
feature extraction. We aim to find a latent gene space that respects spatial
structure.

##  4 Methods

In this section, we build up pipelines and models to achieve this goal. We
adapt matrix decomposition methods to incorporate spatial information. We
first preprocess gene expression to make them amenable for wavelet filtering.
We then use the Daubechies D4 Wavelet Transform [Daubechies, 1992] as wavelet
filter.

###  4.1 Gene Expression Over Location

We leverage a pre-processing step, which we call _input generation_ , to
combine spatial and expression data sources. We have both gene expression
measurements and spatial coordinates for each sample. Each gene has an
expression over each sample, and we can draw the sample over a 2D map by their
coordinates.

Figure 2: The most dense gene has a structured expression pattern.

As shown in Figure 2, the gene with the least sparsity in expression shows
varying levels of expression across different spatial regions. Intuitively, we
may consider the spatial expression pattern from this gene to be spatially
related. However, most genes are sparse and do not show fluctuations in gene
expression. Therefore, we first filter out those genes using a kOverA filter:
we only keep genes that have an expression measure above $A$ in at least $k$
samples [Gentleman et al., 2021]. This has the effect of removing genes that
are rarely active, though it is possible that strongly expressed genes still
show no spatial relationships.

We want gene expression data to have a neat grid structure that can be
expressed as a matrix form. However, at first, expression measurements are
roughly staggered over the spatial locations, as shown in Figure 3. This is a
consequence of the sampling strategy adopted by the 10x genomics Visium
platform. To obtain an evenly resampled version of the expression pattern, we
divide the two-dimensional space into several partitions and compute the local
average of each partition. Detailed are shown in Algorithm 1. Each gene has an
expression over the grid and we now can use a matrix
$Z\in\mathbb{R}^{D^{2}\times P}$ to represent it. The matrix $Z$ becomes the
new input for analysis. We then obtain the same formulation as in Section 3,
where we take $n=D^{2}$ and each location corresponds to a new observation. We
have $\boldsymbol{f}_{k}\in\mathbb{R}^{D^{2}},k=1,\ldots,K$ as the factor
genes, computed by vectorizing the gene expression matrix.

Figure 3: Left: The gene expression over scattered samples over locations.
Right: The local average of gene expression over grid locations, used as input
to our algorithms.

The choice of $D$ depends on the wavelet scale used in the next section. We
choose $D$ to be dyadic to prevent handling edge cases in the wavelet
transformation (recalling Section 2.1). The dyadic size is also a natural
choice in image analysis. In particular, we choose $D=2^{J}$, where $J$ is the
level of scale in the wavelet transformation. Larger $J$ allows finer
recovery, but also requires more parameters and may result in overfitting.

###  4.2 Wavelet Transformation and Shrinkage

We apply a wavelet transformation with shrinkage to denoise the gene
expression matrix and smooth observed spatial expression patterns. In
simulations, we find that this technique gives more accurate recovery as the
signal-to-noise ratio decreases and yields less noisy visualizations for the
spatially-related genes.

We implement wavelet transformation and shrinkage on the processed data matrix
$Z\in\mathbb{R}^{D^{2}\times P}$. Each column of $Z$ is a vectorized
expression matrix. For each column, we first reshape it into a $D\times D$
matrix and apply the wavelet transformation. We have a coefficient list with
$W$ coefficients associated with each gene. Then we conduct wavelet shrinkage
on the coefficient list using the threshold strategies described in Section
2.1. Finally, we vectorize the coefficient list of each gene into a length $W$
vector for each gene. Stacking all vectors together, we have a coefficient
matrix $C\in\mathbb{R}^{W\times P}$. Note that the wavelet scale is specified
by the size $D$ of the input matrix, i.e., we have $D=2^{J}$, where $J$ is the
scale. The details are specified in Algorithm 2.

###  4.3 Matrix Decomposition

The matrix of the shrunk wavelet coefficients $C$ gives a summary of the
denoised matrix. This allows improved reconstruction of gene expression data.
To obtain a low-dimensional approximation, we apply matrix factorization on
the coefficient matrix $C$ after transformation. We use the same notation in
Section 3, with the subscript $c$ to denote the decomposition on coefficient
matrix $C$.

The resulting singular vectors can be used to estimate spatially structured
factor genes (Section 2.2). We use SVD as a frequentist approach to estimating
$\hat{F_{c}}$ and $\hat{L}_{c}$. We also conduct EBMF (Section 2.3), using the
posterior expectation of $F_{c}$ and $L_{c}$ to estimate $\hat{F_{c}}$ and
$\hat{L_{c}}$. EBMF can select the number of factors $K$ by itself, whereas
$K$ must be manually specified in the SVD. The choice of $K$ is informed by
inspecting the scree plot, as in spectral clustering and PCA. We choose the
number of factors when current factors explain sufficient information and the
diminishing returns of additional factors are no longer worth the additional
cost. Given $\hat{F}_{c}$ and $\hat{L}_{C}$, we can compute the estimated
coefficient matrix as $\hat{C}=\hat{F}_{c}\hat{L}_{c}^{T}$.

###  4.4 Inverse Wavelet Transformation

To transfer the coefficient matrix $\hat{C}$ back to the estimated location-
by-gene matrix $\hat{Z}$, we apply the inverse wavelet transformation on
columns of $\hat{C}$, each of which is a vectorized coefficient list for one
gene. This results in a $D^{2}$ expression matrix for each gene. Then we
vectorize the matrix and stack all vectors together. This yields the
reconstructed matrix $\hat{Z}$. Details are given in Algorithm 3. For
visualization, we also conduct a similar process for each column of
$\hat{F_{c}}$. In this case, each column is associated with an factor gene. By
applying the inverse wavelet transformation to each factor gene, we can build
spatial gene expression matrices $M_{1},\ldots,M_{K}$ representing gene
factors.

Algorithm 1 Input Generation

1:Sample by gene matrix $Y\in\mathbb{R}^{N\times P}$, spatial matrix
$S\in\mathbb{R}^{N\times 2}$, size of gene expression matrix $D$

2:Compute the range of $x,y$ coordinates from $S$, compute the coordinates of
vertices of big rectangle map $B$ cover all $N$ samples spatially.

3:Partition interval $x$ and interval $y$ into $D$ equal length interval,
together get $D^{2}$ partitions over rectangle $B$.

4:while $i\leq P$ do $\triangleright$ Consider gene $i$ expression over map
$B$

5: Select $i$-th column of gene matrix $Y$

6: Compute the local average of gene expression of gene $i$ in each partition

7: Get $D^{2}$ matrix $G_{i}$ as gene expression of gene $i$

8: Vectorize $G_{i}$ into a vector $g_{i}$ with length $D^{2}$

9:end while

10:Stacking all $\\{g_{i}\\}_{i=l}^{P}$ together into matrix
$Z\in\mathbb{R}^{D^{2}\times P}$

11:return matrix $Z$ $\triangleright$ A transformed gene expression matrix

Algorithm 2 Wavelet Transformation and Shrinkage

1:Location by gene matrix $Z\in\mathbb{R}^{D^{2}\times P}$, threshold method,
optional threshold parameter $\tau$. $\triangleright$ Apply threshold on
wavelet coefficient if threshold method is specified

2:while $i\leq P$ do $\triangleright$ Consider gene $i$ expression as matrix
$G_{i}$

3: Select $i^{\text{th}}$ $g_{i}$ column of gene matrix $Z$

4: Form $g_{i}$ into expression matrix $G_{i}$ with size $D^{2}$

5: Apply 2-D discrete wavelet transformation over $G_{i}$, get coefficient
list $C_{i}$, with number of coefficients $W$

6: Apply wavelet shrinkage over $C_{i}$ with optional parameter $\tau$

7: Vectorize $C_{i}$ into a long vector $c_{i}$ with length $W$

8:end while

9:Stacking all $\\{c_{i}\\}_{i=l}^{P}$ together into matrix
$C\in\mathbb{R}^{W\times P}$

10:return coefficient matrix $C$ $\triangleright$ Column $i$ of $C$ store the
coefficient of from gene $i$

Algorithm 3 Inverse Wavelet Transformation

1:(reconstructed) Coefficient matrix $C$

2:while $i\leq P$ do

3: Select $i^{\text{th}}$ $c_{i}$ column of gene matrix $C$

4: Form $c_{i}$ into coefficient list $C_{i}$ with number of coefficients $W$

5: Apply 2-D inverse wavelet transformation over $C_{i}$, get post-wavelet
expression matrix $\hat{G_{i}}$

6: Vectorize $\hat{G_{i}}$ into a long vector $\hat{g_{i}}$ with length
$D^{2}$

7:end while

8:Stacking all $\\{\hat{g_{i}}\\}_{i=1}^{P}$ together into matrix
$\hat{Z}\in\mathbb{R}^{D^{2}\times P}$

9:return matrix $\hat{Z}$ $\triangleright$ Post processing reconstructed gene
expression matrix

###  4.5 Evaluation

We evaluate our proposal both quantitatively and qualitatively. Our
quantitative results measures reconstruction error between the estimated
$\hat{Z}$ and $Z$ using the Frobenius norm $\|\hat{Z}-Z\|_{F}^{2}$.

This evaluation is also used for hyperparameter tuning, such as the scale $J$
of the wavelet and the choice of wavelet thresholding method. Our qualitative
result is given by visualization of the estimated factor gene expression
matrices, with emphasis on capturing global spatial structure.

We use cross-validation in computing reconstruction error and calculating
gene-wise errors. This evaluation is helpful in selecting spatially related
genes. We also found a simple connection between the gradient of gene
expression and spatial contribution, discussed in Section 6. This discussion
requires a measure of spatial expression smoothness. To this end, a simple
step computing the fluctuation of gene expression would give a spatial gene
selection that coincides with the reconstruction error selection. For
calculating successive differences, consider $Z$ as the input matrix and let
$\delta_{jk}=(z_{j,k+1}-z_{jk})^{2},k=1,\ldots,D^{2}$. Then we define the
gradient by $\sum_{jk}\delta_{jk}$.

##  5 Simulation

We setup simulations to see whether representing spatial structure with a
wavelet basis supports denoising and visualization of gene expression. We will
see improved recovery in the low signal-to-noise ratio regime. Qualitatively,
we also find factor gene visualizations to be more spatially consistent.

We first describe the simulation mechanism. We the set number of factors to
$K=9$. We generate $K$ gene expression matrices
$M_{1},\ldots,M_{K}\in\mathbb{R}^{D^{2}}$, each representing a factor gene.
These factor gene expression patterns are shown in Figure 4. We vectorize the
patterns $M_{k}$ into factor genes $\boldsymbol{f}_{k}$, which are then scaled
so $\|\boldsymbol{f}_{k}\|=1$. We obtain $F\in\mathbb{R}^{D^{2}\times K}$ by
stacking all the $\boldsymbol{f}_{k}$ together. We generate loadings
$\boldsymbol{l}_{1},\ldots,\boldsymbol{l}_{K}\in\mathbb{R}^{P\times K}$ by
drawing coordinates independently from $N(0,1)$. We set $D=32$ and wavelet
scale $J=5$. We similarly stack $\boldsymbol{l}_{k}$ to obtain
$L\in\mathbb{R}^{P\times K}$. We use $Z=FL^{T}$ as the ground-truth signal
matrix and add noise $E$ to yield data matrix $Z_{d}=FL^{T}+E$, where entries
in $E$ are zero-mean normal noise that corrupt the underlying signal. In our
analysis, we use a Daubechies D4 Wavelet as the wavelet filter, and for
coefficient shrinkage we use hybrid thresholding.

Figure 4: The gene expression pattern for nine factor genes in the simulation
experiment.

We use the pipeline from Section 4. The generated data $Z_{d}$ already
reflects the structure produced by the processing of Section 4.1. For
comparison, we also directly decompose the matrix $Z_{d}$ with the SVD without
any wavelet transformation. We call the resulting factors $\hat{F}_{raw}$ and
reconstruction $\hat{Z}_{raw}$. We also have the same quantitative
reconstruction error and qualitative visualization of factor gene
$\hat{F}_{raw}$. We also measure the size of the gradients across the
estimated gene expression matrix $\hat{M}_{i}$. The gradient is computed by
successive difference of neighboring image pixels. We calculated the sum of
the squares of the gradient. This measurement shows whether the gene
expression matrix has been smoothed. This property is of interest, since
smoother estimates are often more visually appealing.

We denote the signal-to-noise ratio as
$\text{SNR}=\frac{\text{sd}(Z)}{\text{sd}(E)}$, where $sd$ stands for standard
deviation. Let $r=\frac{1}{\text{SNR}}$. We specify 19 evenly spaced settings
of $r$ from 1 to 10. For each setting, we run 100 replicates with different
simulated $Z_{d}$ and apply wavelet and SVD-based dimensionality reduction.
The resulting average errors across $r$ are shown in the Figure 5. The wavelet
and shrinkage technique has better performance when $r$ is larger than 5,
i.e., the low signal-to-noise ratio regime. The gradient of the gene
expression matrix under two methods is shown in Figure 5. The wavelet and
shrinkage approach smooths edges in the factor gene expression image, giving a
more interpretable visualization and lower error in this low SNR setting.

Figure 5: The reconstruction error and gradient for estimated gene-by-location
matrix with different SNR. The $x$-axis shows the rate of
$\frac{1}{\text{SNR}}$ The left subplot shows the gradient changes. The
gradient of gene expression is always lower with the wavelet technique, a
consequence of its smoothing property. The right subplot shows the error – the
wavelet method has lower error as the magnitude of noise increase.

The factor genes are shown in Figure 6. The SVD without wavelet technique has
erratic, outlying pixels, especially in the last three genes. The
visualization is sensitive to outliers. In contrast, the SVD combined with
wavelet has smoother patterns. Like the SVD-based method, it appears to have
mixed several of the true underlying factors in each of the recovered ones.
Moreover, the sharp boundaries visible in the SVD factors become smoothed over
in the wavelet-decomposition. The wavelet method applies a decomposition on
coefficients space after thresholding, while SVD operates on individual
pixels. The SVD capture more information, but also emphasize nuisance
information induced by errors. Wavelet method also reduces model complexity,
improving estimation accuracy. Other non-parametric methods, such as Fourier
transformation and Gaussian Process-based methods also operate on coefficients
space. Still, they would struggle to capture sharp transitions, since their
bases are smooth functions.

Nonetheless, both the SVD and wavelet-based visualizations reflect spatial
trends in the true factor genes (in Figure 4).

(a) Factor genes without wavelet smooth

(b) Factor genes with wavelet smooth

Figure 6: Factor gene visualization. Figure (a) implements SVD on $Z_{d}$, the
visualization of column of $\hat{F}_{raw}$. Figure (b) implements SVD on
coefficient matrix, the visualization of $\hat{F}$.

##  6 Real Data Experiment

In this section, we show that wavelet and shrinkage technique reduces
reconstruction error quantitatively. We ran our method on a public spatially
resolved transcriptomics data [Weber, 2021]. The dataset represents a single
biological sample from the human brain dorsolateral prefrontal cortex (DLPFC)
region, measured with the 10x Genomics Visium platform. Further, we identify a
simple connection between the gradient of gene expression and quantitative
error. A simple step computing the fluctuation of gene expression alone
(calculating successive difference of the gene expression image) selects genes
that have reduced reconstruction error when using wavelet-guided
dimensionality reduction.

We first process the ST data through our pipeline. The ST data contains 4992
observations with 33538 genes. Expressions from most genes are sparse. We
implement the pipeline from Section 4. We pre-process as in Section 4.1 and
then apply a kOverA filter to select genes. We find $k=3,A=7$ gives us 721
genes with average expression as 2.71. Then we apply Algorithm 1 to transfer
the sample-by-gene matrix to a grid-by-gene matrix. We set $D=64$ and wavelet
scale $J=6$. We obtain an image heatmap for each gene like in Figure 7.

Figure 7: Gene expression over the grid.

We then vectorize $P=721$ image expression and stacking vectors together, we
have generated input data $Z\in\mathbb{R}^{64^{2}\times 721}$. We run our
pipelines with and without wavelet transformation for evaluation. We first
implement SVD and EBMF on $Z$ or on coefficient matrix $C$ after applying
wavelet transformation and thresholding each column. We choose the number of
factors $K$ by examining singular values. If the singular value is larger than
500, we keep it as a factor. In EBMF, we set upper bound that $K$ to be
smaller than the corresponding $K$ in SVD, then let algorithm choose $K$
itself.

Quantitative evaluation metric is conducted by cross-validation on non-zero
entries. We set 5 folds cross-validation with replacement. In particular, we
select random $\frac{1}{5}$ non-zero entries in matrix $Z$ and set them to
zero, then we save matrix as the masked matrix (train data) $Z_{train}$. We
store the values and position of the masked entries as test data $Z_{test}$.
We then ran methods on $Z_{train}$. The result shows a difference in whether
to use the wavelet technique. The result from SVD and EBMF are close to each
other coordinate-wise, hence we only show one of them in some results below.We
include the comparison between SVD and EBMF in the Appendix A. As in Section
5, we have $\hat{F}$ and $\hat{Z}$ from the wavelet approach and
$\hat{F}_{raw}$ and $\hat{Z}_{raw}$ on the original data.

###  6.1 Total Error and Parameter Tuning

We compute the reconstruction loss,

$\frac{1}{N}\sum_{i,j\in\text{test}}(\hat{Z}_{ij}-Z_{ij})^{2}$

where $N$ is the number of test entries. We evaluate the loss only on test
entries. We first tune parameters. We set up three settings: decompose the raw
data with SVD or EBMF; wavelet transformation with hybrid thresholding;
wavelet transformation with manual thresholding with threshold
$\tau=10,20,\ldots,100$. In each setting, we ran 100 replicates. The
reconstruction error shown in Figure 8.

Figure 8: The reconstruction error for different wavelet parameters. The upper
plot contains the result from all settings, the leftmost $\tau=-1$ result is
the decomposition without wavelet transformation. $\tau=0$ indicates the
hybrid thresholding. The bottom plot zooms into the experiment under manual
threshold for $\tau\in\left[10,100\right]$.

As shown in the upper panel of Figure 8, wavelet thresholding with manually
set $\tau$ reduces error compared to decomposing raw data. The wavelet has the
most positive effect when $\tau=40$, as shown in the bottom panel. We use
$\tau=40$ in our following analysis.

###  6.2 Genewise Error

We then evaluate how method performance varies for each gene by calculating
the genewise reconstruction error. This reveals genes whose strong spatial
expression structure leads to improved performance when using a wavelet basis.
We still hold out $\frac{1}{5}$ of the entries at random as a test set. We
compute the reconstruction error of test entries on each column of
$\hat{Z}-Z$. We calculate entry-wise loss across 100 replicates and estimate
the average loss. We compare genewise errors with or without wavelet
transformation. The SVD and EBMF show the same result regarding the
decomposition method. We show the result of EBMF in Appendix A Figure 14.

Figure 9: The reconstruction error for each gene. The $x$-axis is the error
from EBMF on raw data, the $y$-axis is the error from the combination of
wavelet thresholding and EBMF.

As we can see, most of the genes have lower reconstruction error when directly
applying the SVD or EBMF. However, for some genes, wavelet smoothing reduces
error. For example, this is seen in genes 712, 713, 716, 719, 715, 710, 711,
721, 717, 596, 373, and 289. We conjecture that these genes have acute
fluctuations as well as clear spatial patterns. The expression matrix of these
genes would have a larger gradient and distinct edges. To verify our
conjecture, we calculated the sum of squares of the gradient of each gene
expression matrix. We show the result in Figure 10.

Figure 10: The sum of squares of the gradient for each gene. The $y$-axis is
the index of each gene. Each bar’s color shows whether that gene has better
performance when using a wavelet basis. The M on the $x$-axis stands for
"millions".

The result verifies our conjecture: genes with larger gradients have better
reconstruction under the wavelet-guided decomposition. This suggests a pre-
processing step for selecting wavelet-suited genes by calculating the gradient
of each gene. The spatially related genes would have a lower reconstruction
error. Among these spatially related genes, one possibility is that we can
divide them into two groups, one for decomposition directly and the other for
the wavelet technique.

###  6.3 Factor Genes

Now we visualize the top factor genes and genes with high loadings on these
factor genes. In Figure 11(a), we decompose $Z$ using the SVD and plot the
first factor gene and matrix slides of genes with the largest 5 loadings on
that factor. The factor gene captured the same patterns as the original genes.
To improve visualization, we find the analogous wavelet-based factors (we use
manual thresholding with $\tau=40$), shown in Figure 11(b).

(a) Factor genes without wavelet smooth

(b) Factor genes with wavelet smooth

Figure 11: The top left figure is first factor gene, the following figure by
row is genes with largest loadings on first factor gene: gene 712, 713, 716,
719, 715.

The wavelet thresholding approach smooths over edges in the original
visualization. The first factor gene captures the global spatial expression.
We have a similar result for the second factor gene, as shown in Figure 12.

(a) Factor genes without wavelet smooth

(b) Factor genes with wavelet smooth

Figure 12: The top left figure is the second factor gene, the following figure
by row is genes with largest loadings on second factor gene: gene 596, 289,
373, 712, 578.

The second factor gene is orthogonal to the first one, capturing different
spatial structures in gene expression. Different factor genes capture global
and local structure, and using a wavelet decomposition provides denoised
spatial expression visualization.

In conclusion, we can use this dimensionality reduction technique for spatial
gene selection and extraction. We select wavelet-suited genes based on the
calculated gradient. Then, we can select spatially related genes based on
reconstruction error via cross-validation. Alternatively, we can extract
factor gene to capture spatial information and use factor genes for
visualization and further analysis.

###  6.4 Code Vignette

We provide a small code block to show vignette of generating Figure 11 and 12,
shown in code block 1.

⬇

1library(waveST)

2res = kOverA_ST(k = 3,A = 7)

3waveST = waveST(data = res$df, spatial = res$viz)

4

5#### Figure 11

6waveST = decompose(waveST, wavemethod = "raw", decom_method = "SVD", K = 5)

7plot(waveST, k = 1, wave = FALSE)

8

9#### Figure 12

10waveST = decompose(waveST, wavemethod = "wave", decom_method = "SVD", K = 5,
tau = 40)

11plot(waveST, k = 1, wave = TRUE)

Listing 1: The basic wavelet-based dimensionality reduction workflow.
kOverA_ST transforms the original spatial expression context into an image.
This is processed by waveST. The final dimensionality reduction step is
performed by decompose.

We have made a new package waveST containing the workflow we developed. The
package is available at https://github.com/OliverXUZY/waveST. The kOverA_ST
function reduces data from an original [Weber, 2021] class input using the
kOverA technique. Then we use waveST function to construct a S4 class waveST,
containing input generated by Algorithm 1. This object-oriented approach
stores properties of the original spatial experiment and simplifies downstream
calls, like decomposition and visualization. We use the decompose function to
apply all decomposition methods. In line 6, we apply the SVD to our original
data, setting the number of factors to 5. In line 11, we apply the wavelet-
based reduction and apply a manual threshold with $\tau=40$. We use plot and
k=1 to visualize the first factor gene and matrix slides of the genes with the
largest 5 loadings on that factor.

##  7 Conclusion

We have proposed a pipeline for dimensionality reduction that respects spatial
structure. Both simulations and real data experiments demonstrate that wavelet
and shrinkage techniques show positive results in spatially resolved
transcriptomics data. We highlight the idea of combining image processing
techniques and statistical methods for application in a spatial genomics
context. One future direction is splitting genes into a groups suited or not
for wavelet-based decomposition, and implementing decomposition with or
without wavelet. Another direction is to focus input generation on only those
genes that are thought to be spatially related. For genes not related to
spatial information, we may perform regular decomposition on original data
$Y\in\mathbb{R}^{N\times P}$, abandoning the spatial information of those
genes. We expect this to improve reconstruction performance. In further
analysis, it is worth considering other wavelet smoothing techniques and
wavelet filters. The current methods incorporate little biological
information. Bringing more domain knowledge will require further techniques,
but is expected to yield better results. The input generation computes local
average over even grids; however, it is possible to apply the wavelet method
for irregularly spaced data [Nason, 2008]. We hope wavelet methods will be
useful in adapting existing methods for statistical genomics to the spatial
setting.

## Acknowledgments

We thank Joseph Arthur for valuable discussions. This research was performed
using the compute resources and assistance of the UW-Madison Center For High
Throughput Computing (CHTC) in the Department of Computer Sciences. The CHTC
is supported by UW-Madison, the Advanced Computing Initiative, the Wisconsin
Alumni Research Foundation, the Wisconsin Institutes for Discovery, and the
National Science Foundation, and is an active member of the OSG Consortium,
which is supported by the National Science Foundation and the U.S. Department
of Energy’s Office of Science.

## References

  * Abu-Jamous and Kelly [2018] B. Abu-Jamous and S. Kelly.  Clust: automatic extraction of optimal co-expressed gene clusters from gene expression data.  _Genome biology_ , 19(1):1–11, 2018.
  * Alter et al. [2000] O. Alter, P. O. Brown, and D. Botstein.  Singular value decomposition for genome-wide expression data processing and modeling.  _Proceedings of the National Academy of Sciences_ , 97(18):10101–10106, 2000.
  * Asp et al. [2019] M. Asp, S. Giacomello, L. Larsson, C. Wu, D. Fürth, X. Qian, E. Wärdell, J. Custodio, J. Reimegård, F. Salmén, et al.  A spatiotemporal organ-wide gene expression and cell atlas of the developing human heart.  _Cell_ , 179(7):1647–1660, 2019.
  * Berglund et al. [2018] E. Berglund, J. Maaskola, N. Schultz, S. Friedrich, M. Marklund, J. Bergenstråhle, F. Tarish, A. Tanoglidi, S. Vickovic, L. Larsson, et al.  Spatial maps of prostate cancer transcriptomes reveal an unexplored landscape of heterogeneity.  _Nature communications_ , 9(1):1–13, 2018.
  * Bishop [1999] C. M. Bishop.  Variational principal components.  1999\.
  * Chen et al. [2015] K. H. Chen, A. N. Boettiger, J. R. Moffitt, S. Wang, and X. Zhuang.  Spatially resolved, highly multiplexed rna profiling in single cells.  _Science_ , 348(6233):aaa6090, 2015.
  * Close et al. [2021] J. L. Close, B. R. Long, and H. Zeng.  Spatially resolved transcriptomics in neuroscience.  _Nature Methods_ , 18(1):23–25, 2021.
  * Cochran et al. [1967] W. T. Cochran, J. W. Cooley, D. L. Favin, H. D. Helms, R. A. Kaenel, W. W. Lang, G. C. Maling, D. E. Nelson, C. M. Rader, and P. D. Welch.  What is the fast fourier transform?  _Proceedings of the IEEE_ , 55(10):1664–1674, 1967.
  * Daubechies [1992] I. Daubechies.  _Ten lectures on wavelets_.  SIAM, 1992.
  * Donoho and Johnstone [1995] D. L. Donoho and I. M. Johnstone.  Adapting to unknown smoothness via wavelet shrinkage.  _Journal of the american statistical association_ , 90(432):1200–1224, 1995.
  * Donoho and Johnstone [1994] D. L. Donoho and J. M. Johnstone.  Ideal spatial adaptation by wavelet shrinkage.  _biometrika_ , 81(3):425–455, 1994.
  * Donoho et al. [1994] D. L. Donoho, I. M. Johnstone, et al.  Ideal denoising in an orthonormal basis chosen from a library of bases.  _Comptes rendus de l’Académie des sciences. Série I, Mathématique_ , 319(12):1317–1322, 1994.
  * Dumitrascu et al. [2021] B. Dumitrascu, S. Villar, D. G. Mixon, and B. E. Engelhardt.  Optimal marker gene selection for cell type discrimination in single cell analyses.  _Nature communications_ , 12(1):1–8, 2021.
  * Gentleman et al. [2021] R. Gentleman, V. J. Carey, W. Huber, and F. Hahne.  _genefilter: genefilter: methods for filtering genes from high-throughput experiments_ , 2021.  R package version 1.74.0.
  * Kiselev et al. [2017] V. Y. Kiselev, K. Kirschner, M. T. Schaub, T. Andrews, A. Yiu, T. Chandra, K. N. Natarajan, W. Reik, M. Barahona, A. R. Green, and M. Hemberg.  Sc3 - consensus clustering of single-cell rna-seq data.  _Nature Methods_ , 2017.  URL http://dx.doi.org/10.1038/nmeth.4236.
  * Kiselev et al. [2018] V. Y. Kiselev, A. Yiu, and M. Hemberg.  scmap: projection of single-cell rna-seq data across data sets.  _Nature methods_ , 15(5):359–362, 2018.
  * Kuppe et al. [2020] C. Kuppe, R. O. R. Flores, Z. Li, M. Hannani, J. Tanevski, M. Halder, M. Cheng, S. Ziegler, X. Zhang, F. Preisker, et al.  Spatial multi-omic map of human myocardial infarction.  _BioRxiv_ , 2020.
  * Mallat [1989] S. G. Mallat.  A theory for multiresolution signal decomposition: the wavelet representation.  _IEEE transactions on pattern analysis and machine intelligence_ , 11(7):674–693, 1989.
  * Moffitt et al. [2016] J. R. Moffitt, J. Hao, G. Wang, K. H. Chen, H. P. Babcock, and X. Zhuang.  High-throughput single-cell gene-expression profiling with multiplexed error-robust fluorescence in situ hybridization.  _Proceedings of the National Academy of Sciences_ , 113(39):11046–11051, 2016.
  * Moncada et al. [2020] R. Moncada, D. Barkley, F. Wagner, M. Chiodin, J. C. Devlin, M. Baron, C. H. Hajdu, D. M. Simeone, and I. Yanai.  Integrating microarray-based spatial transcriptomics and single-cell rna-seq reveals tissue architecture in pancreatic ductal adenocarcinomas.  _Nature biotechnology_ , 38(3):333–342, 2020\.
  * Nason [2008] G. P. Nason.  _Wavelet methods in statistics with R_.  Springer, 2008.
  * Perperoglou et al. [2019] A. Perperoglou, W. Sauerbrei, M. Abrahamowicz, and M. Schmid.  A review of spline function procedures in r.  _BMC medical research methodology_ , 19(1):1–16, 2019.
  * Rödelsperger et al. [2021] C. Rödelsperger, A. Ebbing, D. R. Sharma, M. Okumura, R. J. Sommer, and H. C. Korswagen.  Spatial transcriptomics of nematodes identifies sperm cells as a source of genomic novelty and rapid evolution.  _Molecular biology and evolution_ , 38(1):229–243, 2021.
  * Rodriques et al. [2019] S. G. Rodriques, R. R. Stickels, A. Goeva, C. A. Martin, E. Murray, C. R. Vanderburg, J. Welch, L. M. Chen, F. Chen, and E. Z. Macosko.  Slide-seq: A scalable technology for measuring genome-wide expression at high spatial resolution.  _Science_ , 363(6434):1463–1467, 2019.
  * Shah et al. [2016] S. Shah, E. Lubeck, W. Zhou, and L. Cai.  In situ transcription profiling of single cells reveals spatial organization of cells in the mouse hippocampus.  _Neuron_ , 92(2):342–357, 2016.
  * Shang and Zhou [2022] L. Shang and X. Zhou.  Spatially aware dimension reduction for spatial transcriptomics.  _bioRxiv_ , 2022.
  * Srivatsan et al. [2021] S. R. Srivatsan, M. C. Regier, E. Barkan, J. M. Franks, J. S. Packer, P. Grosjean, M. Duran, S. Saxton, J. J. Ladd, M. Spielmann, et al.  Embryo-scale, single-cell spatial transcriptomics.  _Science_ , 373(6550):111–117, 2021.
  * Ståhl et al. [2016] P. L. Ståhl, F. Salmén, S. Vickovic, A. Lundmark, J. F. Navarro, J. Magnusson, S. Giacomello, M. Asp, J. O. Westholm, M. Huss, et al.  Visualization and analysis of gene expression in tissue sections by spatial transcriptomics.  _Science_ , 353(6294):78–82, 2016.
  * Svensson et al. [2018] V. Svensson, S. A. Teichmann, and O. Stegle.  Spatialde: identification of spatially variable genes.  _Nature methods_ , 15(5):343–346, 2018.
  * Thrane et al. [2018] K. Thrane, H. Eriksson, J. Maaskola, J. Hansson, and J. Lundeberg.  Spatially resolved transcriptomics enables dissection of genetic heterogeneity in stage iii cutaneous malignant melanoma.  _Cancer research_ , 78(20):5970–5979, 2018.
  * Townes and Engelhardt [2021] F. W. Townes and B. E. Engelhardt.  Nonnegative spatial factorization.  _arXiv preprint arXiv:2110.06122_ , 2021.
  * Velten et al. [2022] B. Velten, J. M. Braunger, R. Argelaguet, D. Arnol, J. Wirbel, D. Bredikhin, G. Zeller, and O. Stegle.  Identifying temporal and spatial patterns of variation from multimodal data using mefisto.  _Nature methods_ , pages 1–8, 2022.
  * Wang and Stephens [2021] W. Wang and M. Stephens.  Empirical bayes matrix factorization.  _Journal of Machine Learning Research_ , 22(120):1–40, 2021.
  * Weber [2021] L. M. Weber.  _STexampleData: Collection of spatially resolved transcriptomics datasets in SpatialExperiment Bioconductor format_ , 2021.  URL https://github.com/lmweber/STexampleData.  R package version 1.0.8.
  * Xia et al. [2019] C. Xia, J. Fan, G. Emanuel, J. Hao, and X. Zhuang.  Spatial transcriptome profiling by merfish reveals subcellular rna compartmentalization and cell cycle-dependent gene expression.  _Proceedings of the National Academy of Sciences_ , 116(39):19490–19499, 2019.
  * Zhu and Sabatti [2020] J. Zhu and C. Sabatti.  Integrative spatial single-cell analysis with graph-based feature learning.  _bioRxiv_ , 2020.

##  Appendix A Comparison Between SVD and EBMF

This section we showed the comparison between the reconstruction result
between SVD and EBMF.

###  A.1 Simulation

In this section, we show the element-wise comparison between the
reconstruction of SVD and of EBMF with wavelet technique.

Figure 13: The element-wise comparison of two reconstruction matrices. The
$x$-axis shows the result of SVD decomposition, the $y$-axis shows the result
of EBMF decomposition.

###  A.2 Real Data Result

In this section, we show the reconstuction error per gene result for SVD,
similar to Figure 14.

Figure 14: The reconstruction error for each gene. The $x$-axis is the error
from SVD on raw data, the $y$-axis is the error from the combination of
wavelet thresholding and SVD.

Then we show the result between SVD and EBMF with and without wavelet
technique in Figure 15.

(a) Error per gene without wavelet transformation.

(b) Error per gene with wavelet transformation.

Figure 15: The reconstruction error for each gene.

As we can see, SVD and EBMD have the very similar result.


# Activity Induced Enhanced Diffusion of a Polymer in Poor Solvent

Suman Majumder<EMAIL_ADDRESS>Amity Institute of Applied Sciences, Amity
University Uttar Pradesh, Noida 201313, India  Subhajit Paul
<EMAIL_ADDRESS>International Center for Theoretical Sciences, Tata
Institute of Fundamental Research, Bangalore-560089, India  Wolfhard Janke
<EMAIL_ADDRESS>Institut für Theoretische Physik,
Universität Leipzig, IPF 231101, 04081 Leipzig, Germany

###### Abstract

By means of Brownian dynamics simulations we study the steady-state dynamic
properties of a flexible active polymer in a poor solvent condition. Our
results show that the effective diffusion constant of the polymer $D_{\rm
eff}$ gets significantly enhanced as activity increases, much like in active
particles. The simulation data are in agreement with a theoretically
constructed Rouse model of active polymer, demonstrating that irrespective of
the strength of activity, the long-time dynamics of the polymer chain is
characterized by a universal Rouse-like scaling $D_{\rm eff}\sim N^{-1}$,
where $N$ is the chain length.

Biomolecules are subjected to athermal fluctuations originating from chemical
reactions or other energy conversions, rendering them fall out of equilibrium.
Often that is the underlying cause for a range of biological activities, e.g.,
bacterial motion, shape fluctuations of red-blood cell membranes, enzyme
catalysis Wu and Libchaber (2000); Park  _et al._ (2010); Ben-Isaac  _et al._
(2011); Jee  _et al._ (2018a, b). Hence, given the enormous progress in
understanding of active particles Ramaswamy (2010); Shaebani  _et al._ (2020);
Elgeti  _et al._ (2015), over the years a number of studies have emerged
investigating active polymers as well Winkler and Gompper (2020); Kaiser  _et
al._ (2015); Isele-Holder  _et al._ (2015, 2016); Bianco  _et al._ (2018);
Locatelli  _et al._ (2021); Das  _et al._ (2021); Paul  _et al._ (2021, 2022a,
2022b, 2022c). The motivation for studying active polymers stems from the need
of introducing variety in shape, flexibility and coupling topology in active
entities Paxton  _et al._ (2004); Sanchez  _et al._ (2011); Ebbens  _et al._
(2012); Buttinoni  _et al._ (2013); Kümmel  _et al._ (2013); Zöttl and Stark
(2014). Besides, it is intriguing to check how the knowledge of polymer
physics De Gennes (1979); Doi (1996); Rubinstein and Colby (2003) can be
deployed to understand active matters.

To date active polymer models can be classified into two categories. The
straightforward way is to consider the monomers as active particles and then
connecting them linearly with a bond constraint Kaiser  _et al._ (2015);
Bianco  _et al._ (2018); Das  _et al._ (2021); Paul  _et al._ (2021, 2022a,
2022b, 2022c). In the other approach one takes a passive polymer, i.e.,
without any activity, immersed in a bath of active particles Harder  _et al._
(2014); Kaiser and Löwen (2014); Chaki and Chakrabarti (2019); Liu  _et al._
(2019); Anderson  _et al._ (2022); Goswami  _et al._ (2022). Apart from being
motivated by biological entities, current advanced techniques allows one to
synthesize polymers made of active colloids connected artificially by DNA or
freely jointed droplets Biswas  _et al._ (2017); McMullen  _et al._ (2018).
Theoretically, the constituent monomers can be made active by (i) introducing
a local force tangential to the polymer backbone Bianco  _et al._ (2018), (ii)
considering the monomers having Brownian activity Das  _et al._ (2021); Paul
_et al._ (2022c), or (iii) Vicsek-like activity Paul  _et al._ (2021, 2022a,
2022b).

The conformational and dynamic properties of a passive polymer is
characterized by various well established scaling laws De Gennes (1979); Doi
(1996); Rubinstein and Colby (2003). A polymer undergoes a coil to globule
transition upon changing the solvent condition from good (where monomer-
solvent interaction dominates) to poor (where monomer-monomer interaction
dominates). The spatial extension of the conformations measured in terms of
the radius of gyration $R_{g}$ typically follows the scaling $R_{g}\sim
N^{\nu}$ with respect to the degree of polymerization or number of monomers
$N$. The value of the exponent $\nu\approx 3/5$ and $1/3$ characterizes the
conformations, respectively, in good and bad solvents. The dynamics of a
polymer under a good solvent condition in the free-draining limit, i.e.,
ignoring hydrodynamics, is characterized by Rouse scaling $D\sim N^{-1}$,
where $D$ is the diffusion coefficient of the center of mass of the polymer
Rouse Jr (1953). In presence of hydrodynamics, one expects in the Zimm model
Zimm (1956) (with excluded volume) a scaling $D\sim N^{-\nu}$. In a poor
solvent condition, however, there is no consensus among the available studies
Naghizadeh and Kovac (1987); Milchev  _et al._ (1993). Anomalous behavior with
resemblance to slow glassy dynamics has also been reported Milchev and Binder
(1994).

In case of active polymers too, the focus has been on understanding the
nonequilibrium steady-state conformations and dynamics. In particular,
attempts have been made to adapt scaling theories of passive polymers in good
solvent to study active polymers under the same condition. Recently, Bianco et
al. Bianco  _et al._ (2018) observed an activity induced collapse of a polymer
in good solvent, reminiscent of motility induced phase separation of active
particles Fily and Marchetti (2012). In contrast, polymers made of active
Brownian monomers do not exhibit such collapse in good solvent. On the
dynamics front, in steady state an enhancement of the diffusion coefficient
has been observed for all cases. In all these studies the self-avoidance in
good solvent condition has been mimicked by considering a purely repulsive
interaction among the constituent monomers. Only recently, an interaction
potential with both attractive and repulsive components has been considered
Das  _et al._ (2021); Paul  _et al._ (2021, 2022a, 2022b, 2022c). A passive
polymer having such an interaction exhibits a temperature dependent coil-
globule transition. We have shown that a polymer with active Brownian monomers
in bad solvent condition, exhibits a transition from a globular state at small
activity to coil-like conformations at large activity Paul  _et al._ (2022c).
Hence, it is expected that the dynamics of such active polymers in poor
solvent would also reveal interesting features.

In this Letter, by means of computer simulation supported by analytical
reasoning we investigate the steady-state dynamics of a flexible coarse-
grained model polymer consisting of active Brownian monomers in a poor solvent
condition. To probe the dynamics we have monitored the motions of the center
of mass of the polymer and two different tagged monomers, viz., the central
and end monomers. All these motions exhibit long-time diffusive behaviors
allowing us to calculate the diffusion coefficient $D$ of the polymer. We show
that a universal Rouse-like scaling of $D$ with respect to $N$ is maintained
at all considered strengths of activity, albeit, the polymer could be in a
globular or coiled state.

We consider a flexible polymer consisting of monomers of diameter $\sigma$ at
positions $\vec{r}_{i}$ connected linearly via spring-like bonds. Its dynamics
is governed by the over-damped Langevin equations

$\begin{split}\partial_{t}{\vec{r}}_{i}&=\frac{D_{\rm{tr}}}{k_{B}T}[f_{p}\hat{n}_{i}-\vec{\nabla}U_{i}]+\sqrt{2D_{\rm{tr}}}\,\vec{\Lambda}_{i}^{\rm{tr}},\\\
\partial_{t}{\hat{n}}_{i}&=\sqrt{2D_{\rm{rot}}}(\hat{n}_{i}\times\vec{\Lambda}_{i}^{\rm{rot}}),\end{split}$
(1)

where $f_{p}$ is the strength of the self-propulsion force acting along the
unit vector $\hat{n}_{i}$ that changes stochastically with time, and
$U_{i}=V_{\rm{FENE}}+V_{\rm{LJ}};$ where
$V_{\rm{FENE}}(r)=-0.5KR^{2}{\rm{ln}}[1-((r-r_{0})/R)^{2}]$ takes care of the
bonds between successive monomers and
$V_{\rm{LJ}}(r)=4\epsilon[(\sigma/r)^{12}-(\sigma/r)^{6}]$ is the non-bonded
interaction accounting for the self-avoidance, with a strength $\epsilon$ . In
Eq. (1) $D_{\rm{tr}}$ and $D_{\rm{rot}}$ are, respectively, the translational
and rotational diffusion constants. The components of
$\vec{\Lambda}_{i}^{\rm{tr}}$ and $\vec{\Lambda}_{i}^{\rm{rot}}$ are
independent white Gaussian noises with zero-mean and unit-variance. In our
simulations we have fixed their relative importance
$D_{\rm{tr}}/D_{\rm{rot}}\sigma^{2}\equiv\Delta=1/3$. $D_{\rm{tr}}$ is related
to the temperature $T$ and the drag or friction coefficient $\gamma$ as
$D_{\rm{tr}}=k_{B}T/\gamma$. For convenience we have chosen $\gamma=1$. Time
is measured in units of $\tau_{0}=\sigma^{2}\gamma/\epsilon$ ($\propto
1/D_{\rm{rot}}=\Delta\sigma^{2}\gamma/k_{B}T$ at fixed $k_{B}T/\epsilon$) and
we have used an integration time step of $10^{-5}\tau_{0}$. Here onward, the
activity strength $f_{p}$ will be expressed in terms of the dimensionless
Péclet number $Pe={f_{p}\sigma}/{k_{B}T}$. The case of $Pe=0$ corresponds to a
passive polymer. We perform our simulations for different $Pe$ at
$T=0.1\epsilon/k_{B}$, well below the coil-globule transition temperature of
the passive polymer Majumder  _et al._ (2017), thus mimicking a poor solvent
condition. For further technical details on the model and simulation method we
refer to the Supplemental Material 111See Supplemental Material at […] which
includes Refs. Rouse Jr (1953); Howse  _et al._ (2007). The Supplemental
Material contains a detailed description of the model and method of simulation
used, and a few representative steady-state trajectories of the polymer. It
also includes details of the analogous Rouse model for active polymer in poor
solvent, and a derivation of the corresponding diffusion constant of the
center of mass. Furthermore, we present details of the performed fitting
exercise showing the agreement between theory and simulation data. .

We started our simulations using self-avoiding coils as initial condition.
Then, we allow the system to reach its steady state by running it for
sufficiently long times. Typical steady-state conformations of a polymer of
length $N=128$ obtained for different $Pe$ are presented in Fig. 1. The
conformation of the passive polymer, i.e., $Pe=0$, is a perfectly collapsed
globule. It remains in such a globular state for relatively smaller activities
$Pe\leq 25$ as well. For intermediate values of $Pe$ one may observe a globule
or head-tail-like conformations (as the one presented for $Pe=37.5$). For even
larger $Pe$ the polymer becomes an extended coil. The corresponding
quantitative picture is presented in Fig. 1 in the form of the distribution of
end-to-end distance $R_{e}=|\vec{r}_{1}-\vec{r}_{N}|$ where $\vec{r}_{i}$ is
the position of the $i$-th monomer. It shows that for $Pe<37.5$ the peaks are
at $R_{e}\approx 3$, indicative of collapsed globules. The decrease in peak
height as $Pe$ increases is reflective of the fact that the probability of
getting a collapsed globule decreases and encountering a head-tail-like
conformations increases. For $Pe>37.5$ the distribution broadens and the peak
position shifts towards $R_{e}>15$ suggesting a dominance of coil-like
conformations. The overall picture is reminiscent of the temperature driven
coil-globule transition of a passive polymer. Here, it is driven by the
activity strength. In fact, we have confirmed in Ref. Paul  _et al._ (2022c)
that the conformations obey the scaling law $R_{g}\sim N^{\nu}$ with $\nu=1/3$
and $\approx 3/5$, respectively, at small and large $Pe$.

Figure 1:  Typical steady-state conformations of a polymer of length $N=128$
at different activity strength $Pe$, obtained from simulations at a fixed
temperature $T=0.1\epsilon/k_{B}$. The plots are for corresponding normalized
distributions of the end-to-end distance $R_{e}$.
Figure 2:  Steady-state mean square displacement ${\rm MSD}_{i}(t)$ of the (a)
center of mass, (b) central bead, and (c) end beads of a polymer of length
$N=128$ for different activity strengths $Pe$. Plots in (d), (e), and (f)
present a comparison among the different MSDs for different values of $Pe$.
The dashed lines represent different power laws. All data are obtained from
simulations at temperature $T=0.1\epsilon/k_{B}$.

To probe the dynamics we monitor the motion of the center of mass (cm),
central monomer or bead (cb) and end beads (eb) of the polymer. A bare look at
the typical trajectories over a fixed time period reveals that although
motions are random in general, the distance covered varies significantly for
different $Pe$, clearly suggesting a difference in dynamics (see Fig. S1 in
the Supplemental Material Note (1)). To quantify the differences, from the
obtained trajectories, we calculate the corresponding mean square
displacements

${\rm
MSD}_{i}(t)=\langle\left[\vec{r}_{i}(t)-\vec{r}_{i}(0)\right]^{2}\rangle;~{}i\equiv{\rm
cm,cb,and~{}eb},$ (2)

as a function of time $t$. Figure 2(a) shows that the cm exhibits a typical
long-time diffusive motion $\sim t$, with pronounced short-time ballistic
behavior as the activity increases. A similar long-time behavior is also
observed for the motions of cb and eb, shown respectively, in Figs. 2(b) and
(c). Significantly, different is the appearance of an intermediate regime
which becomes longer as the activity strength $Pe$ increases. In this regime
the behavior of the central bead appears to be $\sim t^{2/3}$ for large $Pe$,
which may lure one to consider it as hydrodynamic Zimm’s scaling of a passive
polymer in good solvent Paul  _et al._ (1991); Dünweg and Kremer (1991, 1993);
Hinczewski  _et al._ (2009). However, this is very unlikely and probably a
mere coincidence since our simulations do not preserve hydrodynamics. The end
beads show an extended intermediate regime, although, the corresponding power-
law exponent seems to be smaller than $2/3$. This rather suggests a Rouse-like
behavior, expected for a passive polymer in good solvent.

Figure 3:  Plots of the time-dependent exponent $\alpha_{i}(t)$ for the data
presented in Fig. 2. The grey shades are introduced to distinguish the early-
time regime from the long-time diffusive regime.

For a better understanding of the time-dependent power-law behavior ${\rm
MSD}_{i}\sim t^{\alpha}$, we calculate the instantaneous exponent

$\alpha_{i}(t)=\frac{d\ln{\rm MSD}_{i}(t)}{d\ln t};~{}i\equiv{\rm
cm,cb,and~{}eb}.$ (3)

Corresponding plots of $\alpha_{i}(t)$ vs. $t$ are presented in Fig. 3 for all
considered values of $Pe$. The exponent $\alpha_{\rm cm}$ for $Pe>0$ starts at
a value $>1$ and quickly [beginning of the darker shade in Fig. 3(a)]
approaches $1$, consistent with the long-time diffusive behavior.

For $Pe\geq 50$, where the polymer is in a coiled state, starting from a value
around $0.9$, the exponent $\alpha_{\rm cb}$ drops significantly before it
climbs up in the diffusive regime [Fig. 3(b)]. However, one can hardly see a
flat intermediate region to consider this as a true scaling regime.
Importantly, the data never really show a steady behavior around the value
$2/3$, thus ruling out the apparent Zimm’s scaling. This drop in $\alpha_{\rm
cb}$ can rather be interpreted as an effect of gradual crossover to the long-
time diffusive regime. The crossover gets delayed with increase in $Pe$, as
evident from Figs. 2(d)-(f) showing that the data for ${\rm MSD}_{\rm cb}(t)$
merge with the one for ${\rm MSD}_{\rm cm}(t)$ at large $t$.

Figure 4:  Chain-length dependence of the mean square displacement of the
center of mass at temperature $T=0.1\epsilon/k_{B}$ for (a) $Pe=0$, (b)
$Pe=25$, and (c) $Pe=62.5$. Regions with darker shades mark the time period
$t\in[10^{2}:10^{3}]$, over which the diffusion coefficient $D_{\rm eff}$ is
calculated using Eq. (7). (d) Scaling of $D_{\rm eff}$ with $N$ for different
$Pe$. The dashed lines represent the prediction in Eq. (6) with $D_{m}=0.42$
for fixed $Pe$. (e) $D_{\rm eff}$ as a function of $Pe$ for different $N$. The
dashed lines there also represent the prediction in Eq. (6) for fixed $N$.

For the end beads $\alpha_{\rm eb}(t)$ shows a similar behavior of approaching
$1$ at late times [Fig. 3(c)]. This implies that the data for ${\rm MSD}_{\rm
eb}(t)$ must coincide with ${\rm MSD}_{\rm cm}(t)$ at large $t$, which can be
verified from the plots in Figs. 2(d), (e), and (f). Similar to $\alpha_{\rm
cb}$, at intermediate times the data for $\alpha_{\rm eb}$ show a drop from
$1$ and tend to become flat before it finally approaches $1$ at large $t$.
This indicates the presence of a true intermediate power-law regime. For
$Pe\geq 50$, the value of $\alpha_{\rm eb}$ in the intermediate flat regime is
less than $2/3$. In absence of hydrodynamics effects, i.e., for a Rouse
polymer with excluded volume, in the intermediate regime, the scaling for the
end monomers is given by Paul  _et al._ (1991); Dünweg and Kremer (1991,
1993); Shusterman  _et al._ (2004)

${\rm MSD}_{\rm eb}(t)\sim t^{2\nu/(1+2\nu)}.$ (4)

For a Gaussian chain having $\nu=1/2$ this provides a $\sim t^{1/2}$ behavior.
In the present case at large $Pe$, the polymer behaves like a self-avoiding
coil with $\nu\approx 3/5$ producing a scaling $\sim t^{6/11}$. Our data is
indeed consistent with such a behavior, shown by the dashed lines in Figs.
2(c) and (f). Thus it can be inferred that the intermediate Rouse scaling
which in general does not hold for a passive polymer in poor solvent Milchev
and Binder (1994), can be recovered in an active polymer at sufficiently large
strength of activity.

For a theoretical understanding, we consider an analog of the Rouse model, the
simplest model describing dynamics of a passive polymer in absence of
hydrodynamics Rouse Jr (1953). In the model successive active monomers along
the chain are connected via harmonic springs. In addition, each monomer
experiences a net random force $\vec{F}_{i}$ resulting from the thermal noise,
active force, and nonbonded interaction in poor solvent condition. Assuming
that $\vec{F}_{i}$ is Delta-correlated over time and space, for the motion of
the cm at long time one can write down (see details in the Supplemental
Material Note (1))

$\begin{split}&{\rm MSD}_{\rm cm}(t)=\langle\left[\vec{r}_{\rm
cm}(t)-\vec{r}_{\rm cm}(0)\right]^{2}\rangle\\\
&=\left\langle\int_{0}^{t}d{t^{\prime}}\int_{0}^{t}d{t^{\prime\prime}\left[\frac{1}{N}\sum_{i=1}^{N}\vec{F}_{i}(t^{\prime})\right]\cdot\left[\frac{1}{N}\sum_{j=1}^{N}\vec{F}_{j}(t^{\prime\prime})\right]}\right\rangle\\\
&=6D_{\rm eff}t=\frac{6(D_{\rm a}/D_{m})}{N}t.\end{split}$ (5)

Here, $D_{\rm eff}$ is effective diffusion constant of the cm of the polymer
which is related to the diffusion constant of the constituent active monomers
$D_{\rm a}$ via a modification factor $D_{m}$, introduced to take into account
the poor solvent condition. From the long-time ${\rm MSD}$ of an active
Brownian particle Howse  _et al._ (2007), one gets $D_{\rm
a}=\left(1+{Pe^{2}}/{18}\right){k_{B}T}/{\gamma}$, which on inserting in Eq.
(5) yields (see Supplemental Material Note (1))

$D_{\rm eff}=\left(1+\frac{Pe^{2}}{18}\right)\frac{k_{B}T}{\gamma D_{m}N},$
(6)

implying a Rouse-like scaling $D_{\rm eff}\sim N^{-1}$ at a fixed $Pe$, and
$D_{\rm eff}\sim Pe^{2}$ for fixed $N$.

To verify the prediction in Eq. (6), we calculate $D_{\rm eff}$ from our
simulation data using the following prescription

$D_{\rm eff}=\frac{1}{6}\lim_{t\rightarrow\infty}\frac{d}{dt}{\rm MSD}_{\rm
cm}(t).$ (7)

The data of ${\rm MSD}_{\rm cm}(t)$ for polymers of different length $N$,
showing a linear behavior in the long-time limit, are presented in Figs.
4(a)-(c). The extracted $D_{\rm eff}$ as a function of $N$ for fixed $Pe$ are
presented in Fig. 4(d). The dashed lines represent Eq. (6) with $D_{m}=0.42$
which was obtained as the most reasonable choice following a rigorous fitting
exercise presented in the Supplemental Material Note (1). The consistency of
our data with the plotted functions irrespective of the value of $Pe$ confirms
the presence of a universal Rouse-like scaling as embedded in the prediction
(6). Furthermore, as $Pe$ increases a significant enhancement of $D_{\rm eff}$
can be noticed. Figure 4(e) shows this enhancement via plots of $D_{\rm eff}$
as a function $Pe$ for fixed $N$. The functional dependence of $D_{\rm eff}$
on $Pe$ for a fixed $N$ is predicted in Eq. (6). Fitting the ansatz (presented
in the Supplemental Material Note (1)) using $D_{\rm eff}$ vs. $Pe$ data yield
the same $D_{m}=0.42$. Plots of Eq. (6) with $D_{m}=0.42$, shown by the dashed
lines in Fig. 4(e) not only depict an unambiguous agreement of the prediction
with the simulation data but also indicates that the modification factor
$D_{m}$ is rather universal, independent of $N$ and $Pe$.

In conclusion, we have presented results for the steady-state dynamics of an
active Brownian polymer in poor solvent condition. In order to explore the
dynamics we have monitored the motions of the center of mass, central monomer,
and the end monomers. The mean square displacement of the end monomers shows
the presence of an intermediate regime. In the large-activity limit, this
intermediate regime exhibits a $\sim t^{6/11}$ scaling, which generally holds
for Rouse dynamics of a passive polymer in good solvent. In the long-time
limit, the mean square displacement of the central and end monomers merge with
that of the diffusive behavior of the center of mass. This allows us to
estimate the long-time effective diffusion coefficient $D_{\rm eff}$ of the
polymer. Analytically, we predict the dependence of $D_{\rm eff}$ on chain
length $N$ and activity strength $Pe$ using a Rouse model of active polymer.
Our numerical results are in agreement with the theoretical prediction showing
a significant enhancement of $D_{\rm eff}$ as a function of $Pe$ obeying a
scaling $D_{\rm eff}\sim Pe^{2}$. Similarly as predicted, the data show that
the universal Rouse-like $D_{\rm eff}\sim N^{-1}$ still holds strongly
irrespective of $Pe$. It would be interesting to explore the robustness of
this Rouse-like behavior for semiflexible polymers with activity 222S.
Majumder, S. Paul, and W. Janke, in progress..

This work is the first to explore the steady-state dynamics of an active
polymer in poor solvent. As a future endeavour, it would be worth to
investigate the same in other active polymer models. Our main result showing
activity induced enhanced diffusion of a polymer in a poor solvent condition
might indulge in design of synthetic active polymers which potentially can be
employed in delivering drugs for a wide variety of media. In connection, it
would also be intriguing to study the effect of hydrodynamics and explicit
solvent on this apparently universal dynamics of active polymer Majumder  _et
al._ (2019); Bera  _et al._ (2022).

###### Acknowledgements.

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Grant No. 189 853 844 – SFB/TRR 102 (Project B04)
and further supported by the Deutsch-Französische Hochschule (DFH-UFA) through
the Doctoral College “$\mathbb{L}^{4}$” under Grant No. CDFA-02-07, and the
Leipzig Graduate School of Natural Sciences “BuildMoNa”. S.M. thanks the
Science and Engineering Research Board (SERB), Govt. of India for a Ramanujan
Fellowship (file no. RJF/2021/000044). S.P. acknowledges ICTS-TIFR, DAE, Govt.
of India for a research fellowship.

## References

  * Wu and Libchaber (2000) X.-L. Wu and A. Libchaber, “Particle diffusion in a quasi-two-dimensional bacterial bath,” Phys. Rev. Lett. 84, 3017–3020 (2000).
  * Park  _et al._ (2010) Y. Park, C.A. Best, K. Badizadegan, R.R. Dasari, M.S. Feld, T. Kuriabova, M.L. Henle, A.J. Levine, and G. Popescu, “Measurement of red blood cell mechanics during morphological changes,” Proc. Nat. Acad. Sci. 107, 6731–6736 (2010).
  * Ben-Isaac  _et al._ (2011) E. Ben-Isaac, Y. Park, G. Popescu, F.L.H. Brown, N.S. Gov, and Y. Shokef, “Effective temperature of red-blood-cell membrane fluctuations,” Phys. Rev. Lett. 106, 238103 (2011).
  * Jee  _et al._ (2018a) A.-Y. Jee, Y.-K. Cho, S. Granick, and T. Tlusty, “Catalytic enzymes are active matter,” Proc. Nat. Acad. Sci. 115, E10812–E10821 (2018a).
  * Jee  _et al._ (2018b) A.-Y. Jee, S. Dutta, Y.-K. Cho, T. Tlusty, and S. Granick, “Enzyme leaps fuel antichemotaxis,” Proc. Nat. Acad. Sci. 115, 14–18 (2018b).
  * Ramaswamy (2010) S. Ramaswamy, “The mechanics and statistics of active matter,” Annu. Rev. Condens. Matter Phys. 1, 323–345 (2010).
  * Shaebani  _et al._ (2020) M.R. Shaebani, A. Wysocki, R.G. Winkler, G. Gompper, and H. Rieger, “Computational models for active matter,” Nat. Rev. Phys. 2, 181–199 (2020).
  * Elgeti  _et al._ (2015) J. Elgeti, R.G. Winkler, and G. Gompper, “Physics of microswimmers—single particle motion and collective behavior: a review,” Rep. Prog. Phys. 78, 056601 (2015).
  * Winkler and Gompper (2020) R.G. Winkler and G. Gompper, “The physics of active polymers and filaments,” J. Chem. Phys. 153, 040901 (2020).
  * Kaiser  _et al._ (2015) A. Kaiser, S. Babel, B. ten Hagen, C. von Ferber, and H. Löwen, “How does a flexible chain of active particles swell?” J. Chem. Phys. 142, 124905 (2015).
  * Isele-Holder  _et al._ (2015) R.E. Isele-Holder, J. Elgeti, and G. Gompper, “Self-propelled worm-like filaments: spontaneous spiral formation, structure, and dynamics,” Soft Matter 11, 7181–7190 (2015).
  * Isele-Holder  _et al._ (2016) R.E. Isele-Holder, J. Jäger, G. Saggiorato, J. Elgeti, and G. Gompper, “Dynamics of self-propelled filaments pushing a load,” Soft Matter 12, 8495–8505 (2016).
  * Bianco  _et al._ (2018) V. Bianco, E. Locatelli, and P. Malgaretti, “Globulelike conformation and enhanced diffusion of active polymers,” Phys. Rev. Lett. 121, 217802 (2018).
  * Locatelli  _et al._ (2021) E. Locatelli, V. Bianco, and P. Malgaretti, “Activity-induced collapse and arrest of active polymer rings,” Phys. Rev. Lett. 126, 097801 (2021).
  * Das  _et al._ (2021) S. Das, N. Kennedy, and A. Cacciuto, “The coil–globule transition in self-avoiding active polymers,” Soft Matter 17, 160–164 (2021).
  * Paul  _et al._ (2021) S. Paul, S. Majumder, and W. Janke, “Motion of a polymer globule with Vicsek-like activity: From super-diffusive to ballistic behavior,” Soft Mater. 19, 306–315 (2021).
  * Paul  _et al._ (2022a) S. Paul, S. Majumder, S.K. Das, and W. Janke, “Effects of alignment activity on the collapse kinetics of a flexible polymer,” Soft Matter 18, 1978–1990 (2022a).
  * Paul  _et al._ (2022b) S. Paul, S. Majumder, and W. Janke, “Role of temperature and alignment activity on kinetics of coil-globule transition of a flexible polymer,” J. Phys.: Conf. Ser. 2207, 012027 (2022b).
  * Paul  _et al._ (2022c) S. Paul, S. Majumder, and W. Janke, “Activity mediated globule to coil transition of a flexible polymer in a poor solvent,” Soft Matter 18, 6392–6403 (2022c).
  * Paxton  _et al._ (2004) W.F. Paxton, K.C. Kistler, C.C. Olmeda, A. Sen, S.K. St. Angelo, Y. Cao, T.E. Mallouk, P.E. Lammert, and V.H. Crespi, “Catalytic nanomotors: Autonomous movement of striped nanorods,” J. Am. Chem. Soc. 126, 13424–13431 (2004).
  * Sanchez  _et al._ (2011) S. Sanchez, A.A. Solovev, S.M. Harazim, and O.G. Schmidt, “Microbots swimming in the flowing streams of microfluidic channels,” J. Am. Chem. Soc. 133, 701–703 (2011).
  * Ebbens  _et al._ (2012) S. Ebbens, M.-H. Tu, J.R. Howse, and R. Golestanian, “Size dependence of the propulsion velocity for catalytic Janus-sphere swimmers,” Phys. Rev. E 85, 020401 (2012).
  * Buttinoni  _et al._ (2013) I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger, and T. Speck, “Dynamical clustering and phase separation in suspensions of self-propelled colloidal particles,” Phys. Rev. Lett. 110, 238301 (2013).
  * Kümmel  _et al._ (2013) F. Kümmel, B. Ten Hagen, R. Wittkowski, I. Buttinoni, R. Eichhorn, G. Volpe, H. Löwen, and C. Bechinger, “Circular motion of asymmetric self-propelling particles,” Phys. Rev. Lett. 110, 198302 (2013).
  * Zöttl and Stark (2014) A. Zöttl and H. Stark, “Hydrodynamics determines collective motion and phase behavior of active colloids in quasi-two-dimensional confinement,” Phys. Rev. Lett. 112, 118101 (2014).
  * De Gennes (1979) P.-G. De Gennes,  _Scaling Concepts in Polymer Physics_ (Cornell University Press, Ithaca and London, 1979).
  * Doi (1996) M. Doi,  _Introduction to Polymer Physics_ (Oxford University Press, New York, 1996).
  * Rubinstein and Colby (2003) M. Rubinstein and R.H. Colby,  _Polymer Physics_ , Vol. 23 (Oxford University Press, New York, 2003).
  * Harder  _et al._ (2014) J. Harder, C. Valeriani, and A. Cacciuto, “Activity-induced collapse and reexpansion of rigid polymers,” Phys. Rev. E 90, 062312 (2014).
  * Kaiser and Löwen (2014) A. Kaiser and H. Löwen, “Unusual swelling of a polymer in a bacterial bath,” J. Chem. Phys. 141, 044903 (2014).
  * Chaki and Chakrabarti (2019) S. Chaki and R. Chakrabarti, “Enhanced diffusion, swelling, and slow reconfiguration of a single chain in non-Gaussian active bath,” J. Chem. Phys. 150, 094902 (2019).
  * Liu  _et al._ (2019) X. Liu, H. Jiang, and Z. Hou, “Configuration dynamics of a flexible polymer chain in a bath of chiral active particles,” J. Chem. Phys. 151, 174904 (2019).
  * Anderson  _et al._ (2022) C.J. Anderson, G. Briand, O. Dauchot, and A. Fernández-Nieves, “Polymer-chain configurations in active and passive baths,” Phys. Rev. E 106, 064606 (2022).
  * Goswami  _et al._ (2022) K. Goswami, S. Chaki, and R. Chakrabarti, “Reconfiguration, swelling and tagged monomer dynamics of a single polymer chain in Gaussian and non-Gaussian active baths,” J. Phys. A: Math. Theor. 55, 423002 (2022).
  * Biswas  _et al._ (2017) B. Biswas, R.K. Manna, A. Laskar, P.B.S. Kumar, R. Adhikari, and G. Kumaraswamy, “Linking catalyst-coated isotropic colloids into “active” flexible chains enhances their diffusivity,” ACS nano 11, 10025–10031 (2017).
  * McMullen  _et al._ (2018) A. McMullen, M. Holmes-Cerfon, F. Sciortino, A.Y. Grosberg, and J. Brujic, “Freely jointed polymers made of droplets,” Phys. Rev. Lett. 121, 138002 (2018).
  * Rouse Jr (1953) P.E. Rouse Jr, “A theory of the linear viscoelastic properties of dilute solutions of coiling polymers,” J. Chem. Phys. 21, 1272–1280 (1953).
  * Zimm (1956) B.H. Zimm, “Dynamics of polymer molecules in dilute solution: Viscoelasticity, flow birefringence and dielectric loss,” J. Chem. Phys. 24, 269–278 (1956).
  * Naghizadeh and Kovac (1987) J. Naghizadeh and J. Kovac, “Cubic-lattice simulation of the dynamics of a collapsed polymer chain,” Phys. Rev. Lett. 59, 1710–1713 (1987).
  * Milchev  _et al._ (1993) A. Milchev, W. Paul, and K. Binder, “Off-lattice Monte Carlo simulation of dilute and concentrated polymer solutions under theta conditions,” J. Chem. Phys. 99, 4786–4798 (1993).
  * Milchev and Binder (1994) A. Milchev and K. Binder, “Anomalous diffusion and relaxation of collapsed polymer chains,” Europhys. Lett. 26, 671–676 (1994).
  * Fily and Marchetti (2012) Y. Fily and M.C. Marchetti, “Athermal phase separation of self-propelled particles with no alignment,” Phys. Rev. Lett. 108, 235702 (2012).
  * Majumder  _et al._ (2017) S. Majumder, J. Zierenberg, and W. Janke, “Kinetics of polymer collapse: Effect of temperature on cluster growth and aging,” Soft Matter 13, 1276–1290 (2017).
  * Note (1) See Supplemental Material at […] which includes Refs. Rouse Jr (1953); Howse  _et al._ (2007). The Supplemental Material contains a detailed description of the model and method of simulation used, and a few representative steady-state trajectories of the polymer. It also includes details of the analogous Rouse model for active polymer in poor solvent, and a derivation of the corresponding diffusion constant of the center of mass. Furthermore, we present details of the performed fitting exercise showing the agreement between theory and simulation data.
  * Paul  _et al._ (1991) W. Paul, K. Binder, D.W. Heermann, and K. Kremer, “Dynamics of polymer solutions and melts. Reptation predictions and scaling of relaxation times,” J. Chem. Phys. 95, 7726–7740 (1991).
  * Dünweg and Kremer (1991) B. Dünweg and K. Kremer, “Microscopic verification of dynamic scaling in dilute polymer solutions: A molecular-dynamics simulation,” Phys. Rev. Lett. 66, 2996 (1991).
  * Dünweg and Kremer (1993) B. Dünweg and K. Kremer, “Molecular dynamics simulation of a polymer chain in solution,” J. Chem. Phys. 99, 6983–6997 (1993).
  * Hinczewski  _et al._ (2009) M. Hinczewski, X. Schlagberger, M. Rubinstein, O. Krichevsky, and R.R. Netz, “End-monomer dynamics in semiflexible polymers,” Macromolecules 42, 860–875 (2009).
  * Shusterman  _et al._ (2004) R. Shusterman, S. Alon, T. Gavrinyov, and O. Krichevsky, “Monomer dynamics in double-and single-stranded DNA polymers,” Phys. Rev. Lett. 92, 048303 (2004).
  * Howse  _et al._ (2007) J.R. Howse, R.A.L. Jones, A.J. Ryan, T. Gough, R. Vafabakhsh, and R. Golestanian, “Self-motile colloidal particles: From directed propulsion to random walk,” Phys. Rev. Lett. 99, 048102 (2007).
  * Note (2) S. Majumder, S. Paul, and W. Janke, in progress.
  * Majumder  _et al._ (2019) S. Majumder, H. Christiansen, and W. Janke, “Dissipative dynamics of a single polymer in solution: A Lowe-Andersen approach,” J. Phys.: Conf. Ser. 1163, 012072 (2019).
  * Bera  _et al._ (2022) A. Bera, S. Sahoo, S. Thakur, and S.K. Das, “Active particles in explicit solvent: Dynamics of clustering for alignment interaction,” Phys. Rev. E 105, 014606 (2022).