colorlinks, linkcolor=red!50!black, citecolor=blue!50!black,
urlcolor=blue!80!black

# Topological inference on brain networks across subtypes of post-stroke
aphasia

Yuan Wang1∗ Jian Yin2 Rutvik H. Desai3 1Department of Epidemiology and
Biostatistics, University of South Carolina, U.S.A.
2Department of Biostatistics, Nanjing Medical University, China.
3Department of Psychology, University of South Carolina, U.S.A.
∗Correspondence<EMAIL_ADDRESS>
###### Abstract

Persistent homology (PH) characterizes the shape of brain networks through the
persistence features. Group comparison of persistence features from brain
networks can be challenging as they are inherently heterogeneous. A recent
scale-space representation of persistence diagram (PD) through heat diffusion
reparameterizes using the finite number of Fourier coefficients with respect
to the Laplace-Beltrami (LB) eigenfunction expansion of the domain, which
provides a powerful vectorized algebraic representation for group comparisons
of PDs. In this study, we advance a transposition-based permutation test for
comparing multiple groups of PDs through the heat-diffusion estimates of the
PDs. We evaluate the empirical performance of the spectral transposition test
in capturing within- and between-group similarity and dissimilarity with
respect to statistical variation of topological noise and hole location. We
also illustrate how the method extends naturally into a clustering scheme by
subtyping individuals with post-stroke aphasia through the PDs of their
resting-state functional brain networks.

###### keywords:

Topological data analysis; Persistent homology; Permutation test; Brain
network.

##  1 Introduction

Brain network modeling based on magnetic resonance imaging (MRI) is an
effective approach to understand the functions and dysfunctions of the brain.
Brain networks have an innate graph structure that have been traditionally
studied through graphical or graph theoretic models based on single-scale
covariance estimation (S. Huang ., 2010) or single-scale graph-theoretic
measures (Sporns, 2002; Rubinov Sporns, 2010). These models effectively
characterize brain network topology and have become the norm for brain network
analysis. However, it has recently been noted that single-scale models may not
be sufficient in capturing the complexity of brain connectivity and multi-
scale models are needed (Betzel Bassett, 2017). On the other hand, a
ubiquitous problem in brain network analysis is selection of threshold on edge
weights to reveal significant connections within and between brain regions.
Arbitrary threshold may cause problem of bias and consistency across studies
(Drakesmith ., 2015; Garrison ., 2015). A multi-scale approach to brain
network modeling has thus become widely adapted through persistent homology
(PH), which captures multi-scale features of data through invariant
topological structures (Edelsbrunner ., 2002). Using the language of
simplicial homology (Hatcher, 2001), PH reveals the underlying topological
structures of data by their persistence through a dynamic assortment of
points, edges, and triangles. The fact that the overall topological changes
hold more significance over fleeting structures in PH makes the algorithm
particularly robust under the presence of noise and artifacts, thus revealing
more insight on network topology than single-scale measures (Carlsson, 2009).
Topological characteristics of the dynamic changes through the PH process are
summarized through persistence features.

Current persistence features consist of barcode and persistence diagram (PD),
the original descriptors proposed by Edelsbrunner . (2002), and persistence
landscape (PL) (Bubenik, 2015) and persistence image (PI) (Adams ., 2017),
both of which were developed when the demand increased for incorporating
persistence features in statistical inference and machine learning models.
Persistence features are inherently heterogeneous for noisy samples, even when
the samples come from homogeneous underlying data objects. The heterogeneous
nature of persistence features means that statistical inference for group
comparison is not straightforward. Parametric inference often requires
stringent distributional assumptions, which are rarely met by persistence
features. So we utilize a nonparametric inference approach. Permutation
testing is a standard nonparametric inference procedure for complex data
objects and features without clear distributional properties. It is known as
the exact test in statistics since the distribution of the test statistic
under the null hypothesis can be exactly computed if we can calculate all the
test statistics under every possible permutation. It is thus one of the most
widely used inference procedures in neuroimaging studies where the data is
typically complex in structure and the underlying distributional properties
are difficult to quantify (Nichols Holmes, 2002; Simpson ., 2013; Winkler .,
2016). However, generating every possible permutation for brain networks with
a large number of nodes is still extremely time consuming even for a modest
sample size. Standard permutation testing through approximations only reaches
a fraction of the exhaustive list of permutations and is computationally
intensive. When the total number of permutations is large, various resampling
techniques have been proposed to speed up the computation in the past (Nichols
Holmes, 2002; Winkler ., 2016). These resampling methods generate a small
fraction of possible permutations and the statistical significance is computed
approximately. Neuroimaging studies typically generate 5,000–1,000,000
permutations, less than a fraction of all possible permutations. A few
approaches have been developed to overcome the computational bottleneck for
permutation testing on persistence features. The exact topological inference
approach allows for fast permutation of monotone functions built on birth or
death times in barcodes with respect to the Komogorov-Smirnov (KS) distance
(Chung, Lee ., 2019). This approach has quadratic run time that beats the
exponential run time of standard permutation tests and has been extended to
compare PLs (Wang ., 2019, 2021). However, the approach is limited to
comparing two features and not applicable for comparing between two sets of
features. Another rapid permutation test based on transpositions does not
require monotonicity and is applicable for comparing two sets of persistence
features (Chung, Xie ., 2019; Songdechakraiwut Chung, 2023). It has allowed us
to develop a unified framework for topological inference through heat kernel
estimation of PDs.

Inference and learning approaches comparing PDs have been built on confidence
band (Fasy ., 2014) to functional representations (Chung ., 2009; Pachauri .,
2011; Bubenik, 2015; Reininghaus ., 2015; Carriére ., 2015; Chen ., 2015;
Adams ., 2017), as comparing raw PDs consisting of planar scatter points
encoding birth and death times of topological structures often require point
matching through, for instance, the Hungarian matching algorithm, which
quickly becomes formidable for large-scale data. It is also unclear how we may
compare two sets of raw PDs. The functional representation approach overcomes
the issue of the points on raw PDs having arbitrary locations and provides an
effective framework for downstream comparison. In this approach, PDs
essentially undergo a smoothing process, in some cases through a scale-space
representation from kernels for heat diffusion of Dirac delta functions
uniquely representing the points of PD (Reininghaus ., 2015). However,
existing kernel features on PD are typically convoluted, which lacks
flexibility when performing resampling-based statistical inference procedures
such as permutation testing. A new scale-space representation of PD was
recently proposed based on the heat kernel (HK) estimation (Kulkarni ., 2020),
where the upper-triangular domain of PDs is represented using a finite number
of Fourier coefficients with respect to the Laplace-Beltrami (LB)
eigenfunction expansion of the domain. It provides a powerful vectorized
algebraic representation for comparisons of PDs at the same coordinates,
foregoing the need for matching across PDs due to their arbitrary point
locations. Motivated by a topology-preserving spectral permutation test (Wang
., 2018), we developed an inference procedure for comparing two sets of PDs
estimated by the new scale-space representation by transposing the PD labels
(Wang ., 2022). By updating only the terms in an $L_{2}$-distance between the
mean HK estimates of two sets of PDs involved in each transposition,
computation becomes much faster than standard permutation testing that
exchanges an arbitrary number of labels in each iteration. This inference
procedure generalizes the method developed by Wang . (2018) for comparing
persistence features of single-trial univariate signals, where the resampling
takes place at the signal level and thus cannot be directly applied to images
and networks. The inference framework now resamples at the feature level,
which allows us to compare PDs of images and networks. We have also extended
it to a new topological ANOVA (T-ANOVA) approach to compare across multiple
groups of PDs without dimensionality reduction, as well as a topological
clustering scheme in application.

In this study, we establish a topological inference framework through
stability of HK estimation on PDs. We evaluate the empirical performance of
the spectral permutation test and T-ANOVA in simulation studies in detecting
heterogeneous topological noise and hole location across multiple images. We
also apply the methods to study topological difference in brain networks
across subtypes of individuals with post-stroke aphasia.

##  2 Methods

Brain networks are typically modeled as a weighted graph, with the edge
weights given by a similarity measure between the measurements on the nodes of
the network (Bassett, 2006; Bien Tibshirani, 2011). Suppose we have a network
represented by the weighted graph $G=(V,w)$ with the node set
$V=\\{1,\dots,p\\}$ and unique positive undirected edge weights $w=(w_{ij})$
constructed from a similarity measure such as Pearson’s correlation. We define
the binary network $G_{\epsilon}=(V,w_{\epsilon})$ as a subgraph of $G$
consisting of the node set $V$ and the binary edge weights $w_{\epsilon}$
defined by

$w_{ij,\epsilon}=\left\\{\begin{array}[]{cc}1&\text{if}~{}w_{ij}<\epsilon;\\\
0&\text{otherwise}.\end{array}\right.$ (1)

As we increase $\epsilon$, which we call the filtration value, more edges are
included in the binary network $G_{\epsilon}$ and so the size of the edge set
increases. Since edges connected in the network do not get disconnected again,
we observe a sequence of nested subgraphs

$G_{\epsilon_{0}}\subset G_{\epsilon_{1}}\subset
G_{\epsilon_{2}}\subset\cdots,$ (2)

for any

$\epsilon_{0}\leq\epsilon_{1}\leq\epsilon_{2}\leq\cdots.$

This sequence of nested subgraphs make up a Rips filtration where two nodes
with a weight $w_{ij}$ smaller than $\epsilon$ are connected, and the birth
and death of homological features in the form of clusters of nodes and holes
formed by more than 3 edges are tracked through the filtration (Lee ., 2011,
2014). We pair the birth and death times of clusters and holes as the
coordinates of scatter points on a planar graph
$\\{(a_{i},b_{i})\\}_{i=1}^{L}$ in the persistence diagram (PD). The
persistence of clusters and holes is measured by the drop from their
corresponding points to the $y=x$ line on the PD. Long persistence indicates
that the corresponding cluster or hole is more likely to be an underlying
feature in the network. As an illustration in Figure 1, we see how a point
that corresponds to a hole in a key shape stands out with high persistence in
the PD from the Rips filtration constructed on a 100-point point cloud sampled
from a key shape with a hole.

Figure 1: Left Three: The evolving 1-skeleton of a 100-point point cloud
sampled from a key shape with a hole. Right: PD from the Rips filtration
constructed on the 1-skeletons of the point cloud. The point in the PD that
corresponds to the key hole stands out with high persistence - much further
away from the diagonal ($y=x$) line than the rest of the points.

###  2.1 Heat kernel representation of persistence diagram

Since PDs do not form a vector space, they do not possess a natural
statistical framework (Bubenik, 2015) and requires additional manipulation
such as kernel smoothing for downstream statistical analysis. As with all
noisy data, smoothing is needed for reducing noise (typically random, often
artifactual) to better reveal the underlying data structure. We could either
smooth data used to construct the networks or smooth persistence descriptors
such as PD. In principle, topological noise and artifacts should be better
removed with the latter approach as persistence descriptors are designed to
capture topological structures, be they inherent or transient. Another
important reason for smoothing PDs is that the heterogenous nature of raw PDs
makes it difficult to perform various algebraic operations for statistical
inference. Various smoothing methods have been applied to PDs such that
statistical inference can be directly performed on them. Beginning with the
work of Chung . (2009), each PD is discretized using the the uniform square
grid and a concentration map is then obtained by counting the number of points
in each pixel, which is equivalent to smoothing PD with a uniform kernel. This
approach is somewhat similar to the voxel-based morphometry (Ashburner
Friston, 2000), where brain tissue density maps are used as a shapeless metric
for characterizing concentration of the amount of tissue. Pachauri . (2011)
followed up the approach by smoothing the PD by a Gaussian kernel centered at
every point. Later, Bubenik (2015) proposed the persistence descriptor PL by
representing the PD as a function in the Banach space $L_{p}(\mathbb{R}^{2})$
aimed at statistical analysis. PL is easily invertible to a PD, but
overemphasizes the high-persistence features. To account for the overall
pattern of persistence features, a persistence scale-space (PSS) kernel
approach was then proposed by Reininghaus . (2015), where the points in PD are
treated as heat sources modeled as Dirac-delta functions and used as an
initial condition for a heat diffusion problem with a Dirichlet boundary
condition on the diagonal. The closed-form solution of the diffusion problem
is an $L_{2}(\Omega)$ function obtained by convolving the initial condition
with a Gaussian kernel, with $\Omega=\\{(x,y)\in\mathbb{R}^{2}:y\geq x\\}$
being the closed half plane above the diagonal line $y=x$, and the feature map
from the PDs to $L_{2}(\Omega)$ at a fixed scale yields the PSS kernel. The
Hilbert space structure of $L_{2}(\mathbb{R}^{2})$ can be used to construct a
PL kernel similar to PSS (Reininghaus ., 2015). The relatively new persistence
descriptor PI sampled at discrete uniform grid to produce homogenous
vectorized data out of PDs (Adams ., 2017). PIs live in Euclidean space and
are therefore amenable to a broader range of learning techniques than PLs
(Adams ., 2017). A new heat kernel representation for PDs has recently been
proposed by Kulkarni . (2020), which not only allows a non-convoluted
vectorized representation for comparisons at the same coordinates of PDs but
also smoothing PD at different scales. It has also been combined with
transposition test, a novel permutation testing approach, for fast inference
on PDs (Wang ., 2022). We provide in the next two sections a detailed
description of heat kernel representation and transposition test on PDs.

Heat kernel representation has been established as a smoothing framework for
noisy measurements on a general manifold $\mathcal{M}\subset\mathbb{R}^{d}$
(Chung ., 2007, 2014). We assume the fundamental stochastic model

$f(p)=h(p)+\varepsilon(p),p\in\mathcal{M},$ (3)

where $f$ is the noisy measurement, $h$ is the unknown signal, and
$\varepsilon$ is a zero-mean Gaussian random field. We make the general enough
assumptions that $f\in L_{2}(\mathcal{M})$, the space of square integrable
functions on $\mathcal{M}$ with the inner product $\langle
f_{1},f_{2}\rangle=\int_{\mathcal{M}}f_{1}(p)f_{2}(p)d\mu(p)$, where $\mu$ is
the Lebesgue measure. A self-adjoint operator $\mathcal{L}$, i.e. $\langle
f_{1},\mathcal{L}f_{2}\rangle=\langle\mathcal{L}f_{1},f_{2}\rangle$ for all
$f_{1},f_{2}\in L_{2}(\mathcal{M})$, induces orthonormal eigenvalues
$\lambda_{k}$ and eigenfunctions $\psi_{k}$ on $\mathcal{M}$:

$\mathcal{L}\psi_{k}=\lambda_{k}\psi_{k},k=0,1,\dots,$ (4)

where, without loss of generality, we can sort the eigenvalues $\lambda_{k}$
such that

$0=\lambda_{0}\leq\lambda_{1}\leq\cdots,$

and the eigenfunctions $\psi_{k}$ can be numerically computed by solving a
generalized eigenvalue problem. Then, by Mercer’s Theorem, any symmetric
positive definite kernel can be written as

$K(p,q)=\sum_{k=0}^{\infty}\tau_{k}\psi_{k}(p)\psi_{k}(q)$ (5)

Now consider the diffusion-like Cauchy problem

$\frac{\partial
h(\sigma,p)}{\partial\sigma}+\mathcal{L}(\sigma,p)=0,p\in\mathcal{M},$ (6)

with the initial condition $h(\sigma=0,p)=f(p)$. The partial differential
equation (6) diffuses the noisy data $h(p)$ over $\sigma$. For the self-
adjoint operator $\mathcal{L}$, (6) has the unique solution (Chung ., 2007)

$h(\sigma,p)=\sum_{k=0}^{\infty}e^{-\lambda_{k}\sigma}\langle
h,\psi_{k}\rangle\psi_{k}(p),$ (7)

which provides an estimate $\hat{h}_{\sigma}(p)$ of the unknown signal $h(p)$.
The bandwidth $\sigma$ controls the amount of smoothing in the estimate; as
$\sigma$ increases, $\hat{h}_{\sigma}(p)$ becomes smoother. When $\mathcal{L}$
is the Laplace-Beltrami (LB) operator, the diffusion equation (6) becomes the
isotropic heat diffusion equation and the kernel (5) becomes the heat kernel
(HK)

$K_{\sigma}(p,q)=\sum_{k=0}^{\infty}e^{-\lambda_{k}\sigma}\psi_{k}(p)\psi_{k}(q),p,q\in\mathcal{M},$
(8)

where the $\psi_{k}$ are the eigenfunctions of the LB operator $\Delta$
satisfying

$\Delta\psi_{k}(p)=\lambda_{k}\psi_{k}(p)$

for $p\in\mathcal{M}$. The HK framework has been shown to be equivalent to
kernel regression and wavelet (Chung ., 2014).

To construct a HK representation of PD, we restrict the domain of diffusion to
$\mathcal{M}=\mathcal{T}=\\{(x,y)\in\mathbb{R}^{2}:y>x\\}$, i.e. the upper
triangular region above the diagonal line $y=x$ where the scatter points of
the PD $D=\\{(a_{i},b_{i})\\}_{i=1}^{P}$ are located. We constrain
$\mathcal{T}$ within a certain range, such as standardizing the coordinates of
the PD, so that $\mathcal{T}$ is bounded. Consider heat diffusion equation

$\frac{\partial h(\sigma,p)}{\partial\sigma}=\Delta h(\sigma,p)$ (9)

with the initial condition

$\quad h(\sigma=0,p)=\sum_{i=1}^{P}\delta_{(a_{i},b_{i})}(p),$

where $\delta_{(a_{i},b_{i})}$ is the Dirac-delta function at $(a_{i},b_{i})$.
The scatter points in the PD serve as the heat sources of the diffusion
process. To simplify notation, we will refer to any series $h(\sigma,p)$ as
$h_{\sigma}(p)$ as the bandwidth $\sigma$ is fixed. A unique solution to (9)
is given by the HK expansion

$\displaystyle h_{\sigma}(p)$ $\displaystyle=$
$\displaystyle\int_{\mathcal{T}}K_{\sigma}(p,q)h_{0}(q)\;d\mu(q)$ (10)
$\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}e^{-\lambda_{k}\sigma}f_{k}\psi_{k}(p),$

where

$K_{\sigma}(p,q)=\sum_{k=0}^{\infty}e^{-\lambda_{k}\sigma}\psi_{k}(p)\psi_{k}(q),p,q\in\mathcal{T},$
(11)

is the HK with respect to the eigenfunctions $\psi_{k}$ of the LB operator
$\Delta$ satisfying $\Delta\psi_{k}(p)=\lambda_{k}\psi_{k}(p)$ for
$p\in\mathcal{T}$, and

$f_{k}=\big{\langle}h_{0},\psi_{k}\big{\rangle}=\int_{\mathcal{T}}h_{0}(q)\psi_{k}(q)\;d\mu(q)\\\
=\sum_{i=1}^{P}\psi_{k}(a_{i},b_{i})$ (12)

are the Fourier coefficients with respect to the the LB eigenfunctions. The
first eigenvalue $\lambda_{0}=0$ of the LB operator corresponds to
eigenfunction $\psi_{0}=\frac{1}{\sqrt{\mu(\mathcal{T})}}$, where
$\mu(\mathcal{T})$ is the area of the triangular region $\mathcal{T}$ and
$\sigma$ is the bandwidth of the HK.

The HK expansion (10) provides a vectorized representation of the PD $D$ so
that we can compare across PDs at the same coordinates. In practice, we
include sufficiently large $\kappa$ terms to approximate the HK expansion:

$h^{\kappa}_{\sigma}(p)=\sum_{k=0}^{\kappa}e^{-\lambda_{k}\sigma}f_{k}\psi_{k}(p),$
(13)

which we refer to as the degree-$\kappa$ HK estimate of the given PD. When
$\sigma=0$, we can completely recover the initial scatter points. As
$\sigma\to\infty$, it is essentially smoothing the PD with a uniform kernel on
$\mathcal{T}$. Figure 2 shows the HK smoothing of a PD with respect to the
bandwidths $\sigma=0,0.1,1,10$. Note that the Fourier coefficients $f_{k}$
remain the same for all $k$ when constructing the HK expansion at different
diffusion scale $\sigma$.

Figure 2:  Heat kernel (HK) smoothing of the PD, constructed in Figure 1,
through Laplace-Beltrami (LB) eigenfunctions with respect to the bandwidths
$\sigma=0$ (original PD), 0.1, 1, 10. Top Row: Smoothed PDs. Bottom Row:
Corresponding Fourier coefficients with respect to the LB eigenfunctions
presented in matrix form.

As a distance measure for the HK-estimated PDs, we use the $L_{2}$-distance
between the functions $h^{1},h^{2}\in L_{2}(\mathcal{T})$ defined as

$\|h^{1}-h^{2}\|^{2}_{2}=\sum_{k=0}^{\infty}e^{-\lambda_{k}\sigma}(h^{1}_{k}-h^{2}_{k})^{2},$
(14)

where the $h^{1}_{k}$ and $h^{2}_{k},k=0,\dots,\infty,$ are the respective
Fourier coefficients of $h^{1}$ and $h^{2}$ as defined in (12) with respect to
the LB eigenfunctions.

In the standard kernel setup, we have the feature map

$\Phi_{\sigma}:\mathcal{D}\rightarrow L_{2}(\mathcal{T}),$

where $L_{2}(\mathcal{T})$ is the space of square integrable functions on
$\mathcal{T}$ with the $L_{2}$-distance between the functions $g^{1},g^{2}\in
L_{2}(\mathcal{T})$ defined as

$\|g^{1}-g^{2}\|^{2}_{2}=\sum_{k=0}^{\infty}e^{-\lambda_{k}\sigma}(g^{1}_{k}-g^{2}_{k})^{2},$
(15)

where the $g^{1}_{k}$ and $g^{2}_{k},k=0,\dots,\infty,$ are the respective
Fourier coefficients of $g^{1}$ and $g^{2}$ as defined in (12) with respect to
the LB eigenfunctions $\psi_{k},k=0,\dots,\infty$. Given bandwidth $\sigma>0$,

$\Phi_{\sigma}(D)=h_{\sigma}=\sum_{k=0}^{\infty}e^{-\lambda_{k}\sigma}f_{k}\psi_{k}(p),p\in
D,$

as defined in (10) for a PD $D\in\mathcal{D}$. This feature map corresponds to
the kernel

$K_{\sigma}(D_{1},D_{2})=\langle\Phi_{\sigma}(D_{1}),\Phi_{\sigma}(D_{2})\rangle_{L_{2}(\mathcal{T})},$

an explicit form of which is given by (11):

$K_{\sigma}(p,q)=\sum_{k=0}^{\infty}e^{-\lambda_{k}\sigma}\psi_{k}(p)\psi_{k}(q),p\in
D_{1},q\in D_{2}.$ (16)

We can show stability of the heat kernel

$\|K_{\sigma}*g^{1}-K_{\sigma}*g^{2}\|_{2}\leq\|g^{1}-g^{2}\|_{2}$ (17)

as follows: The integral version of Jensen’s inequality is

$\phi\left(\int w(x)\;dx\right)\leq\int\phi(w(x))\;dx$

for convex function $\phi$ (Matkowski Nikodem, 1994). Following Jensen’s
inequality,

$\displaystyle\|K_{\sigma}*g(p)\|_{2}^{2}$ $\displaystyle=$
$\displaystyle\int_{\mathcal{T}}\left|\int_{\mathcal{T}}K_{\sigma}(p,q)g(q)\;d\mu(q)\right|^{2}\;d\mu(p)$
(18) $\displaystyle\leq$
$\displaystyle\int_{\mathcal{T}}\int_{\mathcal{T}}K_{\sigma}(p,q)\|g(q)\|^{2}\;d\mu(q)\;d\mu(p)$
(19) $\displaystyle=$
$\displaystyle\int_{\mathcal{T}}\left|g(q)\right|^{2}\int_{\mathcal{T}}K_{\sigma}(p,q)\;d\mu(p)\;d\mu(q)$
(20) $\displaystyle=$
$\displaystyle\int_{\mathcal{T}}\left|g(q)\right|^{2}\;d\mu(q).$ (21)

We used the fact heat kernel $K_{\sigma}(p,q)$ is a probability distribution
such that

$\int_{\mathcal{T}}K_{\sigma}(p,q)\;d\mu(p)=1.$

Hence

$\|K_{\sigma}*g(p)\|_{2}\leq\|g(p)\|_{2}$

showing HK smoothing on PD is a contraction map (Chung ., 2018). Letting
$g=g^{1}-g^{2}$, we have the stability results. The HK smoothing reduces the
topological variability in PD.

We use a simple example with each of two PDs containing one of the two points
$(-\lambda,\lambda)$ and $(-\lambda+1,\lambda+1)$ (Reininghaus ., 2015), as an
illustration of the stability of the kernel smoothing procedures. When
comparing two PDs, the $L_{2}$-distance induced by the HK does not weigh over
any points in the PDs, as the distance between the two points is
$\sum_{k=0}^{\infty}e^{-\lambda_{k}\sigma}(\psi_{k}(-\lambda,\lambda)-\psi_{k}(-\lambda+1,\lambda+1))^{2}$,
which remains constant as $\lambda\rightarrow\infty$. In contrast, the PL-
induced kernel distance is dominated by variations in the points of high
persistence in the PDs, as the distance between the two points grows in the
order of $\sqrt{\lambda}$ and is unbounded, whereas the Wasserstein distance
and PSS-induced kernel distance do not over emphasize the high-persistence
points, as the distance between the two points asymptotically approach a
constant as $\lambda\rightarrow\infty$ (Reininghaus ., 2015). While the PSS
kernel representation, like our HK representation of PD, also uses an exact
solution to the heat diffusion problem with the original PD as the initial
condition (Figure 3), the implicit form of the solution is difficult to
manipulate for cost-effective resampling-based statistical inference. It is
likewise difficult to manipulate the Wasserstein distance and PL-induced
distance for the same purpose.

Figure 3: An example of a PD and PSS- and HK-estimated versions.

###  2.2 Permutation test on HK-estimated PDs

Existing kernel features on PD have been shown theoretically and empirically
to work well with machine learning frameworks (Reininghaus ., 2015; Adams .,
2017) but are typically convoluted, which lacks flexibility when performing
resampling-based statistical inference procedures such as permutation testing.
Our past studies have shown powerful applications of the series representation
of the heat diffusion problem, such as comparing the persistence features of
brain signals through built permutation test based on HK estimates of signals
(Wang ., 2018), where we studied how topology of signals is preserved by
permuting Fourier coefficients of sine and cosine basis functions. The
approach provides a ground for permutation testing based on spectral
components. The downside, however, is the computational load, with spectral
permutation of single-trial signals requiring hours on end to converge.

Here, we use the HK for PD smoothing and subsequent statistical inference
based on the HK-estimated PDs. Once we have the HK estimates of PD, we can use
them as the basis for statistical inference. Suppose we want to permute the
elements of two ordered sets with sizes $m$ and $n$

$\displaystyle{\bf x}$ $\displaystyle=$
$\displaystyle(x_{1},x_{2},\dots,x_{m}),$ $\displaystyle{\bf y}$
$\displaystyle=$ $\displaystyle(y_{1},y_{2},\dots,y_{n}),$

in a permutation test with the test statistic $f({\bf x},{\bf y})$. Under the
null hypothesis, we assume exchangeability of ${\bf x}$ and ${\bf y}$. Each
permutation is an unrestricted rearrangement of the combined ordered set ${\bf
z}=(x_{1},\dots,x_{m},y_{1},\dots,y_{n})$ and we denote all possible $(m+n)!$
permutations $\mathbb{S}_{m+n}$, which is a symmetric groups of order $m+n$.
The standard approximate permutation test typically used in practice is built
on on uniform sampling from the full set of permutations. The required number
of permutations for convergence increases exponentially as the sample sizes
increase. Even with sample sizes like $m=n=20$, the random permutation test
requires significant computational resources if we compute the test statistic
for each exchange of group labels.

A transposition is defined as a permutation $\pi_{ij}$ that exchanges the
$i$-th and $j$-th elements between ${\bf x}$ and ${\bf y}$ while keeping all
the other elements fixed, i.e.

$\displaystyle\pi_{ij}({\bf x})$ $\displaystyle=$
$\displaystyle(x_{1},\dots,x_{i-1},y_{j},x_{i+1},\dots,x_{m}),$
$\displaystyle\pi_{ij}({\bf y})$ $\displaystyle=$
$\displaystyle(y_{1},\dots,y_{j-1},x_{i},y_{j+1},\dots,y_{n}).$

Any permutation in $\mathbb{S}_{m+n}$ can be reached by a sequence of
transpositions (Chung, Xie ., 2019). The random transposition is a random walk
related to card shuffling problems and it is a special case of walk in
symmetric groups (Aldous, 1983; Aldous Diaconis, 1986). The walk between
elements within ${\bf x}$ or ${\bf y}$ is also allowed but will not affect the
computation a symmetric test functions. Instead of performing uniform random
sampling in $\mathbb{S}_{m+n}$, we can perform a sequence of random walks and
compute the test statistic at each walk. Consider walks in the two sample
setting. We will determine how test statistic changes over each walk. Over
random walk or transposition $\pi_{ij}$, the statistic changes from $L({\bf
x},{\bf y})$ to $L(\pi_{ij}({\bf x}),\pi_{ij}({\bf y}))$. Instead of computing
$L(\pi_{ij}({\bf x}),\pi_{ij}({\bf y}))$ directly, we can compute it from
$L({\bf x},{\bf y})$ incrementally in constant run time by updating the value
of $L({\bf x},{\bf y})$. If $L$ is an algebraic function that only involves
addition, subtraction, multiplication, division, integer exponents, there must
exists a function $M$ such that

$L(\pi_{ij}({\bf x}),\pi_{ij}({\bf y}))=M(L({\bf x},{\bf y}),x_{i},y_{i}),$

where the computational complexity of $g$ is constant (Chung, Xie ., 2019).
For instance, basic test statistics such as the two-sample $t$-statistic and
$F$-statistic are algebraic functions. If we take computation involving
fractional exponents as constant run time, then a much wider class of
statistics such as correlations can all have iterative formulation with
constant run time. In the case of computing two-sample $t$-statistic with $m$
and $n$ samples directly, we need to compute the sample means, which takes
$O(m)$ and $O(n)$ algebraic operations each. Then need to compute the sample
variances and pool them together, which requires $O(3m+2)$ and $O(3n+2)$
operations each. Combining the numerator and denominator in $t$-statistic
takes $O(16)$ operations. Thus, it takes total $O(4(m+n)+20)$ operations to
compute the $t$-statistic at each permutation. In general, by only updating
the terms in the test statistic affected by each transposition, the
transposition test would require considerably less computational resources
than the standard approximate permutation test.

When we compare two groups of PDs with sample sizes $m$ and $n$, we assume
under the null hypothesis that the functional means of the HK expansion of PDs
are the same for both groups, for a fixed bandwidth $\sigma>0$. The Fourier
coefficients in the HK expansion of population PDs in the two groups are
unknown. We estimate them with the HK expansion of sample PDs $\\{f^{i}\\}$
and $\\{g^{j}\\}$ from the groups approximated by their degree-$\kappa$
estimates:

$\displaystyle f^{i}(p)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\kappa}e^{-\lambda_{k}\sigma}f^{i}_{k}\psi_{k}(p),i=1,\dots,m,$
(22) $\displaystyle g^{j}(p)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\kappa}e^{-\lambda_{k}\sigma}g^{j}_{k}\psi_{k}(p),j=1,\dots,n,$
(23)

where $f^{i}_{k}$ and $g^{j}_{k}$, $k=0,\dots,\kappa$, are the Fourier
coefficients with respect to the $k$-th LB eigenfunction
$\boldsymbol{\psi}_{k}$. Their functional means are

$\displaystyle\bar{f}(p)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\kappa}e^{-\lambda_{k}\sigma}\bar{f}_{k}\psi_{k}(p),$
(24) $\displaystyle\bar{g}(p)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\kappa}e^{-\lambda_{k}\sigma}\bar{g}_{k}\psi_{k}(p),$
(25)

where $\bar{f}_{k}=\frac{1}{m}\sum_{i=1}^{m}f^{i}_{k}$ and
$\bar{g}_{k}=\frac{1}{n}\sum_{j=1}^{n}g^{j}_{k}$ are the mean Fourier
coefficients. We then use the $L_{2}$-norm difference
$\left\|\bar{f}-\bar{g}\right\|^{2}_{2}$ between the functional means as a
test statistic for measuring the group difference in HK expansion of the PDs.
We can algebraically show that

$\left\|\bar{f}-\bar{g}\right\|^{2}_{2}=\sum_{k=0}^{\kappa}e^{-\lambda_{k}\sigma}(\bar{f}_{k}-\bar{g}_{k})^{2}.$
(26)

In a standard approximate permutation test, the subject labels of the two
groups are randomly exchanged. Here, we build the permutation test on
transposition $\pi_{ij}$ that only exchanges the $i$-th and $j$-th subject
labels between $\\{f^{i},i=1,\dots,m\\}$ and $\\{g^{j},j=1,\dots,n\\}$ and
keeps all the other PDs fixed, i.e.

$\displaystyle\pi_{ij}(f^{1},\dots,f^{m})$ $\displaystyle=$
$\displaystyle(f^{1},\dots,g^{j},\dots,f^{m}),$ (27)
$\displaystyle\pi_{ij}(g^{1},\dots,g^{n})$ $\displaystyle=$
$\displaystyle(g^{1},\dots,f^{i},\dots,g^{n}),$ (28)

which we call a spectral transposition. Any permutation of the two groups of
$m$ and $n$ subjects is reachable by a sequence of transpositions, which has
been shown to be computationally much more efficient than the standard
permutation testing procedure of exchanging all labels at once (Chung, Xie .,
2019). We generate the empirical distribution for the permutation test through
the spetral transpositions. In one spectral transposition $\pi_{ij}$, we
obtain the $L_{2}$-distance between the functional means of the
degree-$\kappa$ HK estimates of PDs based on transposed labels:

$L_{2}(f,g)=\left\|\bar{f}^{\prime}-\bar{g}^{\prime}\right\|^{2}_{2}=\sum_{k=0}^{\kappa}e^{-\lambda_{k}\sigma}(\bar{f}^{\prime}_{k}-\bar{g}^{\prime}_{k})^{2},$
(29)

where

$\bar{f}^{\prime}_{k}=\bar{f}_{k}+\frac{1}{m}(g^{j}_{k}-f^{i}_{k})\mbox{ and
}\bar{g}^{\prime}_{k}=\bar{g}_{k}+\frac{1}{n}(f^{i}_{k}-g^{j}_{k})$

are the means of transposed Fourier coefficients. Since we know $\bar{f}_{k}$
and $\bar{g}_{k}$ already, we simply update the terms
$\frac{1}{m}(g^{j}_{k}-f^{i}_{k})$ and $\frac{1}{n}(f^{i}_{k}-g^{j}_{k})$
affected by the transposition. The $p$-value of the spectral permutation test
is then calculated as the proportion of $L_{2}$-distances in the empirical
distribution exceeding the $L_{2}$-distance between the observed PDs. To
ensure convergence, we perform upward of 100,000 permutations until the
$p$-value stabilizes.

###  2.3 Topological analysis of variance via transpositions on HK-estimated
PDs

Topological analysis of variance allows us to assess within- and between-group
similarity and dissimilarity in PDs across multiple groups. The challenge of
applying an ANOVA procedure to raw PDs is that they do not have unique means
(Mileyko ., 2011). Thus, Heo . (2012) applied the standard ANOVA procedure to
raw PDs reduced in dimensionality via Isomap. In contrast, our HK-estimates of
PDs have well-defined functional means and $L_{2}$-distance through Fourier
coefficients, which provides a natural framework for topological analysis of
variance on PDs without any dimensionality reduction beforehand.

To describe our heuristics in constructing an effective topological ANOVA
framework, suppose the $K$ groups of HK-estimated PDs are expressed as
follows:

$\begin{array}[]{ccccc}\text{Group 1}:&f^{11}&f^{12}&\cdots&f^{1n_{1}}\\\
\text{Group 2}:&f^{21}&f^{22}&\cdots&f^{2n_{2}}\\\ \vdots&\vdots\\\
\text{Group K}:&f^{K1}&f^{K2}&\cdots&f^{Kn_{K}}\end{array}$

Motivated by the standard ANOVA procedure, we could try and build an
$F$-statistic comparing $K$ groups of HK-estimated PDs through the
$L_{2}$-distance in (26). A topological between-group sum of squares could
take the form of

$\sum_{i=1}^{K}n_{i}||\bar{f}^{i}-\bar{f}||^{2}_{2},$ (30)

and a topological within-group sum of squares the form of

$\sum_{i=1}^{K}\sum_{j=1}^{n_{i}}||f^{ij}-\bar{f}^{i}||^{2}_{2},$ (31)

where $f^{ij}$ is the HK-estimate of the $j$-th PD of the $i$-th group,
$\bar{f}^{i}$ is the functional mean of the HK-estimates of PDs in the $i$-th
group, and $\bar{f}$ is the grand functional mean over the HK-estimates of all
PDs. The functional means would serve as the topological centroids. Ideally
the $F$-statistic would follow $F$-distribution under some mild normality
assumptions on the HK-estimated PDs, such as

$\frac{\sum_{i=1}^{K}n_{i}||\bar{f}^{i}-\bar{f}||^{2}_{2}/K-1}{\sum_{i=1}^{K}\sum_{j=1}^{n_{i}}||f^{ij}-\bar{f}^{i}||^{2}_{2}/N-K}\sim
F_{K-1,N-K},$

with $N=\sum_{i=1}^{K}n_{i}$ and

$\sum_{i=1}^{K}n_{i}||\bar{f}^{i}-\bar{f}||^{2}_{2}\sim\chi^{2}_{K-1},$
$\sum_{i=1}^{K}\sum_{j=1}^{n_{i}}||f^{ij}-\bar{f}^{i}||^{2}_{2}\sim\chi^{2}_{N-K}.$

However, normality assumptions for heterogeneous features like PDs may be too
strong to satisfy on multivariate data.

Instead of fiddling with parametric constraints, we use a permutational ANOVA
approach that bypasses the distributional issue and has found significant
applications on multivariate data in response to complex experimental designs
of ecological studies, where variables usually consist of counts of counts,
percentage cover, frequencies, or biomass for a large number of species, and
many other fields including chemistry, social sciences, agriculture, medicine,
genetics, psychology, economics (Anderson, 2001, 2017). Here we build our test
statistic for the permutational ANOVA based on pre-calculated pairwise
distances between PDs so that no recalculation of distances is required after
each transposition. We will only need to update the within- and between-group
sums of distances after each transposition. We will refer to our topological
ANOVA procedure as T-ANOVA, where we define the topological between-group sum
of squares (TSSB) and topological within-group sum of squares (TSSW) based on
sums of pairwise $L_{2}$-distances:

TSSB $\displaystyle=$ $\displaystyle\sum_{\begin{subarray}{c}i,i^{\prime}=1\\\
i<i^{\prime}\end{subarray}}^{K}\sum_{j,j^{\prime}}||f^{ij}-f^{i^{\prime}j^{\prime}}||^{2}_{2}$
(32) TSSW $\displaystyle=$
$\displaystyle\sum_{i=1}^{K}\sum_{j<j^{\prime}}||f^{ij}-f^{ij^{\prime}}||^{2}_{2}.$
(33)

We measure the between- and within-group disparity with the ratio statistic

$\phi=\frac{\text{TSSB}}{\text{TSSW}}.$ (34)

In each transposition, we randomly sample the group labels $i_{1}$ and $i_{2}$
out of the $K$ groups with respect to the proportions of the group sizes
$n_{i}/N$. We then uniformly sample the subject labels $j_{1}$ and $j_{2}$ out
of the $i_{1}$-th and $i_{2}$-th group respectively for transposition. We can
prove by induction that any permutation between the groups can be reached by a
sequence of transpositions through Theorem 1 in (Chung, Xie ., 2019) showing
any permutation between two groups can be reached by a sequence of
transpositions.

In a transposition, we only update the pairwise $L_{2}$-distances in TSSB and
TSSW affected by the transposition:

$\displaystyle\text{TSSW}^{\prime}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{K}\sum_{j<j^{\prime}}||f^{ij}-f^{ij^{\prime}}||^{2}_{2}$
(35) $\displaystyle+$ $\displaystyle\sum_{j^{\prime}\neq
j_{2}}||f^{i_{1}j_{1}}-f^{i_{2}j^{\prime}}||^{2}_{2}-\sum_{j^{\prime}\neq
j_{1}}||f^{i_{1}j_{1}}-f^{i_{1}j^{\prime}}||^{2}_{2}$ $\displaystyle+$
$\displaystyle\sum_{j^{\prime}\neq
j_{1}}||f^{i_{2}j_{2}}-f^{i_{1}j^{\prime}}||^{2}_{2}-\sum_{j^{\prime}\neq
j_{2}}||f^{i_{2}j_{2}}-f^{i_{2}j^{\prime}}||^{2}_{2}$ $\displaystyle=$ TSSW
$\displaystyle+$ $\displaystyle\sum_{j^{\prime}\neq
j_{2}}||f^{i_{1}j_{1}}-f^{i_{2}j^{\prime}}||^{2}_{2}-\sum_{j^{\prime}\neq
j_{1}}||f^{i_{1}j_{1}}-f^{i_{1}j^{\prime}}||^{2}_{2}$ $\displaystyle+$
$\displaystyle\sum_{j^{\prime}\neq
j_{1}}||f^{i_{2}j_{2}}-f^{i_{1}j^{\prime}}||^{2}_{2}-\sum_{j^{\prime}\neq
j_{2}}||f^{i_{2}j_{2}}-f^{i_{2}j^{\prime}}||^{2}_{2},$

where we adjust terms involving only groups $i_{1}$ and $i_{2}$ with (35) and
(2.3).

$\displaystyle\text{TSSB}^{\prime}$ $\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}i,i^{\prime}=1\\\
i<i^{\prime}\end{subarray}}^{K}\sum_{j,j^{\prime}}||f^{ij}-f^{i^{\prime}j^{\prime}}||^{2}_{2}$
(37) $\displaystyle-$ $\displaystyle\sum_{j^{\prime}\neq
j_{2}}||f^{i_{1}j_{1}}-f^{i_{2}j^{\prime}}||^{2}_{2}+\sum_{j^{\prime}\neq
j_{1}}||f^{i_{1}j_{1}}-f^{i_{1}j^{\prime}}||^{2}_{2}$ $\displaystyle-$
$\displaystyle\sum_{j^{\prime}\neq
j_{1}}||f^{i_{2}j_{2}}-f^{i_{1}j^{\prime}}||^{2}_{2}+\sum_{j^{\prime}\neq
j_{2}}||f^{i_{2}j_{2}}-f^{i_{2}j^{\prime}}||^{2}_{2}$ (38) $\displaystyle+$
$\displaystyle\mathbb{1}(i^{\prime}\neq
i_{1},i_{2})\sum_{j^{\prime}=1}^{n_{i^{\prime}}}(||f^{i_{2}j_{2}}-f^{i^{\prime}j^{\prime}}||^{2}_{2}-||f^{i_{1}j_{1}}-f^{i^{\prime}j^{\prime}}||^{2}_{2})$
(39) $\displaystyle+$ $\displaystyle\mathbb{1}(i^{\prime}\neq
i_{1},i_{2})\sum_{j^{\prime}=1}^{n_{i^{\prime}}}(||f^{i_{1}j_{1}}-f^{i^{\prime}j^{\prime}}||^{2}_{2}-||f^{i_{2}j_{2}}-f^{i^{\prime}j^{\prime}}||^{2}_{2})$
$\displaystyle=$ TSSB $\displaystyle-$ $\displaystyle\sum_{j^{\prime}\neq
j_{2}}||f^{i_{1}j_{1}}-f^{i_{2}j^{\prime}}||^{2}_{2}+\sum_{j^{\prime}\neq
j_{1}}||f^{i_{1}j_{1}}-f^{i_{1}j^{\prime}}||^{2}_{2}$ $\displaystyle-$
$\displaystyle\sum_{j^{\prime}\neq
j_{1}}||f^{i_{2}j_{2}}-f^{i_{1}j^{\prime}}||^{2}_{2}+\sum_{j^{\prime}\neq
j_{2}}||f^{i_{2}j_{2}}-f^{i_{2}j^{\prime}}||^{2}_{2},$

where we adjust terms involving only groups $i_{1}$ and $i_{2}$ with (37) and
(38), terms involving groups other than $i_{2}$ that are affected by $i_{1}$
with (39), and terms involving groups other than $i_{1}$ that are affected by
$i_{2}$ with (2.3). The ratio statistic is then updated to

$\phi^{\prime}=\frac{\text{TSSB}^{\prime}}{\text{TSSW}^{\prime}}.$ (41)

The $p$-value of the T-ANOVA test is then calculated as the proportion of
$\phi^{\prime}$ in the empirical distribution exceeding the $\phi$ between the
observed PDs. We keep the transposed labels as the current labels on which we
build the next transposition and randomize all labels every 500 transpositions
to improve convergence rate.

##  3 Performance Evaluation

We conduct two sets of simulation studies to evaluate performance of the two-
sample transposition test and T-ANOVA.

###  3.1 Performance of two-sample transposition test

We investigate how the spectral transposition test detects underlying
topological similarity and dissimilarity at the presence of topological noise
and artifact.

Figure 4: Left: We randomly sample 100 points from the image with an innate
shape of a key. Top Right Row: Underlying key shape and possible locations of
topological noise in the form of a small hole. Bottom Right Row: Variants of
the key shape in Group 2. They could appear in the 4 pre-specified forms or
randomly out of the variants.

####  3.1.1 Power of detecting hole in structure

We evaluate the power of the transposition test in detecting a key shape with
a distinct hole (Figure 4 left), under different heterogeneity conditions. In
each simulation, two groups of five 100-point point clouds are generated: the
100 points in each point cloud of the first group are generated randomly from
the part of the rectangular image, whereas the 100 points in each point cloud
of the second group are generated randomly with a varied percentage (90%, 95%,
100%) of points from the shape of the key. Rips filtration is constructed on
each point cloud. The proposed spectral permutation test is then applied to
compare the PDs of the Rips filtrations in the two groups. When there are
respectively 90%, 95%, and 100% points sampled from the shape of the key in
the second group, the spectral permutation test rejects ($p$-value $<$ 0.05)
the null hypothesis of no group difference in 91, 100, and 100% of 100
simulations (corresponding means $\pm$ standard deviations of $p$-values:
0.0124$\pm$0.0327, 0.0041$\pm$0.0125, 0.0008$\pm$0.0057, showing that the test
stays sensitive in detecting the group shape difference when points in the
second group are not entirely sampled from the shape of the key.

####  3.1.2 Robustness of performance under variation of topological noise and
hole location

We conduct two studies to assess the robustness of the test when the
underlying topological structure is ’contaminated’ with heterogeneous
topological noise and when the underlying structure undergoes non-topological
changes.

We first evaluate the robustness of performance under heterogeneity of
topological noise. In each of 100 simulations, we use the spectral
transposition test to compare Group 1 of $m$ random samples with a varied
percentage (90%, 95%, 100%) of 100 points from the original key shape with
Group 2 of $n$ random samples from the key shape ’contaminated’ with
topological noise in the form of a much smaller hole next to the keyhole with
pre-specified (in such case $m=4$ vs $n=4$) or random locations (in such case
$m=5$ vs $n=5$, $m=20$ vs $n=20$, or $m=100$ vs $n=100$). Figure 4 (top right
row) shows the 4 possible locations of the topological noise in Group 2. We
expect the test to stay robust to this topological noise. Table 1 summarizes
the results for different percentage of points when the topological noise
appears at pre-specified vs random locations. The spectral transposition test
stays robust to the topological noise in fixed and random locations.

Robustness under Variation of Topological Noise
---
Percentage | Pre-specified Loc. ($m=4$ vs $n=4$) | Random Loc. ($m=5$ vs $n=5$)
100% | 0.4567$\pm$0.2874 | 0.4133$\pm$0.2482
95% | 0.4777$\pm$0.2843 | 0.4498$\pm$0.2844
90% | 0.4455$\pm$0.2791 | 0.5214$\pm$0.2983
Percentage | Random Loc. ($m=20$ vs $n=20$) | Random Loc. ($m=100$ vs $n=100$)
100% | 0.5060$\pm$0.3163 | 0.4328$\pm$0.2764
95% | 0.5016$\pm$0.2998 | 0.4193$\pm$0.2863
90% | 0.4827$\pm$0.2919 | 0.5260$\pm$0.2812
Robustness under Variation of Hole Location
Percentage | Pre-specified Loc. ($m=4$ vs $n=4$) | Random Loc. ($m=5$ vs $n=5$)
100% | 0.2917$\pm$0.2624 | 0.5005$\pm$0.2883
95% | 0.2973$\pm$0.2407 | 0.5342$\pm$0.2775
90% | 0.3065$\pm$0.2505 | 0.4434$\pm$0.3050
Percentage | Radom Loc. ($m=20$ vs $n=20$) | Random Loc. ($m=100$ vs $n=100$)
100% | 0.4998$\pm$0.2901 | 0.3845$\pm$0.2620
95% | 0.4608$\pm$0.2999 | 0.3924$\pm$0.2777
90% | 0.4810$\pm$0.2568 | 0.4550$\pm$0.2835
Table 1: Summary of mean$\pm$standard deviation of $p$-values from the
spectral transposition test in 100 simulations. Top half: In each simulation,
the test is used to compare a group of $m$ random samples with a varied
percentage (90%, 95%, 100%) of 100 points from the original key shape with a
group of $n$ random samples with the same percentage of 100 points from the
key shape with topological noise in the form of a much smaller hole next to
the keyhole. The location of the smaller hole in each random sample of the
second group can be pre-specified or randomly chosen from the pre-specified
options. Bottom half: In each simulation, the test is used to compare a group
of $m$ random samples with a varied percentage (90%, 95%, 100%) of 100 points
from the original key shape with only the top left quarter of the keyhole
left, with a group of $n$ random samples with the same percentage of 100
points from the original key shape with a random quarter of the keyhole left.

We then evaluate the robustness of performance under variation of hole
location. In each of 100 simulations, we use the spectral transposition test
to compare Group 1 of $m$ random samples with a varied percentage (90%, 95%,
100%) of 100 points from the original key shape with only the top left quarter
of the keyhole left, with Group 2 of $n$ random samples with the same
percentage of 100 points from the original key shape with a pre-specified (in
such case $m=4$ vs $n=4$) or random (in such case $m=5$ vs $n=5$, $m=20$ vs
$n=20$, or $m=100$ vs $n=100$) quarter of the keyhole left. Figure 4 (bottom
right half) shows the 4 possible variants of the keyhole in Group 2. We expect
the test to stay robust to this change in structure, which is not topological
in nature. Table 1 summarizes the results for different percentage of points
when the variants appears at pre-specified vs random locations. The spectral
transposition test stays robust to the structural variants in fixed and random
locations.

####  3.1.3 Computational time

The computational time of the spectral transposition test grows steadily as
the group sample sizes grow. The mean time for each simulation run for $m=4,5$
vs $n=4,5$ between 7 and 10 seconds and standard deviation within 3 seconds.
For $m=20$ vs $n=20$, the mean time for each simulation run is between 8 and
10 seconds and standard deviation within 3 seconds. For $m=100$ vs $n=100$,
the mean time for each simulation run is between 9 and 11 seconds and standard
deviation within 3 seconds.

###  3.2 Performance of T-ANOVA

In each of the simulation studies in this section, we test the performance of
the T-ANOVA in comparing three groups of point clouds simulated under
different settings. The performance is compared against the standard PERMANOVA
test (Anderson, 2001), as well as the topological analysis of variance test
proposed by Heo . (2012) that runs the univariate ANOVA on dimensionality-
reduced PDs by Isomap.

Figure 5: An example of $n_{1}=n_{2}=n_{3}=5$ 100-point point clouds where the
100 points in each point cloud of the first two groups are generated randomly
from the part of the rectangular image, whereas the 100 points in each point
cloud of the third group are generated randomly with 95% of points from the
shape of the key.

####  3.2.1 Sensitivity in detecting differential hole presence among multiple
groups

In each simulation, three groups of $n_{1},n_{2},n_{3}$ 100-point point clouds
are generated, where the 100 points in each point cloud of the first two
groups are generated randomly from the part of the rectangular image, whereas
the 100 points in each point cloud of the third group are generated randomly
with a varied percentage (90%, 95%, 100%) of points from the shape of the key
(Figure 5). Table 2 shows the results of the T-ANOVA test in comparison with
the other tests.

Sensitivity in Detecting Differential Hole Presence Among Multiple Groups
---
$n_{1}=n_{2}=n_{3}=5$
Percentage of Points in Key Shape | T-ANOVA | Heo’s ANOVA | PERMANOVA
100% | $0.0009\pm 0.0016$ | $0.0038\pm 0.0101$ | $0.0025\pm 0.0033$
95% | $0.0016\pm 0.0033$ | $0.0109\pm 0.0258$ | $0.0045\pm 0.0064$
90% | $0.0060\pm 0.0146$ | $0.0436\pm 0.0791$ | $0.0192\pm 0.0438$
$n_{1}=n_{2}=n_{3}=20$
Percentage of Points in Key Shape | T-ANOVA | Heo’s ANOVA | PERMANOVA
100% | $0.0000\pm 0.0000$ | $0.0000\pm 0.0000$ | $0.0000\pm 0.0000$
95% | $0.0000\pm 0.0000$ | $0.0000\pm 0.0000$ | $0.0000\pm 0.0000$
90% | $0.0000\pm 0.0000$ | $0.0000\pm 0.0001$ | $0.0000\pm 0.0000$
$n_{1}=5,n_{2}=20,n_{3}=100$
Percentage of Points in Key Shape | T-ANOVA | Heo’s ANOVA | PERMANOVA
100% | $0.0253\pm 0.0699$ | $0.0353\pm 0.0101$ | $0.0011\pm 0.0055$
95% | $0.0485\pm 0.0827$ | $0.0511\pm 0.1255$ | $0.0028\pm 0.0117$
90% | $0.0998\pm 0.1377$ | $0.1281\pm 0.2395$ | $0.0091\pm 0.0238$
Table 2:  Summary of mean$\pm$standard deviation of $p$-values from the
T-ANOVA, Heo’s ANOVA, and PERMANOVA in 100 simulations. In each simulation,
three groups of $n_{1},n_{2},n_{3}$ 100-point point clouds are generated,
where the 100 points in each point cloud of the first two groups are generated
randomly from the part of the rectangular image, whereas the 100 points in
each point cloud of the third group are generated randomly with a varied
percentage (90%, 95%, 100%) of points from the shape of the key. Robustness
under Variation of Topological Noise
---
Pre-specified Location: $n_{1}=n_{2}=n_{3}=4$
Percentage of Points in Key Shape | T-ANOVA | Heo’s ANOVA | PERMANOVA
100% | $0.4860\pm 0.2607$ | $0.5267\pm 0.2675$ | $0.4901\pm 0.2965$
95% | $0.4897\pm 0.2608$ | $0.5211\pm 0.2822$ | $0.5050\pm 0.2851$
90% | $0.4974\pm 0.2963$ | $0.5163\pm 0.3125$ | $0.4619\pm 0.2667$
Random Location: $n_{1}=n_{2}=n_{3}=5$
Percentage of Points in Key Shape | T-ANOVA | Heo’s ANOVA | PERMANOVA
100% | $0.5275\pm 0.3052$ | $0.4988\pm 0.2882$ | $0.5518\pm 0.2952$
95% | $0.5166\pm 0.3138$ | $0.5062\pm 0.2830$ | $0.5308\pm 0.2870$
90% | $0.4970\pm 0.2975$ | $0.5205\pm 0.2689$ | $0.4923\pm 0.2928$
Random Location: $n_{1}=n_{2}=n_{3}=20$
Percentage of Points in Key Shape | T-ANOVA | Heo’s ANOVA | PERMANOVA
100% | $0.4915\pm 0.2624$ | $0.4812\pm 0.2640$ | $0.4727\pm 0.2720$
95% | $0.4860\pm 0.2862$ | $0.4864\pm 0.3072$ | $0.5141\pm 0.2734$
90% | $0.4349\pm 0.2759$ | $0.5139\pm 0.2791$ | $0.4706\pm 0.2996$
Random Location: $n_{1}=5,n_{2}=20,n_{3}=100$
Percentage of Points in Key Shape | T-ANOVA | Heo’s ANOVA | PERMANOVA
100% | $0.5374\pm 0.2937$ | $0.5001\pm 0.3028$ | $0.4593\pm 0.2836$
95% | $0.5413\pm 0.2747$ | $0.5016\pm 0.2738$ | $0.4692\pm 0.3037$
90% | $0.4938\pm 0.2855$ | $0.5014\pm 0.2955$ | $0.5085\pm 0.3040$
Robustness under Variation of Hole Location
Pre-specified Location: $n_{1}=n_{2}=n_{3}=4$
Percentage of Points in Key Shape | T-ANOVA | Heo’s ANOVA | PERMANOVA
100% | $0.4887\pm 0.2873$ | $0.4664\pm 0.2537$ | $0.5057\pm 0.2862$
95% | $0.4948\pm 0.2776$ | $0.4479\pm 0.2889$ | $0.5109\pm 0.2779$
90% | $0.4364\pm 0.2572$ | $0.4651\pm 0.2916$ | $0.4463\pm 0.2767$
Random Location: $n_{1}=n_{2}=n_{3}=5$
Percentage of Points in Key Shape | T-ANOVA | Heo’s ANOVA | PERMANOVA
100% | $0.4785\pm 0.2737$ | $0.5505\pm 0.2706$ | $0.4932\pm 0.3036$
95% | $0.5039\pm 0.3082$ | $0.4872\pm 0.3055$ | $0.4572\pm 0.3079$
90% | $0.4887\pm 0.2919$ | $0.5302\pm 0.2784$ | $0.4095\pm 0.2784$
Random Location: $n_{1}=n_{2}=n_{3}=20$
Percentage of Points in Key Shape | T-ANOVA | Heo’s ANOVA | PERMANOVA
100% | $0.5249\pm 0.3118$ | $0.4839\pm 0.2928$ | $0.5048\pm 0.3159$
95% | $0.5292\pm 0.3082$ | $0.4528\pm 0.2790$ | $0.5398\pm 0.3044$
90% | $0.5183\pm 0.3009$ | $0.5419\pm 0.2801$ | $0.5198\pm 0.2870$
Random Location: $n_{1}=5,n_{2}=20,n_{3}=100$
Percentage of Points in Key Shape | T-ANOVA | Heo’s ANOVA | PERMANOVA
100% | $0.4411\pm 0.2541$ | $0.4738\pm 0.3197$ | $0.5481\pm 0.2707$
95% | $0.4498\pm 0.2774$ | $0.4623\pm 0.2825$ | $0.5271\pm 0.2884$
90% | $0.4670\pm 0.2849$ | $0.4953\pm 0.2868$ | $0.4747\pm 0.2870$
Table 3:  Summary of mean$\pm$standard deviation of $p$-values of the T-ANOVA,
Heo’s ANOVA, and PERMANOVA in 100 simulations. Top half: In each simulation,
the test is used to compare Group 1, 2, 3 of respective $n_{1},n_{2},n_{3}$
random samples are generated with a pre-specified percentage (90%, 95%, 100%)
of 100 points from the original key shape ’contaminated’ with topological
noise in the form of a much smaller hole next to the keyhole with pre-
specified (in such case $n_{1}=n_{2}=n_{3}=4$) or random locations (in such
case $n_{1}=n_{2}=n_{3}=5$, $n_{1}=n_{2}=n_{3}=20$, or
$n_{1}=5,n_{2}=20,n_{3}=100$). Bottom half: In each simulation, the test is
used to compare Group 1, 2, 3 of respective $n_{1},n_{2},n_{3}$ random samples
are generated with a pre-specified percentage (90%, 95%, 100%) of 100 points
from the original key shape with a pre-specified (in such case
$n_{1}=n_{2}=n_{3}=4$) or random (in such case $n_{1}=n_{2}=n_{3}=5$,
$n_{1}=n_{2}=n_{3}=20$, or $n_{1}=5,n_{2}=20,n_{3}=100$) quarter of the
keyhole left

####  3.2.2 Robustness under variation of noise and hole location

We conduct two studies to assess the robustness of the test when the
underlying topological structure is ’contaminated’ with heterogeneous
topological noise and when the hole location shifts as in Section 3.1.2.

We first evaluate the robustness of T-ANOVA under heterogeneity of topological
noise. In each of 100 simulations, Group 1, 2, 3 of respective
$n_{1},n_{2},n_{3}$ random samples are generated with a pre-specified
percentage (90%, 95%, 100%) of 100 points from the original key shape
’contaminated’ with topological noise in the form of a much smaller hole next
to the keyhole with pre-specified (in such case $n_{1}=n_{2}=n_{3}=4$) or
random locations (in such case $n_{1}=n_{2}=n_{3}=5$, $n_{1}=n_{2}=n_{3}=20$,
or $n_{1}=5,n_{2}=20,n_{3}=100$). The 4 possible locations of the topological
noise in each group are the same as Figure 4. We use the T-ANOVA test to
compare PDs across the three groups. We expect the test to stay robust to the
topological noise. Table 3 (top half) shows the results of the T-ANOVA test in
comparison with standard PERMANOVA and the topological ANOVA proposed by Heo .
(2012).

We then evaluate the robustness of T-ANOVA under variation of hole location.
In each of 100 simulations, Group 1, 2, 3 of respective $n_{1},n_{2},n_{3}$
random samples are generated with a pre-specified percentage (90%, 95%, 100%)
of 100 points from the original key shape with a pre-specified (in such case
$n_{1}=n_{2}=n_{3}=4$) or random (in such case $n_{1}=n_{2}=n_{3}=5$,
$n_{1}=n_{2}=n_{3}=20$, or $n_{1}=5,n_{2}=20,n_{3}=100$) quarter of the
keyhole left. The 4 possible variants of the key shape in each group are the
same as Figure 4 . We use the T-ANOVA test to compare PDs across the three
groups. We expect the test to stay robust to this change in structure, which
is not topological in nature. Table 3 (bottom half) shows the results of the
T-ANOVA test in comparison with the other tests.

Pre-specified Location: $n_{1}=n_{2}=n_{3}=4$
---
Percentage of Points in Key Shape | T-ANOVA | PERMANOVA
100% | $113.92\pm 1.58$ | $8.36\pm 0.10$
95% | $114.06\pm 1.53$ | $8.39\pm 0.08$
90% | $114.34\pm 1.82$ | $8.34\pm 0.06$
Random Location: $n_{1}=n_{2}=n_{3}=5$
Percentage of Points in Key Shape | T-ANOVA | PERMANOVA
100% | $116.46\pm 1.39$ | $9.02\pm 0.07$
95% | $115.82\pm 1.57$ | $9.04\pm 0.13$
90% | $115.91\pm 1.45$ | $9.16\pm 0.26$
Random Location: $n_{1}=n_{2}=n_{3}=20$
Percentage of Points in Key Shape | T-ANOVA | PERMANOVA
100% | $127.74\pm 1.49$ | $36.82\pm 0.34$
95% | $128.07\pm 1.68$ | $36.80\pm 0.22$
90% | $128.66\pm 2.66$ | $36.76\pm 0.26$
Random Location: $n_{1}=5,n_{2}=20,n_{3}=100$
Percentage of Points in Key Shape | T-ANOVA | PERMANOVA
100% | $159.93\pm 2.59$ | $132.14\pm 0.58$
95% | $159.65\pm 2.74$ | $132.54\pm 0.82$
90% | $160.45\pm 2.77$ | $133.32\pm 2.56$
Table 4:  Summary of mean$\pm$standard deviation of computational time in
seconds for a million transpositions by the T-ANOVA and PERMANOVA in 100
simulations under the topological noise setting.

####  3.2.3 Computational time

Table 4 shows the means and standard deviations of computational times for one
million transpositions under the topological noise setting (the hole location
and sensitivity studies have similar computational times, so we only present
one setting here). Just like the two-sample test, T-ANOVA shows steady growth
of computational time as group sample sizes increase, in comparison with the
sharp time growth of PERMANOVA. Heo’s ANOVA is fast as it runs a univariate
ANOVA on the dimensionality-reduced PDs.

###  3.3 Summary

The results show that the performance of our T-ANOVA test is comparable with
the two baseline methods in terms of robustness under variation of topological
noise and hole location, as well as sensitivity in detecting differential hole
presence among multiple groups. In comparison with PERMANOVA, the advantage of
the transposition approach of T-ANOVA shows up in the steady growth of
computational time as group sample sizes increase. Although T-ANOVA is
comparable in performance as Heo’s ANOVA, it does not require dimensionality
reduction of PDs. More importantly, it has a natural framework for distance-
based clustering, which we illustrate in the Application section.

##  4 Application

Stroke is the leading cause of severe adult disability in the United States
(Tsao ., 2022). A left-hemisphere stroke commonly leads to aphasia, a speech-
language disorder often classified into subtypes according to behavioral
symptoms. Traditional subtypes of aphasia are determined through the Aphasia
Quotient (AQ) subtest scores of the Revised Western Aphasia Battery (WAB-R)
(Kertesz, 2007) that assess speech and language abilities such as spontaneous
speech fluency, auditory comprehension, repetition, and naming performance.
These scores binarize the patients into categories. For instance, the
spontaneous speech fluency score ($\geq 5$ vs. $\leq 4$) is a rating based on
subjective evaluation mostly about quantity and grammaticality of output along
with other features, such as word-finding difficulty, paraphasias, and
hesitations. It separates individuals into fluent and non-fluent categories.
Eight traditional subtypes thus arise from the binarized categories of
fluency, comprehension, and repetition (Figure 6).

Figure 6: Traditional aphasia subtypes according to binary categories of
fluency, comprehension, and repetition.

Studies over the years have addressed WAB subtyping issues since its initial
version in 1982 and proposed new ways identifying coherent clusters of aphasia
subtypes (Ferro Kertesz, 1987; Crary ., 1992; John ., 2017; Fromm ., 2022).
Unsupervised learning approaches such as K-means clustering has been applied
to behavioral scores beyond WAB-R to redefine aphasia subtypes (Fromm .,
2022). There is, however, a lack of exploration on aphasia subtyping via
clustering brain network features. Our goal is aimed at identifying patterns
of damage in the brain networks that lead to overlapping behavioral deficits.
This study takes a topological angle at the clustering and inference of the
resting-state functional brain networks of aphasic individuals, and
summarizing basic statistical characteristics of the WAB-R AQ subtest scores
of the clusters.

###  4.1 Data acquisition and preprocessing

The rs-fMRI data were acquired from 103 participants with aphasia resulting
from a single ischemic or hemorrhagic stroke involving the left hemisphere on
a Siemens Prisma 3T scanner with a 20-channel head coil located at the Center
for the Study of Aphasia Recovery at the University of South Carolina. The
following imaging parameters of images were used: a multiband sequence (x2)
with a $216\times 216$ mm field of view, a $90\times 90$ matrix size, and a
72-degree flip angle, 50 axial slices (2 mm thick with $20\%$ gap yielding 2.4
mm between slice centers), repetition time TR =1650 ms, TE=35 ms, GRAPPA=2, 44
reference lines, interleaved ascending slice order. During the scanning
process, the participants were instructed to stay still with eyes closed. A
total of 370 volumes were acquired.

The preprocessing procedures of the rs-fMRI data include motion correction,
brain extraction and time correction using a novel method developed for stroke
patients (Yourganov ., 2018). The Realign and Unwarp procedure in SPM12 with
default settings was used for motion correction. Brain extraction was then
performed using the SPM12 script pm_brain_mask with default settings. Slice
time correction was also done using SPM12. The mean fMRI volume for each
participant was then aligned to the corresponding T2-weighted image to compute
the spatial transformation between the data and the lesion mask. The fMRI data
were then spatially smoothed with a Gaussian kernel with FWHM= 6 mm. To
eliminate artifacts driven by lesions, a pipeline proposed by Yourganov .
(2018) was applied on the the rs-fMRI. The FSL MELODIC package was used to
decompose the data into independent components (ICs) and to compute the
Z-scored spatial maps for the ICs. The spatial maps were thresholded at
$p<0.05$ and compared with the lesion mask for the participant. The Jaccard
index, computed as the ratio between the numbers of voxels in the intersection
and union, was used to quantify the amount of spatial overlap between the
lesion mask and thresholded IC maps, both of which were binary. ICs
corresponding to Jaccard index greater than $5\%$ were deemed significantly
overlapping with the lesion mask and then regressed out of the fMRI data using
the fsl_regfilt script from the FSL package. By applying the automated
anatomical labelling (AAL) atlas, 116 regions of interest (ROIs) were created
and used as nodes in the brain networks subsequently constructed.

The Aphasia Quotient, a score strongly related to the overall lesion damage in
brain, was measured in the participants. In terms of behavioral measures, the
following WAB-R subscores were used to measure performance of participants in
fluency, comprehension, repetition, object naming, and sentence completion:
Information Content, Fluency Rating, Spontaneous Speech Rating, Comprehension
Yes/No Questions, Comprehension Auditory Words, Comprehension Sequential
Commands, Comprehension Subscore, Repetition Subscore, Object Naming, Word
Fluency, Sentence Completion, Responsive Speech, and Naming Subscore.

###  4.2 Resting-state functional brain network and filtration

We first constructed resting-state functional brain networks from the rs-fMRI
described above. The 116 AAL ROIs served as the nodes of the resting-state
functional network of each individual and Pearson’s correlation between the
BOLD signals at two ROIs serve as their edge weight. A Rips filtration was
built on the resting-state functional correlation matrix of each individual.
The PDs decoding the birth and death times of 1-cycles in the individual Rip
filtrations were then smoothed with the HK representation.

###  4.3 Aphasia subtyping via topological clustering of brain networks

Topological clustering has been applied in different angles to studies of
resting-state functional brain networks (Stolz ., 2017; Chung ., 2023). To the
best of our knowledge, this is the first study to explore aphasia subtyping
through topological clustering of resting-state functional brain networks.
Here we take advantage of the HK representation of PDs and extend the T-ANOVA
into a topological clustering scheme where clusters were identified with
respect to topological centroids calculated as the functional means of the HK
estimates of the PDs representing the brain networks of individuals in the
study. We compared the statistical characteristics of the topological clusters
to baseline clusters obtained through K-means clustering of the WAB-R
subscores. We repeated the clustering process 100 times in each instance and
checked for consistency across repetitions. Three topological clusters had the
overall best fit so we compared the results of three baseline clusters with
them.

Figure 7: The lesion map (left two columns) and average absolute connectivity
(right two columns) of three topological and baseline clusters. Figure 8: Box
plots of subscores for all participants and those in each of the three
clusters.

The overall lesion map and average absolute connectivity of three baseline and
topological clusters are shown in Figure 7. The lesion map was created by
augmenting stroke lesion damage in the brain of all subjects within each
cluster. Note that the three baseline clusters appear to be confounded by the
overall lesion extent of the subjects as they show distinctly different lesion
extent (Cluster 1 ¿ Cluster 2 ¿ Cluster 3). This is confirmed by the AQ and
subscore distributions summarized in Figure 8, where the AQ score is known to
positively correlate with lesion extent and the subscore distributions show a
distinct monotone pattern consistent with that of AQ across clusters. On the
other hand, the topological clusters do not appear to be confounded by lesion
extent as the lesion extent do not vary significantly across the clusters and
the subscore distributions do not follow a specific trend with reference to
the AQ score. As of the average connectivity, we see different connectivity
patterns in the three topological clusters, whereas the baseline clusters show
similar connectivity pattens. To confirm that the topological clusters did
capture significant statistical difference in brain networks, we also compared
the brain networks across different clusters through the T-ANOVA on their HK-
estimated PDs. Figure 9 shows the empirical distribution of the ratio
statistic based on $L_{2}$-distances of HK-estimated 1-dimensional PDs within
and between the three clusters over 1 million transpositions. The observed
value of the ratio statistic was 5.4728, yielding a $p$-value of 0 and the
conclusion of significant topological difference between the one-dimensional
hole presence in the three clusters of brain networks.

Figure 9: Empirical distribution of the ratio statistic based on $L_{2}$-distances of HK-estimated 1-dimensional PDs within and between the three clusters over 1 million transpositions. WAB-R Subscore(s) | Cluster 1 | Cluster 2 | Cluster 3
---|---|---|---
Comprehension Yes/No Questions | Low | Medium | Medium
Comprehension Auditory Words | Low | Medium | High
Comprehension Sequential Commands | Low | Medium | High
Comprehension Subscore | Low | Medium | High
Repetition Subscore | Low | Medium | High
Fluency Rating | Low | High | Medium
Word Fluency | Low | Medium | High
Information Content | Low | Medium | Medium
Spontaneous Speech Rating | Low | High | Medium
Naming Subscore | Low | Medium | High
Object Naming | Low | Medium | High
Sentence Completion | Low | Medium | Medium
Responsive Speech | Low | Medium | Medium
Table 5: Pattern of median and interquartile range of WAB-R subscores across
three topological clusters/subtypes.

Now, using the three topological clusters as a basis for exploring aphasia
subtypes, behavioral measures in the form of WAB-R subscores across the three
clusters/subtypes have pattern of median and interquartile range summarized in
Table 5. In terms of the three categories (fluency, comprehension, repetition)
used for traditional aphasia subtyping, the comprehension subscores
(Comprehension Yes/No Questions, Comprehension Auditory Words, Comprehension
Sequential Commands, Comprehension Subscore), Repetition Subscore, and Word
Fluency show an overall pattern low-medium-high in medians across Cluster 1, 2
& 3, with the exception of Comp. Yes/No Qs which sees some leveling off in
Cluster 2 & 3, whereas Fluency Rating shows a low-high-medium pattern across
the three clusters. The low, medium, and high are all in comparison to the
median over all subjects.

##  5 Discussion

In this study, we established a topological inference framework based on HK
representation of PDs. Although it does not require the PDs to be extracted
from a specific type of data, we centered the application of the methods
around group comparison of PDs from brain networks. But simulating brain
networks with holes is not straightforward, so we used point clouds from
images with an underlying shape in illustration and simulation studies. We
also extended the framework to topological clustering of brain networks with
application to subtyping individuals with post-stroke aphasia.

Methodologically, the topological inference framework filled a few gaps left
from our previous works. As we pointed out in Section 2, the spectral
transposition test generalizes the permutation test proposed by Wang . (2018)
that compares single-trial signals by permuting coefficients respective of
Fourier basis functions. We can now permute the Fourier coefficients of the HK
estimates of two groups of PDs, which may come from multi-trial univariate or
multivariate signals. Thus the framework is now not only applicable for
single-trial univariate signals, but also for multi-trial univariate and
multivariate signals. Furthermore, we now have T-ANOVA that can compare the
topological features of multi-group univariate and multivariate signals
without further reducing the dimensionality of the features. The multi-group
transposition approach can also be used for speeding up ANOVA procedures in
non-topological settings. Resampling is also of high relevance to deep
learning. Since the power of deep learning is constrained in small sample
schemes, data augmentation methods are needed to increase the training data by
resampling (SG. Huang ., 2021). In future studies, the proposed spectral
permutation method can be easily adapted for deep learning where the input is
augmented persistence features reconstructed from resampled HK coefficients of
PDs.

Since the analytical paradigm proposed in the methods section was already
complicated, we featured topological clustering only in application. In the
application of topological clustering and inference to subtyping individuals
with post-stroke aphasia, one would argue that three clusters may be too
sparse for actual clinical interpretation even though three clusters were
empirically determined to have the best fit. Future studies can refine the
approach by exploring more clusters, e.g. matching the number of traditional
subtypes to see if they have any consistency.

## Acknowledgments

The authors would like to thank Dr. Roger Newman-Norlund for facilitating
access to the fMRI dataset used in this study and Dr. Moo Chung for helpful
discussions on the early versions of the manuscript. Funding sources: NIH
R01DC017162 and R01DC01716202S1 (PI: RHD). Author contributions:
conceptualization (YW, RHD), statistical analysis (YW, JY), interpretation of
results (YW, RHD), writing and editing (All). Compliance with ethical
standards: The neuroimaging scans were approved by the Institutional Review
Board (IRB) at the University of South Carolina. Conflict of interest and
disclosure: None.

## References

  * Adams . (2017) Adams2017Adams, H., Emerson, T., Kirby, M., Neville, R., Peterson, C., Shipman, P.Ziegelmeier, L.  2017\.  Persistence images: A stable vector representation of persistent homology Persistence images: A stable vector representation of persistent homology.  The Journal of Machine Learning Research18218–252.
  * Aldous (1983) Aldous1983Aldous, D.  1983\.  Random walks on finite groups and rapidly mixing Markov chains Random walks on finite groups and rapidly mixing Markov chains.  Séminaire de Probabilités XVII 1981/82 Séminaire de probabilités xvii 1981/82 ( 243–297).  Springer.
  * Aldous Diaconis (1986) Aldous1986Aldous, D. Diaconis, P.  1986\.  Shuffling cards and stopping times Shuffling cards and stopping times.  The American Mathematical Monthly93333–348.
  * Anderson (2001) Anderson2001Anderson, M.  2001\.  A new method for non-parametric multivariate analysis of variance A new method for non-parametric multivariate analysis of variance.  Austral Ecology2632 - 46.
  * Anderson (2017) Anderson2017Anderson, M.  2017\.  Permutational multivariate analysis of variance (PERMANOVA) Permutational multivariate analysis of variance (permanova).  Wiley StatsRef : Statistics Reference Online.
  * Ashburner Friston (2000) Ashburner2000Ashburner, J. Friston, K.  2000\.  Voxel-based morphometry - the methods Voxel-based morphometry - the methods.  NeuroImage11805–821.
  * Bassett (2006) Bassett2006Bassett, D.  2006\.  Adaptive reconfiguration of fractal small-world human brain functional networks Adaptive reconfiguration of fractal small-world human brain functional networks.  Proceedings of the National Academy of Sciences of the United States of America20319518–23.
  * Betzel Bassett (2017) Betzel2017Betzel, R. Bassett, D.  2017\.  Multi-scale brain networks Multi-scale brain networks.  NeuroImage16073 - 83.
  * Bien Tibshirani (2011) Bien2011Bien, J. Tibshirani, R.  2011\.  Hierarchical clustering with prototypes via minimax linkage Hierarchical clustering with prototypes via minimax linkage.  Journal of the American Statistical Association1061075 - 1084.
  * Bubenik (2015) Bubenik2015Bubenik, P.  2015\.  Statistical topological data analysis using persistence landscapes Statistical topological data analysis using persistence landscapes.  Journal of Machine Learning Research1677 - 102.
  * Carlsson (2009) Carlsson2009Carlsson, G.  2009\.  Topology and data Topology and data.  Bulletin of the American Mathematical Society46255 – 308.
  * Carriére . (2015) Carriere2015Carriére, M., Oudot, S. Ovsjanikov, M.  2015\.  Stable topological signatures for points on 3D shapes Stable topological signatures for points on 3d shapes.  Eurographics Symposium on Geometry Processing34.
  * Chen . (2015) Chen2015Chen, YC., Wang, D., Rinaldo, A. Wasserman, L.  2015\.  Statistical analysis of persistence intensity functions Statistical analysis of persistence intensity functions.  arXiv:1510.02502v1.
  * Chung . (2009) Chung2009Chung, M., Bubenik, P. Kim, P.  2009\.  Persistence diagrams of cortical surface data Persistence diagrams of cortical surface data.  Proceedings of Information Processing in Medical Imaging (IPMI)386 - 397.
  * Chung . (2007) Chung2007Chung, M., Dalton, K., Shen, L., Evans, A. Davidson, R.  2007\.  Weighted Fourier representation and its application to quantifying the amount of gray matter Weighted Fourier representation and its application to quantifying the amount of gray matter.  IEEE Transactions on Medical Imaging26566–581.
  * Chung, Lee . (2019) Chung2019aChung, M., Lee, H., Ombao, H. Solo, V.  2019\.  Exact topological inference of the resting-state brain networks in twins Exact topological inference of the resting-state brain networks in twins.  Network Neuroscience.
  * Chung . (2023) Chung2023Chung, M., Ramos, C., de Paiva, F., Mathis, J., Prabharakaren, V., Nair, V.Struck, A.  2023\.  Unified topological inference for brain networks in temporal lobe epilepsy using the Wasserstein distance Unified topological inference for brain networks in temporal lobe epilepsy using the wasserstein distance.  arXiv: 2302.06673.
  * Chung . (2014) Chung2014.MICCAIChung, M., Schaefer, S., van Reekum, CM., Schmitz, L., Sutterer, M. Davidson, R.  2014\.  A unified kernel regression on manifolds detects aging-related changes in the amygdala and hippocampus A unified kernel regression on manifolds detects aging-related changes in the amygdala and hippocampus.  MICCAI, Lecture Notes in Computer Science (LNCS)8674789-796.
  * Chung . (2018) Chung2018.EMBCChung, M., Wang, Y. Wu, G.  2018\.  Heat kernel smoothing in irregular image domains Heat kernel smoothing in irregular image domains.  International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)5101-5104.
  * Chung, Xie . (2019) Chung2019bChung, M., Xie, L., Huang, SG., Wang, Y., Yan, J. Shen, L.  2019\.  Rapid acceleration of the permutation test via transpositions Rapid acceleration of the permutation test via transpositions.  International Workshop on Connectomics in NeuroImaging, Lecture Notes in Computer Science1184842-53.
  * Crary . (1992) Crary1992Crary, MA., Wertz, RT. Deal, JL.  1992\.  Classifying aphasias: Cluster analysis of Western aphasia battery and Boston diagnostic aphasia examination results Classifying aphasias: Cluster analysis of western aphasia battery and boston diagnostic aphasia examination results.  Aphasiology629–36.
  * Drakesmith . (2015) Drakesmith2015Drakesmith, M., Caeyenberghs, K., Dutt, A., Lewis, G., David, A. Jones, D.  2015\.  Overcoming the effects of false positives and threshold bias in graph theoretical analyses of neuroimaging data Overcoming the effects of false positives and threshold bias in graph theoretical analyses of neuroimaging data.  NeuroImage118313–333.
  * Edelsbrunner . (2002) Edelsbrunner2002Edelsbrunner, H., Letscher, D. Zomorodian, A.  2002\.  Topological persistence and simplification Topological persistence and simplification.  Discrete Computational Geometry511 - 533.
  * Fasy . (2014) Fasy2014Fasy, B., Lecci, F., Rinaldo, A., Wasserman, L., Balakrishnan, S. Singh, A.  2014\.  Confidence sets for persistence diagrams Confidence sets for persistence diagrams.  The Annals of Statistics422301 - 2339.
  * Ferro Kertesz (1987) Ferro1987Ferro, JM. Kertesz, A.  1987\.  Comparative classification of aphasic disorders Comparative classification of aphasic disorders.  Journal of clinical and experimental Neuropsychology9365–375.
  * Fromm . (2022) Fromm2022Fromm, D., Greenhouse, J., Pudil, M., Shi, Y. MacWhinney, B.  2022\.  Enhancing the classification of aphasia: a statistical analysis using connected speech Enhancing the classification of aphasia: a statistical analysis using connected speech.  Aphasiology36Fromm, Davida Greenhouse, Joel Pudil, Mitchell Shi, Yichun MacWhinney, Brian.
  * Garrison . (2015) Garrison2015Garrison, K., Scheinost, D., Finn, E., Shen, X. Constable, RT.  2015\.  The (in)stability of functional brain network measures across thresholds The (in)stability of functional brain network measures across thresholds.  NeuroImage118651 - 661.
  * Hatcher (2001) Hatcher2001Hatcher, A.  2001\.  Algebraic Topology Algebraic topology.  Cambridge University Press.
  * Heo . (2012) Heo2012Heo, G., Gamble, J. Kim, P.  2012\.  Topological analysis of variance and the maxillary complex Topological analysis of variance and the maxillary complex.  Journal of the American Statistical Association107477 - 492.
  * S. Huang . (2010) Huang2010Huang, S., Li, J., Sun, L., Ye, J., Fleisher, A., Wu, T.Reiman, E.  2010\.  Learning brain connectivity of Alzheimer’s disease by sparse inverse covariance estimation Learning brain connectivity of alzheimer’s disease by sparse inverse covariance estimation.  NeuroImage50935 - 949.
  * SG. Huang . (2021) Huang2021Huang, SG., Chung, M., Qiu, A. ADNI.  2021\.  Fast mesh data augmentation via Chebyshev polynomial of spectral filtering Fast mesh data augmentation via chebyshev polynomial of spectral filtering.  Neural Networks143198 - 208.
  * John . (2017) John2017John, AA., Javali, M., Mahale, R., Mehta, A., Acharya, PT. Srinivasa, R.  2017\.  Clinical impression and Western aphasia battery classification of aphasia in acute ischemic stroke: Is there a discrepancy? Clinical impression and western aphasia battery classification of aphasia in acute ischemic stroke: Is there a discrepancy?  Journal of Neurosciences in Rural Practice8074–078.
  * Kertesz (2007) Kertesz2007Kertesz, A.  2007\.  The Western Aphasia Battery - Revised The western aphasia battery - revised.
  * Kulkarni . (2020) Kulkarni2020Kulkarni, A., Chung, M., Bendlin, B. Prabhakaran, V.  2020\.  Investigating heritability across resting state brain networks via heat kernel smoothing on persistence diagrams Investigating heritability across resting state brain networks via heat kernel smoothing on persistence diagrams.  Proceedings of the IEEE International Symposium on Biomedical Imaging Workshops1–4.
  * Lee . (2011) Lee2011bLee, H., Chung, M., Kang, H., Kim, BN. Lee, D.  2011\.  Computing the shape of brain networks using graph filtration and Gromov-Hausdorff Metric Computing the shape of brain networks using graph filtration and gromov-hausdorff metric.  Proceedings of the International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI)14302 - 309\.
  * Lee . (2014) Lee2014Lee, H., Chung, M., Kang, H. Lee, D.  2014\.  Hole detection in metabolic connectivity of Alzheimer’s disease using k-Laplacian Hole detection in metabolic connectivity of alzheimer’s disease using k-laplacian.  International Conference on Medical Imaging Computing and Computer-Assisted Intervention (MICCAI). Lecture Notes in Computer Science (LNCS).8675297 - 304.
  * Matkowski Nikodem (1994) Matkowski1994Matkowski, J. Nikodem, K.  1994\.  An integral Jensen inequality for convex multifunctions An integral jensen inequality for convex multifunctions.  Results in Mathematics26348–353.
  * Mileyko . (2011) Mileyko2011Mileyko, Y., Mukherjee, S. Harer, J.  2011\.  Probability measures on the space of persistence diagrams Probability measures on the space of persistence diagrams.  Inverse Problems271 - 21.
  * Nichols Holmes (2002) Nichols2002Nichols, T. Holmes, A.  2002\.  Nonparametric permutation tests for functional neuroimaging: a primer with examples Nonparametric permutation tests for functional neuroimaging: a primer with examples.  Human Brain Mapping151-25.
  * Pachauri . (2011) Pachauri2011Pachauri, D., Hinrichs, C., Chung, MK., Johnson, SC. Singh, V.  2011\.  Topology-based kernels with application to inference problems in Alzheimer’s disease Topology-based kernels with application to inference problems in alzheimer’s disease.  IEEE Transactions on Medical Imaging301760 - 1770\.
  * Reininghaus . (2015) Reininghaus2015Reininghaus, J., Huber, S., Bauer, U. Kwitt, R.  2015\.  A stable multi-scale kernel for topological machine learning A stable multi-scale kernel for topological machine learning.  Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)4741–4748.
  * Rubinov Sporns (2010) Rubinov2010Rubinov, M. Sporns, O.  2010\.  Complex network measures of brain connectivity: uses and interpretations Complex network measures of brain connectivity: uses and interpretations.  NeuroImage521059 - 1069.
  * Simpson . (2013) Simpson2013bSimpson, S., Lyday, R., Hayasaka, S., Marsh, A. Laurienti, P.  2013\.  A permutation testing framework to compare groups of brain networks A permutation testing framework to compare groups of brain networks.  Frontiers in Computational Neuroscience71-13.
  * Songdechakraiwut Chung (2023) Song2023Songdechakraiwut, T. Chung, M.  2023\.  Topological learning for brain networks Topological learning for brain networks.  Annals of Applied Statistics17403 - 433.
  * Sporns (2002) Sporns2002Sporns, O.  2002\.  Network analysis, complexity, and brain function Network analysis, complexity, and brain function.  Complexity856 - 60.
  * Stolz . (2017) Stolz2017Stolz, B., Harrington, H. Porter, M.  2017\.  Persistent homology of time-dependent functional networks constructed from coupled time series Persistent homology of time-dependent functional networks constructed from coupled time series.  Chaos27.
  * Tsao . (2022) Tsao2022Tsao, CW., Aday, AW., Almarzooq, ZI., Alonso, A., Beaton, AZ., Bittencourt, MS.others  2022\.  Heart disease and stroke statistics—2022 update: a report from the American Heart Association Heart disease and stroke statistics—2022 update: a report from the american heart association.  Circulation145e153-e639.
  * Wang . (2021) Wang2021Wang, Y., Behroozmand, R., Phillip Johnson, L., Bonilha, L. Fridriksson, J.  2021\.  Topological signal processing and inference of event-related potential response Topological signal processing and inference of event-related potential response.  Journal of Neuroscience Methods363109324.
  * Wang . (2022) Wang2022Wang, Y., Chung, M. Fridriksson, J.  2022\.  Spectral permutation test on persistence diagrams Spectral permutation test on persistence diagrams.  Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP)1461-1465.
  * Wang . (2018) Wang2018Wang, Y., Ombao, H. Chung, M.  2018\.  Topological data analysis of single-trial electroencephalographic signals Topological data analysis of single-trial electroencephalographic signals.  Annals of Applied Statistics121506-1534.
  * Wang . (2019) Wang2019Wang, Y., Ombao, H. Chung, M.  2019\.  Statistical persistent homology of brain signals Statistical persistent homology of brain signals.  Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP)1125 - 1129.
  * Winkler . (2016) Winkler2016Winkler, A., Ridgway, G., Douaud, G., Nichols, T. Smith, S.  2016\.  Faster permutation inference in brain imaging Faster permutation inference in brain imaging.  NeuroImage141502 - 516.
  * Yourganov . (2018) Yourganov2018Yourganov, G., Fridriksson, J., Stark, B. Rorden, C.  2018\.  Removal of artifacts from resting-state fMRI data in stroke Removal of artifacts from resting-state fmri data in stroke.  NeuroImage: Clinical17297-305.