# A selective review of recent developments in spatially variable gene
detection for spatial transcriptomics

Sikta Das Adhikari1,2,  Jiaxin Yang2 ,  Jianrong Wang2∗ ,  Yuehua Cui1
To whom correspondence should be addressed. Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
1Department of Statistics and Probability, Michigan State University, East
Lansing, MI 48824, USA. and 2Department of Computational Mathematics, Science
and Engineering, Michigan State University, East Lansing, MI 48824, USA.

###### Abstract

_With the emergence of advanced spatial transcriptomic technologies, there has
been a surge in research papers dedicated to analyzing spatial transcriptomics
data, resulting in significant contributions to our understanding of biology.
The initial stage of downstream analysis of spatial transcriptomic data has
centered on identifying spatially variable genes (SVGs) or genes expressed
with specific spatial patterns across the tissue. SVG detection is an
important task since many downstream analyses depend on these selected SVGs.
Over the past few years, a plethora of new methods have been proposed for the
detection of SVGs, accompanied by numerous innovative concepts and
discussions. This article provides a selective review of methods and their
practical implementations, offering valuable insights into the current
literature in this field._

_Keywords_  _Spatial transcriptomics; Spatially variable genes; Spatially
resolved transcriptomics; Single cell RNA sequencing_

##  1 Introduction

Recent advancements in Spatially-resolved transcriptomics (SRT) technology
have provided comprehensive gene expression data for thousands of genes across
multiple samples or spatial spots, accompanied by their respective spatial
coordinates. Depending on the specific technology utilized, a sample could
represent a single cell (as in the case of STARmap technology), a cell-sized
local region (as with HDST technology[1]), or a localized region comprising
dozens of cells (as seen in Slide-seq[2, 3] and Visium technologies). The
latest SRT platforms, such as 10x Genomics Visium and Slide-seqV2, encompass
thousands of spatial locations within each tissue sample, with future
developments poised to achieve even higher resolutions. As technology
progresses, the demand for more robust statistical frameworks to effectively
analyze spatial data intensifies.

Although spatial transcriptomic (ST) data permit addressing a range of
distinct questions, a fundamental initial step in the downstream analysis of
spatial data is the identification of spatially variable genes (SVGs). These
are genes that exhibit variations in expression levels either across the
entire tissue or within predefined spatial domains. These genes can
potentially unveil tissue heterogeneity and the underlying structural factors
that drive distinct expression patterns across spatial locations, thus
offering valuable insights into biology.

Numerous methods have been developed for the identification of SVGs. These
methods encompass a spectrum of approaches, including the utilization of
standard spatial statistics measures like Moran’s I statistic[4] and Geary’s C
statistic[5] to rank genes based on their spatial autocorrelation. More
advanced methods employ model-based approaches such as SpatialDE[6],
SpatialDE2[7], SPARK and its extensions[8], nnSVG[9], BOOST-GP[10], marked
point process frameworks like Trendsceek[11] and scGCO[12], or model-free
frameworks like sepal[13] and GLISS[14]. Additionally, there are toolboxes,
such as MERINGUE[15], Giotto[16], Seurat[17] that integrate some of these
methods into comprehensive end-to-end analysis frameworks.

Downstream analysis involving SVGs encompasses various tasks, such as spatial
clustering, deciphering spatial domains, and identifying spatial domain-
specific SVGs. Additionally, there are numerous other downstream analyses that
leverage additional information like scRNASeq data, histological images, and
more, for tasks such as spatial decomposition of spots, gene imputation, the
inference of cell-cell and gene-gene interactions and spatial location
reconstruction for scRNA-seq data. However, this review primarily concentrates
on SVG detection frameworks and does not delve into these other downstream
analyses.

Thus, it is the primary focus of this paper to discuss selected frameworks for
SVG identification, serving as a valuable resource for researchers new to this
field, enabling them to become acquainted with existing SVG identification
frameworks, their unique characteristics, novelty, as well as their pros and
cons.

##  2 Overview of the frameworks for detecting SVGs

Generally, in a spatial transcriptomics setup, the available spatial dataset
contains gene expression measures/counts for $m$ genes distributed across $N$
known spatial coordinates or spots. This section establishes the key symbols
that will be frequently utilized. Specifically, $y=(y_{1},y_{2},...,y_{N})$ is
defined as the gene expression profiles/counts for a given gene across spatial
coordinates (referred to as samples or spots), denoted by
$s=(s_{1},...,s_{N})$. The coordinates of the spatial locations are typically
two-dimensional, i.e., $s_{i}=(s_{i1},s_{i2})$, but any dimensional
coordinates can be employed. The primary objective of these SVG detection
models is to ascertain which genes, out of the $m$ genes, are spatially
variable across the tissue. In other words, the main goal is to determine
whether the gene expression measure $y$ depends on or relates to the spatial
locations where the gene expression measures are collected.

###  2.1 Gene expression data and pre-processing step

The gene expression measure $y$ are generally of count data type from
sequencing reads. Various SVG detection models have been developed to
specifically model count data following some mandatory filtering and quality
control steps. Some examples of these models include SPARK-X[18], BOOST-
GP[10], SINFONIA[19], and GPcounts[20]. The gene expression count data often
exhibit over-dispersion and contain numerous zero values, mainly due to the
technology employed for data generation or simply because many genes are
poorly expressed for biological reasons. These particular issues in count data
are generally taken care of by using negative binomial models which handle
over-dispersion well. For the issue of zero-inflation, Zhao et al, 2022 [21]
showed that modeling zero inflation is not necessary in spatial
transcriptomics, thus is not a concern in many method development. On the
other hand, some methods, for example SpatialDE[6], nnSVG[9], and BOOST-
MI[22], use normalized gene expression data in the model for easy
implementation, where in most of cases, the data is modeled using multivariate
normal distribution after transformation. Authors in SPARK[8] proposed two
different data models, SPARK and SPARK-G which uses count data and normalized
data, respectively. The data normalization method is not unique for these
methods. The normalization step generally removes the bias due to differences
in sequencing depth using size factors and normalizes the data using log
transformation(log10 or log2 transformations after adding a pseudo-count value
$c$, preferably 1). The method sepal[13] uses a slightly different
normalization procedure which involves mapping the log-transformed values to
the interval [0,1] and using a pseudocount 2. Other normalization methods,
such as scran, scuttle, and scater R/Bioconductor packages[23, 24], can also
be applied. Table 1 provides information on some selective methods together
with their required input data type and the implemented model:

Table 1: A selective list of methods for SVG detection in ST data analysis along with the required input data type and the implemented model. Method | Input data type | Data model
---|---|---
Trendsceek[11] | Normalized | Marked point processes
SpatialDE[6] | Normalized | Multivariate Normal
SpatialDE2[7] | Count | Poisson
SPARK[8] | Count | Overdispersed poisson
SPARK-G[8] | Normalized | Multivariate Normal
SPARK-X[18] | Count | Non-parametric model
nnSVG[9] | Normalized | Multivariate Normal
BOOST-GP[10] | Count | Zero-inflated negative binomial
SOMDE[25] | Normalized | Multivariate Normal
sepal[13] | Normalized | Model-free
GLISS[14] | Normalized | Model-free
SINFONIA[19] | Count | Model-free
BOOST-MI[22] | Normalized | Modified Ising model
scGCO[12] | Normalized | Marked point process
CTSV[26] | Count | Zero-inflated negative binomial
HEARTSVG[27] | Count | Model-free
MULTILAYER[28] | Normalized | Model-free
GPcounts[20] | Count | Negative binomial

###  2.2 Overview of model-based frameworks

####  2.2.1 Gaussian process(GP) regression based and similar models

The majority of the methods, including some of the state-of-the-art algorithms
to detect SVG, are based on Gaussian process (GP) regression models. For
example, one of the first published SVG detection methods, SpatialDE[6],
models the normalized gene expression $y$ for a given gene using the following
multivariate normal model:

$p(y|\mu,\sigma_{s}^{2},\delta,R)\sim N(y|\mu 1,\sigma_{s}^{2}K+\delta I),$
(1)

where the covariance term is decomposed into a spatial and a non-spatial part,
where $\delta I$ represents the non-spatial part and $\sigma_{s}^{2}K$ is the
spatial covariance matrix, whose $(i,j)^{th}$ element in the kernel matrix $K$
denotes the spatial similarity between the $i^{th}$ and $j^{th}$ spot
calculated based on the corresponding coordinates $s_{i}$ and $s_{j}$. The
choice of the kernel function plays a very important role in detecting the
spatial correlation present in the gene expressions. More discussion about
kernel function is provided in the next subsection.

Other methods like SPARK-G[8] (the Gaussian version of SPARK), nnSVG[9], and
SOMDE[25] implement similar GP models for modeling normalized gene expression
data with some extra features or added level of complexity. SPARK-G and nnSVG
provide the option to include extra covariate terms in the model. The
covariates or the explanatory variables could contain batch information, cell-
cycle information, or other information that is important to adjust for during
the analysis. SOMDE is a two-step procedure. This approach involves first
utilizing a self-organizing map to cluster neighboring cells into nodes.
Subsequently, it employs a Gaussian process to model and analyze the spatial
gene expression patterns at the node level.

Table 1 shows that methods like SPARK[8], SpatialDE2[7], BOOST-GP[10],
CTSV[26], and GPcounts[20] model count data directly. SPARK models the count
data using an overdispersed poisson model where the logarithm of the unknown
Poisson rate parameter is assumed to follow a stationary Gaussian process with
similar spatial and non-spatial covariance components. BOOST-GP presents a
novel Bayesian hierarchical model to analyze spatial transcriptomic data,
which models the count data using a zero-inflated negative binomial(ZINB)
model. The logarithm of the normalized expression level, which is included in
the expectation term in NB, can be seen as a GP with a spatial covariance term
representing the spatial variability in case there is a spatial pattern.
GPcounts also uses negative binomial distribution to model the UMI(Unique
Molecular Identifier) data. SpatialDE2 employs a Generalized Linear Mixed
Model (GLMM) for count data modeling. In contrast to GP-based techniques that
typically separate covariance into a spatial and a non-spatial component,
SpatialDE2 dissects the covariance into several spatial components along with
a non-spatial random component. CTSV implements a slightly different technique
and does not use the GP model. In CTSV, the gene specific, spot specific and
cell-type specific relative mean expression level in the ZINB model is a
linear combination of $h_{1}(s_{.1})$ and $h_{2}(s_{.2})$ where the functions
$h_{1}(\cdot)$ and $h_{2}(\cdot)$ represents the underlying true spatial
effect modeled with the kernel function in GP model.

####  2.2.2 Statistical inference and selecting kernel function in GP-based
frameworks

Typically, when evaluating the existence of spatial patterns within the data,
an assessment is made by testing the alternative hypothesis, which suggests
the presence of a spatial covariance term in the model, against the null
hypothesis, where the spatial covariance term is set to zero, indicating the
absence of spatial variability. This comparison between the model fitted under
the alternative hypothesis and the null model forms the basis of a
significance testing procedure. This often involves conducting significance
tests and drawing conclusions based on p-values in frequentist approaches. For
example, in model (1), testing SVG is equivalent to testing
$H_{0}:\sigma_{s}^{2}=0$.

As previously mentioned, selecting the appropriate kernel function for
computing the spatial covariance matrix is a critically important step in
identifying spatial patterns within the data. Ideally, the kernel function
should accurately capture the true underlying relationship between the $y$
values and the spatial coordinates $s$. In practice, the actual underlying
function remains unknown, and the closer the chosen kernel function
approximates the true functional form, the more precise the model
specification becomes, rendering the test more robust and powerful.

SpatialDE employs a squared exponential covariance function(a.k.a Gaussian
kernel function or radial basis kernel function) to compute the spatial
covariance matrix:

$K_{i,j}=k(s_{i},s_{j})=exp(-\frac{|s_{i}-s_{j}|^{2}}{2l^{2}})$

The hyperparameter $l$, recognized as the characteristic length scale or
bandwidth, determines how rapidly the covariance decays as a function of
distance and is typically chosen by grid search. SOMDE also uses the squared
exponential (Gaussian) kernel in their model with ten different length scales
and chooses the one that achieves the highest log-likelihood ratio value.
GPcounts uses linear or periodic kernel based on BIC values. SPARK asserts
that relying on a single kernel restricts the ability to robustly identify
spatially variable genes across diverse spatial patterns. Therefore, SPARK
(and SPARK-G) adopts an approach involving a total of ten distinct spatial
kernels. These comprise five periodic kernels (e.g., Cosine kernels) with
varying periodicity parameters and five Gaussian kernels with different
smoothness parameters. SPARK proceeds to compute ten p-values, each derived
from a different test employing these various kernel functions. These p-values
are subsequently combined using the Cauchy combination rule [29, 30]. Similar
to SPARK, SpatialDE2 incorporates a variety of pre-defined kernels with
varying structures and length scales. It also offers the flexibility to
conduct an omnibus test as an alternative to independently testing each kernel
and subsequently merging the p-values. nnSVG posits that genes can potentially
display a vast spectrum of spatial patterns, and using the same set of kernel
functions for all of the genes might lead to less powerful tests. They
consider the use of an exponential covariance function as a kernel function
where the length scale parameter of the kernel function is fitted for each
gene, which allows capturing the unique spatial variability pattern of the
gene. CTSV uses five different sets of functional forms for $h_{1}(s_{.1})$
and $h_{2}(s_{.2})$, which includes linear functions, squared exponential
functions, and periodic functions with different sets of scaling parameters
and the five p-values calculated from five different forms are combined using
the Cauchy combination method.

Although different models discussed here have some similarities in testing
procedures, the model fitting techniques implemented and the testing
procedures utilized are different and are summarized in Table 2.

Table 2: List of some popular SVG detection methods with model-fitting and testing information. Method | Bayesian/ | Model fitting and | Hypothesis testing
---|---|---|---
Frequentist | parameter estimation | method
SpatialDE | Frequentist | Maximizing marginal | LogLikelihood ratio test
|  | log likelihood |
SpatialDE2 | Frequentist | Only null model parameters needs | Score test based
|  | to be estimated by BLUP | on Zhang and Lin[31]
SPARK | Frequentist | Approximate-inference algorithm | Satterthwaite method
|  | based on the PQL approach | on the basis of score statistics
SPARK-G | Frequentist | Maximum likelihood | Score test
nnSVG | Frequentist | Fast optimization algorithms for | Likelihood ratio test
|  | NNGP models (BRISC R package) |
BOOST-GP | Bayesian | Sampling from posterior | Bayes Factor or posterior
|  | using MCMC | probabilities of inclusion (PPI)
SOMDE | Frequentist | Gradient optimization | Likelihood ratio test
CTSV | Frequentist | Approx. maximum likelihood | Wald tests (R
|  | using conjugate gradient(CG) | package pscl)
|  | algorithm |
GPcounts | Frequentist | Optimization of log marginal | Likelihood ratio test
|  | likelihood by variational |
|  | approximation |

The statistical power of GP-based methods hinges on the selection of kernel
functions[8], which can complicate the model selection and limit SVG detection
power. To address this challenge, the authors in [22] introduced BOOST-MI.
This novel approach employs Bayesian modeling of spatial transcriptomics data
via a modified Ising model to identify SV genes. As an initial step, BOOST-MI
takes normalized gene expression data as input and dichotomizes the normalized
expression levels into a binary spatial pattern. Subsequently, BOOST-MI
proceeds to identify a wide spectrum of spatial patterns displayed by the
genes by inferring the Ising model interaction parameter within a Bayesian
framework. It achieves this by generating samples from the posterior
distribution of the parameters through a double Metropolis-Hastings (DMH)
algorithm[32]. Subsequently, it computes the Bayes factor based on these
posterior samples, which are then used for selecting SV genes.

Trendsceek[11], one of the earliest published SVG detection methods, models
data as marked point processes, where they assign points to represent the
spatial locations of spots and marks on each point to represent expression
levels. The pivotal objective of Trendsceek revolves around evaluating the
dependency between the spatial distribution of points and their respective
marks through pairwise analyses as a function of the inter-point distances.
The underlying premise is that if there exists no dependency between marks and
point locations, the resulting scores should remain constant across various
distances. A resampling procedure is executed to gauge the significance of a
gene’s spatial expression pattern, involving permutations of expression values
that create a null model with no spatial expression dependency.

Similar to Trendsceek, ScGCO[12](single-cell graph cuts optimization) method
also models gene expression data as a marked point process where points
represent the spatial locations of measured cells or spots, and marks are
discrete gene expression states (such as, down-regulated or up-regulated)
associated with points. It analyzes the dependency of points with a specific
mark on spatial locations using a hypothesis test. Under the null hypothesis
(i.e., no spatial dependency), it assumes that points with a specific mark in
a 2D space are distributed in a completely random fashion and can be described
by a homogeneous spatial Poisson process. Genes with spatial regions whose
number of cells/spots of specific marks are associated with statistically
significant low probabilities under the null model are selected as SVG.

###  2.3 Overview of model-free frameworks

There are other SVG detection methods such as SPARK-X[18], sepal[13],
GLISS[14], and SINFONIA[19] which do not attempt to model the data generation
process or rely on distributional assumptions. Instead, they use model-free
techniques to detect SVGs. The authors introduced sepal[13] (Spatial
Expression Pattern Locator), an innovative method that leverages transcript
diffusion simulations to identify genes exhibiting spatial patterns. It
simulates transcript diffusion within the spatial domain and measures the time
required for convergence. The core idea is that transcripts with random
spatial distributions will converge more quickly or reach a homogeneous state
faster compared to those with distinct spatial patterns. Consequently, the
diffusion time serves as an indicator of a gene’s degree of spatial
variability. Genes with longer diffusion times exhibit less spatial
randomness. Therefore, ranking genes based on this indicator and selecting the
top-ranked genes as SVGs is a logical approach.

SINFONIA[19] offers a scalable approach to initially identify spatially
variable genes through ensemble strategies as part of its spatial
transcriptomic data analysis, with the ultimate goal of deciphering spatial
domains. SINFONIA initially identifies the $k$ nearest neighbors in Euclidean
space for each spot and builds a Spatial Neighbor Graph (SNG) using the weight
matrix where the ($i,j$)th element is determined by a function of the distance
between the $i$th and $j$th spot. Next, SINFONIA calculates Moran’s I and
Geary’s C statistics for each gene based on the weight matrix $W$ to assess
spatial autocorrelation. The underlying concept is that genes with more
pronounced spatial autocorrelation exhibit more organized spatial expression
patterns.

HEARTSVG[27] utilizes a unique, distribution-free, test-based approach that
focuses on identifying non-SVGs first and then infers the presence of SVGs
using this information. The process involves assessing serial autocorrelations
within the marginal expressions across the global spatial context to pinpoint
non-SVGs. This, in turn, enables the automatic recognition of all other genes
as SVGs, regardless of their spatial patterns. HEARTSVG asserts its
superiority in terms of robustness and computational efficiency by abstaining
from assumptions about specific underlying spatial patterns for these variable
genes.

SPARK-X[18] is a nonparametric method grounded in the following insight: if
$y$ is independent of $s$, then the spatial distance between two locations $i$
and $j$ would also be unrelated to the gene-expression difference between
those two locations. SPARK-X constructs two $N\times N$ projection covariance
matrices: (1) The expression covariance matrix based on gene expression
levels; and (2) the distance covariance matrix based on all spatial locations.
It employs a test statistic derived from the product of these two covariance
matrices to evaluate the independence between the gene expression ($y$) and
the spatial coordinates ($s$). In simpler terms, if gene expressions are
indeed independent of spatial coordinates, the product of these covariance
matrices will yield a small value. Conversely, if gene expressions are
dependent on the spatial coordinates, the product of the matrices will yield a
large value.

Similar to the kernel matrix used in methods like SpatialDE or SPARK, the
statistical power of the SPARK-X test inevitably hinges on how the distance
covariance matrix is constructed and how well it aligns with the true
underlying spatial patterns exhibited by the gene of interest. To ensure
robust identification of spatially varying genes across diverse spatial
expression patterns, SPARK-X explores various transformations of the spatial
coordinates ($s$) and subsequently generates distinct distance covariance
matrices. Specifically, the algorithm applies five Gaussian transformations
with varying smoothness parameters and five cosine transformations to the
spatial coordinates ($s$). This process results in the creation of eleven
distinct p-values, corresponding to the ten transformed distance covariance
matrices and the original one constructed using the original coordinates.
These individual p-values are then combined using the Cauchy combination
method.

MULTILAYER[28] treats spatially transcriptomics data as a raster image and
uses digital image strategies to resolve tissue substructures. The basic unit
in MULTILAYER is the "gexel", gene expression element analogous to a pixel in
a digital image. The gene expression levels per gexel relative to the average
gene expression are computed within the tissue. Genes are considered
upregulated or downregulated when their normalized read counts per gene are
above or below the average behavior, respectively. Differentially expressed
genes are ranked based on the number of related gexels, providing a rapid view
of genes that are overrepresented on the digital map based on their relative
expression.

GLISS[14] (Graph Laplacian-based Integrative Single-cell Spatial Analysis)
utilizes a graph-based feature learning framework to detect and discover SVGs
and recover cell locations in scRNA-seq data by leveraging spatial
transcriptomic and scRNA-seq data. The workflow involves multiple steps.
First, SV genes are identified from ST data using graph-based feature
selection. Next, it determines the cells of interest in the scRNA-seq data
based on unsupervised learning methods and leverage these selected SVGs to
discover new SVGs in scRNA-seq data. The final goal of this workflow is to
cluster genes based on their spatial patterns.

The BPS (Big-Small Patch)[33] method, introduced in a recent publication,
utilizes a non-parametric model for the identification of spatially variable
genes in 2D or 3D spatial transcriptomics data. The approach involves taking
normalized spatial transcriptomics data as input. It defines big and small
patches for each spatial spot based on neighboring spots with larger or
smaller radii, respectively. The method then calculates local means of gene
expression for both big and small patches. Following this, it calculates the
ratio between the variances of local means for each gene, approximating a log-
normal distribution for the distribution of these ratios. Subsequently, a
p-value is determined for each gene based on this approximated distribution.

##  3 Statistical Inference with Multiple Testing Control

We have previously discussed both model-based and model-free methods for
detecting SVGs. The mathematical models employed for capturing the data
generation process and the innovative model-free SVG detection technique have
proven valuable for uncovering significant SVGs that offer critical biological
insights. However, from a statistical perspective, concerns arise regarding
the potential for false discoveries of genes that lack genuine spatial
variability. This concern becomes more pronounced when a large number of genes
are simultaneously tested across most frameworks. If the false discovery rate
or type 1 error is not adequately controlled, it may lead to incorrect
conclusions and the selection of numerous genes that exhibit false spatial
variability.

Various methods have been developed for multiplicity correction (MC) to
address this concern. Some methods analytically constrain the false discovery
rate (FDR) to remain below a predetermined threshold, while others do not
analytically control the FDR and simply select a user-specified number of top
genes as SVGs. Researchers may choose a method that aligns better with their
research goals and the type of downstream analysis they intend to perform. In
Table 3, we present an overview of these methods, organized around these
critical questions. The permutation-based method is usually considered as the
golden standard method as it is purely data-driven and distribution free.
However, it is the least scalable one since it is computationally more
demanding. The FDR-based methods have been the commonly applied ones since
they offer type I error control while maintaining high power compared to the
Bonferroni method. Nevertheless, depending on the downstream analysis goal, it
is not necessary to strictly enforce the MC rule. For example, when the goal
is to find the low dimensional embedding of genes, such as in spatial PCA
analysis [34] people usually choose top ranked genes for further analysis. In
such cases, strictly enforcing MC is not needed.

Table 3: List of methods and the procedures used for multiple testing control. Method | If framework analytically | How SVGs are selected
---|---|---
controls FDR
Trendsceek | Yes | Permutation based p-values,
|  | Benjamini–Hochberg procedure[35] for MC
SpatialDE | Yes | Analytically estimated p-values,
|  | q-value method[36] for MC
SpatialDE2 | Yes | Analytically estimated p-values,
|  | Benjamini–Yekutieli procedure[37] for MC
SPARK | Yes | Analytically estimated p-values,
|  | Benjamini–Yekutieli procedure for MC
SPARK-G | Yes | same as SPARK
SPARK-X | Yes | same as SPARK
nnSVG | Yes | Analytical approximate p-values,
|  | Benjamini–Hochberg method for MC
BOOST-GP | Yes | Based on Bayesian FDR controlled PPI threshold
SOMDE | No | Top ranked genes based on
|  | spatial variability score
sepal | No | Top $k$ genes with highest ranks
GLISS | Yes | Permutation based p-values,
|  | Benjamini–Hochberg procedure for MC
SINFONIA | No | Top $k$ genes with highest score
|  | and an ensemble technique
BOOST-MI | No | Based on specific Bayes Factor threshold
scGCO | Yes | Analytically estimated p-values,
|  | Benjamini–Hochberg procedure for MC
CTSV | Yes | Analytically estimated p-values,
|  | q-value method for MC
HEARTSVG | Yes | Analytically estimated p-values,
|  | MC by Bonferroni/Holm/Hochberg
MULTILAYER | No | Based on the two-fold threshold
|  | of a test statistic
GPcounts | Yes | Analytical or permuted p-values,
|  | q-value method for MC

##  4 Exploring Performance, Advantages, and Limitations

In the preceding sections, we have explored the complexities associated with
spatial count data. In many instances, these count data are characterized by
sparsity and overdispersion. Section 2 of this review classifies modeling
frameworks based on whether they directly model the count data or opt for
modeling the normalized data. Some literatures [18, 10] argue against modeling
normalized data with a Gaussian distribution due to concerns that such a
parametric approximation may result in reduced statistical power and
difficulties in controlling type 1 errors, especially when dealing with small
p-values.

On the other hand, methods that employ normalized count data, such as
SpatialDE, SPARK-G, and nnSVG, offer advantages, including simpler model
structures and reduced computational challenges. Notably, SPARK employs a dual
modeling approach, encompassing both an overdispersed Poisson model (SPARK)
and a Gaussian model (SPARK-G) for count data analysis. They declare that
SPARK-G exhibits significantly improved computational efficiency compared to
the Poisson version of SPARK. Moreover, SPARK-G may demonstrate greater
resilience to model misspecification, potentially enhancing its effectiveness
in specific data applications.

Although many researchers prefer to model count data directly, there is no
consensus on the preferred approach for directly modeling count data either.
While some opt for Poisson distribution models, others argue that it may be
insufficient to address issues of overdispersion, suggesting that a negative
binomial distribution is more suitable in such cases. Furthermore, when data
exhibit extreme sparsity, the utilization of a zero-inflated Poisson or
negative binomial model may be more logical, although it tends to introduce
greater complexity into the model. But we need to note that direct modeling of
sparse count data with a negative binomial distribution or other over-
dispersed Poisson distributions incurs algorithm stability issues [18, 38,
20].

With the continuous evolution of spatial transcriptomic technologies,
researchers now have access to increasingly vast and high-resolution spatial
datasets. Analyzing these extensive datasets demands the use of efficient and
scalable methods for downstream analysis. Notably, approaches like Trendsceek
and BOOST-GP impose substantial computational demands. In a study referenced
from SRTsim[39], it was observed that when applying these methods to synthetic
data, Trendsceek (v.1.0.0) required approximately 10 hours, while BOOST-GP
needed about 8 hours to analyze a single synthetic dataset containing 1000
genes and 673 locations. In the same research context, SOMDE (v.0.1.8)
struggled, failing to process nearly 90 percent of the genes and yielding NA
values.

Another comprehensive comparison, outlined in a review paper[40], assessed the
performance of various SVG detection methods. The evaluation considered
computational time and memory usage across 20 diverse spatial datasets, each
varying in the number of spots or samples. Among the methods examined,
including SpatialDE, SPARK-X, nnSVG, SOMDE, Giotto KM, and Giotto Rank (both
are implemented in the Giotto package), SPARK-X emerged as the swiftest, with
SOMDE following as the second-best option, albeit notably slower than SPARK-X.
SpatialDE exhibited poorer performance in larger datasets, while nnSVG proved
faster than SpatialDE for larger datasets but relatively slower for datasets
with fewer spatial locations. In particular, SPARK-X [18] scales linearly with
the number of spatial locations, while other methods scale cubically (e.g.,
SpatialDE and SPARK) or quadratically (SpatialDE2).

In terms of peak memory usage, study [40] revealed that SOMDE consumed the
least memory, with SPARK-X ranking second. Conversely, SpatialDE demonstrated
high peak memory consumption. Considering the trade-off between speed and
memory usage, SPARK-X and SOMDE emerged as the two most efficient methods, as
determined by the experiment. Furthermore, the evaluation included other
methods such as Giotto KM, Giotto Rank, and Moran’s I , but none of these
alternatives matched the efficiency of SPARK-X or SOMDE based on the
experimental findings.

In summary, each modeling framework comes with its own set of pros and cons,
necessitating careful consideration of the trade-off between computational
efficiency/cost and performance when selecting the most suitable approach for
analyzing spatial count data. The model-free or nonparametric approaches do
not try to capture the data generation process and offer alternative
frameworks to detect SVG. Most of the method frameworks are very intuitive but
each comes with its own sets of restrictions or assumptions. For example,
Trendsceek is a resampling-based method, which incurs a substantial
computational load, rendering its application impractical for extensive ST
datasets. SPARK-X exhibits impressive performance for high dimensional data,
but the authors recommend using it with large sample (e.g., spot) size, say
3,000 or more.

##  5 Assessing Input Data and Model Outputs

For the various methodologies we reivewed so far, some of these approaches
primarily focus on identifying genes that exhibit spatial variability across
the entire tissue, exemplified by methods like SpatialDE and SPARK. In
contrast, others are additionally equipped to detect genes with spatial
variability within predefined spatial domains, as seen in nnSVG. Other methods
like SpaGCN[41] and STAMarker[42] are designed to identify Spatial domains and
detect SVGs within spatial domains.

Additionally, certain methods aim to identify SVGs to facilitate downstream
analysis. For instance, SINFONIA, as cited in this work, provides a scalable
approach for the initial identification of spatially variable genes using
ensemble strategies within the context of spatial transcriptomic data
analysis. The ultimate objective of this method is to decipher distinct
spatial domains within the tissue.

Furthermore, some of these methods leverage additional information as input,
such as single-cell RNA sequencing (scRNA-seq) data, spatial domain
information, or tissue-specific markers, in conjunction with spatial
transcriptomic data. For instance, CTSV requires scRNA-seq data and a set of
marker genes as input alongside the spatial transcriptomic data. It employs
deconvolution techniques like SPOTlight[43], RCTD[44], or SpatialDWLS[45] to
estimate cell-type proportions for each spatial spot. Ultimately, this
approach identifies spatially variable genes specific to different cell types.
Similarly, Trendsceek identifies genes with significant spatial trends and
subsequently determines the subset of cells occupying spatial regions of
interest.

Given the distinct ultimate objectives and input criteria for each method, it
would be unfair to evaluate their performance solely based on a single
parameter. Rather, the utility or superiority of these frameworks depends on
the researcher’s specific goals and the nature of their research inquiries. In
this context, we present a table that combines these selective frameworks,
including details about their typical inputs and primary research objectives
(see Table 4).

Table 4: List of selective methods with input data type and main goal. Method(publication) | Input data | Main goal
---|---|---
SpatialDE(2018) | ST data | Finding SVG
|  | Spatial gene-clustering
SpatialDE2(Archived,2021) | ST data | Tissue region segmentation
|  | Finding SVG
|  | Spatial gene-clustering
Trendsceek(2018) | ST data | Finding SVG
| scRNA-Seq data | Identifying cells in spatially
|  | significant gene expression regions
SPARK(2020) | ST data | Finding SVG
SPARK-X(2021) | ST data | Finding SVG
nnSVG(2023) | ST data | Finding SVGs across tissue
|  | or within spatial domains
BOOST-GP(2021) | ST data | Finding SVG
SOMDE(2021) | ST data | Finding SVG
sepal(2021) | ST data | Finding SVG
|  | Spatial gene-clustering
GLISS(Archived,2020) | ST data | Finding SVG, recovering cell
| scRNA-seq data | locations in scRNA-seq data
|  | and gene-clustering
SINFONIA(2023) | ST data | Finding SVG for
|  | deciphering spatial domains
BOOST-MI(2022) | ST data | Finding SVG
scGCO(2022) | ST data | Finding SVG
CTSV(2022) | ST data | Detecting cell-type-specific
| scRNA-seq | SVG
| set of marker genes |
HEARTSVG (Archived,2023) | ST data | Detecting SVG and spatial domain
MULTILAYER(2021) | ST data | Detecting SVG, dimensionality
|  | reduction, spatial clustering and more
GPcounts(2021) | ST data | Finding SVG, identifying gene-specific
| scRNA-seq data | branching locations and more
SpaGCN(2021) | ST data | Identifying spatial domains
| histology image data | and SVG in domain
STAMarker(Archived,2022) | ST data | Spatial domain-specific variable genes

##  6 Publicly Accessible Code for Major SVG Detection Methods

Every prominent SVG detection method featured in this paper has made its code
publicly available. Certain frameworks have packages published in CRAN or
available as Python modules, while others have shared their code on Github,
and the package can be installed directly from Github. Here, we have compiled
a list of the packages and repositories associated with these techniques,
along with the coding language they have used (see Table 5). This compilation
aims to facilitate convenient access to their respective code bases, making it
easier for researchers to choose a method based on their preferred programming
language.

Table 5: List of methods with implementing code language and package site. Method | Code language | Package or GitHub or vignette
---|---|---
Trendsceek | R | https://github.com/edsgard/trendsceek
SpatialDE | Python | https://github.com/Teichlab/SpatialDE
SpatialDE2 | Python | https://github.com/PMBio/SpatialDE
SPARK | R | https://github.com/xzhoulab/SPARK
SPARK-G |  | https://xzhoulab.github.io/SPARK/01_about/
SPARK-X |  |
nnSVG | R | https://bioconductor.org/packages/release/bioc/html/nnSVG.html
|  | https://github.com/lmweber/nnSVG
BOOST-GP | R/C++ | https://github.com/Minzhe/BOOST-GP
SOMDE | Python | https://pypi.org/project/somde
|  | https://github.com/XuegongLab/somde
sepal | Python | https://github.com/almaan/sepal
|  | 10.5281/zenodo.4573237
GLISS | Python | https://github.com/junjiezhujason/gliss
SINFONIA | Python | https://github.com/BioX-NKU/SINFONIA
BOOST-MI | R | https://github.com/Xijiang1997/BOOST-MI
ScGCO | Python | https://github.com/WangPeng-Lab/scGCO
CTSV | R | https://bioconductor.org/packages/devel/bioc/html/CTSV.html
HEARTSVG | R | https://github.com/cz0316/HEARTSVG.git
MULTILAYER | Python | https://github.com/SysFate/MULTILAYER
GPcounts | Python | https://github.com/ManchesterBioinference/GPcounts
BPS | Python | https://github.com/juexinwang/BSP/

##  7 Concluding Remarks

We systematically reviewed recently developed frameworks for identifying
spatially variable genes and grouped them into different categories and delved
into the unique aspects of their models and underlying principles. Here, we
provide a brief discussion encompassing various facets, including
preprocessing steps, modeling frameworks, inference techniques, scalability,
and practical applicability of these frameworks. We explored the performance
of select methods as reported in previously published papers. Nevertheless, it
is essential to note that we refrained from conducting evaluations based
solely on the number of SVGs detected or the trade-off between statistical
power and FDR. This decision arises from the fact that the methods discussed
in this paper often serve different research objectives, each tailored to
specific research questions. For example, a method primarily focused on
spatial clustering may yield similar outcomes when considering the top 100
genes versus the top 110 genes. In contrast, a method geared toward accurately
identifying genuine SVGs and scrutinizing individual SVGs to glean deeper
insights into biological mechanisms may prioritize stringent control of false
discovery rates, making it a pivotal concern in their evaluation. The
evaluation criteria must align with the unique goals and nuances of each
method, akin to comparing apples to oranges when attempting to gauge their
performance solely based on the number of SVGs selected.

This paper[40] has previously investigated several methods for detecting
spatial gene expression variations and benchmarked their performance based on
different measures. It reported that, although each SVG detection method
successfully identifies a significant number of SVGs, there is limited overlap
in the SVGs detected when a significance cutoff is applied to filter the SVGs.
The study’s simulation analysis revealed that, in most cases, the estimated
FDRs do not accurately reflect the true FDRs. These findings indicate that
there is room for improvement in the commonly used methods for SVG detection
and their associated FDR control approaches.

In the context of Gaussian process based methods, one potential issue could be
related to the selection of kernel function. For instance, as an improvement
to spatialDE, approaches like SPARK and SPARK-X employ a variety of different
kernels to robustly identify various traits. However, they apply the same set
of parameter values to all genes, even when these genes may exhibit vastly
different spatial patterns. While nnSVG offers improvements by allowing gene-
specific kernel function parameter selection, it relies on a single type of
kernel function. This opens room for further methodological development for
optimal kernel selection when kernel-based methods are applied for SVG
detection.

Furthermore, model-free techniques, in many cases, do not analytically control
FDR, making it challenging to establish a specific cutoff for selecting SVGs.
Many methods claim to detect more SVGs than others, often undetected by
alternative methods. However, the mere detection of more SVGs does not
necessarily indicate the superiority of a framework if it does not effectively
control the FDR. If the goal is to pinpoint the top $k$ (say 1000) SVGs for
subsequent analysis without the necessity of precisely quantifying detection
uncertainty, these methods can be employed. However, for a more rigorous
approach, it is crucial to implement stringent FDR control measures to prevent
false discoveries. In our empirical analysis, we observed that numerous
methods exhibit elevated false positive rates with inflated p-values (data not
shown). There is an urgent demand for the development of more rigorous
statistical approaches to enhance false positive control.

In summary, we have provided a selective survey of recently published and
archived literature on SVG detection, offering an analysis of their practical
utility, adaptability, innovation, and constraints from various practical
perspectives. This effort aims to facilitate new researchers in gaining a
holistic understanding of the available methods and assist them in selecting a
framework aligned with their specific research needs and questions.

## Acknowledgement

This work was supported, in part, by awards R01GM131398 from the National
Institutes of Health and NSF1942143 from the National Science Foundation.

## References

  * [1] Sanja Vickovic, Gökcen Eraslan, Fredrik Salmén, Johanna Klughammer, Linnea Stenbeck, Denis Schapiro, Tarmo Äijö, Richard Bonneau, Ludvig Bergenstråhle, José Fernandéz Navarro, Joshua Gould, Gabriel K Griffin, Åke Borg, Mostafa Ronaghi, Jonas Frisén, Joakim Lundeberg, Aviv Regev, and Patrik L Ståhl.  High-definition spatial transcriptomics for in situ tissue profiling.  Nat Methods, 16(10):987–990, September 2019.
  * [2] Samuel G Rodriques, Robert R Stickels, Aleksandrina Goeva, Carly A Martin, Evan Murray, Charles R Vanderburg, Joshua Welch, Linlin M Chen, Fei Chen, and Evan Z Macosko.  Slide-seq: A scalable technology for measuring genome-wide expression at high spatial resolution.  Science, 363(6434):1463–1467, March 2019.
  * [3] Robert R Stickels, Evan Murray, Pawan Kumar, Jilong Li, Jamie L Marshall, Daniela J Di Bella, Paola Arlotta, Evan Z Macosko, and Fei Chen.  Highly sensitive spatial transcriptomics at near-cellular resolution with Slide-seqV2.  Nat Biotechnol, 39(3):313–319, December 2020.
  * [4] Patrick AP Moran.  Notes on continuous stochastic phenomena.  Biometrika, 37(1/2):17–23, 1950.
  * [5] Robert C Geary.  The contiguity ratio and statistical mapping.  The incorporated statistician, 5(3):115–146, 1954.
  * [6] Valentine Svensson, Sarah A. Teichmann, and Oliver Stegle.  Spatialde: identification of spatially variable genes.  Nature Methods, 15(5):343–346, May 2018.
  * [7] Ilia Kats, Roser Vento-Tormo, and Oliver Stegle.  Spatialde2: Fast and localized variance component analysis of spatial transcriptomics.  bioRxiv, 2021.
  * [8] Shiquan Sun, Jiaqiang Zhu, and Xiang Zhou.  Statistical analysis of spatial expression patterns for spatially resolved transcriptomic studies.  Nature Methods, 17(2):193–200, Feb 2020.
  * [9] Lukas M. Weber, Arkajyoti Saha, Abhirup Datta, Kasper D. Hansen, and Stephanie C. Hicks.  nnsvg for the scalable identification of spatially variable genes using nearest-neighbor gaussian processes.  Nature Communications, 14(1):4059, Jul 2023.
  * [10] Qiwei Li, Minzhe Zhang, Yang Xie, and Guanghua Xiao.  Bayesian modeling of spatial molecular profiling data via Gaussian process.  Bioinformatics, 37(22):4129–4136, 06 2021.
  * [11] Daniel Edsgärd, Per Johnsson, and Rickard Sandberg.  Identification of spatial expression trends in single-cell gene expression data.  Nature Methods, 15(5):339–342, May 2018.
  * [12] Ke Zhang, Wanwan Feng, and Peng Wang.  Identification of spatially variable genes with graph cuts.  Nature Communications, 13(1):5488, Sep 2022.
  * [13] Alma Andersson and Joakim Lundeberg.  sepal: identifying transcript profiles with spatial patterns by diffusion-based modeling.  Bioinformatics, 37(17):2644–2650, 03 2021.
  * [14] Junjie Zhu and Chiara Sabatti.  Integrative spatial single-cell analysis with graph-based feature learning.  bioRxiv, 2020.
  * [15] Brendan F Miller, Dhananjay Bambah-Mukku, Catherine Dulac, Xiaowei Zhuang, and Jean Fan.  Characterizing spatial gene expression heterogeneity in spatially resolved single-cell transcriptomic data with nonuniform cellular densities.  Genome Res, 31(10):1843–1855, May 2021.
  * [16] Ruben Dries, Qian Zhu, Rui Dong, Chee-Huat Linus Eng, Huipeng Li, Kan Liu, Yuntian Fu, Tianxiao Zhao, Arpan Sarkar, Feng Bao, et al.  Giotto: a toolbox for integrative analysis and visualization of spatial expression data.  Genome biology, 22:1–31, 2021.
  * [17] Yuhan Hao, Stephanie Hao, Erica Andersen-Nissen, William M Mauck, Shiwei Zheng, Andrew Butler, Maddie J Lee, Aaron J Wilk, Charlotte Darby, Michael Zager, et al.  Integrated analysis of multimodal single-cell data.  Cell, 184(13):3573–3587, 2021.
  * [18] Jiaqiang Zhu, Shiquan Sun, and Xiang Zhou.  Spark-x: non-parametric modeling enables scalable and robust detection of spatial expression patterns for large spatial transcriptomic studies.  Genome Biology, 22(1):184, Jun 2021.
  * [19] Rui Jiang, Zhen Li, Yuhang Jia, Siyu Li, and Shengquan Chen.  SINFONIA: Scalable identification of spatially variable genes for deciphering spatial domains.  Cells, 12(4), February 2023.
  * [20] Nuha BinTayyash, Sokratia Georgaka, S T John, Sumon Ahmed, Alexis Boukouvalas, James Hensman, and Magnus Rattray.  Non-parametric modelling of temporal and spatial counts data from RNA-seq experiments.  Bioinformatics, 37(21):3788–3795, November 2021.
  * [21] Peiyao Zhao, Jiaqiang Zhu, Ying Ma, and Xiang Zhou.  Modeling zero inflation is not necessary for spatial transcriptomics.  Genome Biology, 23(1):118, May 2022.
  * [22] Xi Jiang, Guanghua Xiao, and Qiwei Li.  A bayesian modified ising model for identifying spatially variable genes from spatial transcriptomics data.  Statistics in Medicine, 41(23):4647–4665, 2022.
  * [23] Aaron TL Lun, Davis J McCarthy, and John C Marioni.  A step-by-step workflow for low-level analysis of single-cell rna-seq data with bioconductor.  F1000Research, 5, 2016.
  * [24] Davis J McCarthy, Kieran R Campbell, Aaron TL Lun, and Quin F Wills.  Scater: pre-processing, quality control, normalization and visualization of single-cell rna-seq data in r.  Bioinformatics, 33(8):1179–1186, 2017.
  * [25] Minsheng Hao, Kui Hua, and Xuegong Zhang.  SOMDE: a scalable method for identifying spatially variable genes with self-organizing map.  Bioinformatics, 37(23):4392–4398, 06 2021.
  * [26] Jinge Yu and Xiangyu Luo.  Identification of cell-type-specific spatially variable genes accounting for excess zeros.  Bioinformatics, 38(17):4135–4144, 07 2022.
  * [27] Xin Yuan, Yanran Ma, Ruitian Gao, Shuya Cui, Yifan Wang, Botao Fa, Shiyang Ma, Ting Wei, Shuangge Ma, and Zhangsheng Yu.  Heartsvg: a fast and accurate method for spatially variable gene identification in large-scale spatial transcriptomic data.  bioRxiv, 2023.
  * [28] Julien Moehlin, Bastien Mollet, Bruno Maria Colombo, and Marco Antonio Mendoza-Parra.  Inferring biologically relevant molecular tissue substructures by agglomerative clustering of digitized spatial transcriptomes with multilayer.  Cell Systems, 12(7):694–705.e3, 2021.
  * [29] Yaowu Liu, Sixing Chen, Zilin Li, Alanna C Morrison, Eric Boerwinkle, and Xihong Lin.  Acat: a fast and powerful p value combination method for rare-variant analysis in sequencing studies.  The American Journal of Human Genetics, 104(3):410–421, 2019.
  * [30] Natesh S. Pillai and Xiao-Li Meng.  An unexpected encounter with cauchy and lévy, 2015.
  * [31] Daowen Zhang and Xihong Lin.  Hypothesis testing in semiparametric additive mixed models.  Biostatistics, 4(1):57–74, January 2003.
  * [32] Faming Liang.  A double metropolis–hastings sampler for spatial models with intractable normalizing constants.  Journal of Statistical Computation and Simulation, 80:1007–1022, 09 2010.
  * [33] Juexin Wang, Jinpu Li, Skyler T. Kramer, Li Su, Yuzhou Chang, Chunhui Xu, Michael T. Eadon, Krzysztof Kiryluk, Qin Ma, and Dong Xu.  Dimension-agnostic and granularity-based spatially variable gene identification using bsp.  Nature Communications, 14(1):7367, Nov 2023.
  * [34] Lulu Shang and Xiang Zhou.  Spatially aware dimension reduction for spatial transcriptomics.  Nature Communications, 13(1):7203, Nov 2022.
  * [35] Yoav Benjamini and Yosef Hochberg.  Controlling the false discovery rate: A practical and powerful approach to multiple testing.  Journal of the Royal Statistical Society: Series B (Methodological), 57(1):289–300, 1995.
  * [36] John D Storey and Robert Tibshirani.  Statistical significance for genomewide studies.  Proceedings of the National Academy of Sciences, 100(16):9440–9445, 2003.
  * [37] Yoav Benjamini and Daniel Yekutieli.  The control of the false discovery rate in multiple testing under dependency.  The Annals of Statistics, 29(4):1165 – 1188, 2001.
  * [38] William N Venables and Brian D Ripley.  Modern applied statistics with S-PLUS.  Springer Science & Business Media, 2013.
  * [39] Jiaqiang Zhu, Lulu Shang, and Xiang Zhou.  Srtsim: spatial pattern preserving simulations for spatially resolved transcriptomics.  Genome Biology, 24(1):39, Mar 2023.
  * [40] Carissa Chen, Hani Jieun Kim, and Pengyi Yang.  Evaluating spatially variable gene detection methods for spatial transcriptomics data.  bioRxiv, 2022.
  * [41] Jian Hu, Xiangjie Li, Kyle Coleman, Amelia Schroeder, Nan Ma, David J. Irwin, Edward B. Lee, Russell T. Shinohara, and Mingyao Li.  Spagcn: Integrating gene expression, spatial location and histology to identify spatial domains and spatially variable genes by graph convolutional network.  Nature Methods, 18(11):1342–1351, Nov 2021.
  * [42] Chihao Zhang, Kangning Dong, Kazuyuki Aihara, Luonan Chen, and Shihua Zhang.  Stamarker: Determining spatial domain-specific variable genes with saliency maps in deep learning.  bioRxiv, 2022.
  * [43] Marc Elosua-Bayes, Paula Nieto, Elisabetta Mereu, Ivo Gut, and Holger Heyn.  SPOTlight: seeded NMF regression to deconvolute spatial transcriptomics spots with single-cell transcriptomes.  Nucleic Acids Research, 49(9):e50–e50, 02 2021.
  * [44] Dylan M Cable, Evan Murray, Luli S Zou, Aleksandrina Goeva, Evan Z Macosko, Fei Chen, and Rafael A Irizarry.  Robust decomposition of cell type mixtures in spatial transcriptomics.  Nat Biotechnol, 40(4):517–526, February 2021.
  * [45] Rui Dong and Guo-Cheng Yuan.  Spatialdwls: accurate deconvolution of spatial transcriptomic data.  Genome biology, 22(1):145, 2021.


# A generalized voter model on complex networks

Casey M. Schneider-Mizell, Leonard M. Sander  Department of Physics,
University of Michigan, Ann Arbor, MI 48105

###### Abstract

We study a generalization of the voter model on complex networks, focusing on
the scaling of mean exit time. Previous work has defined the voter model in
terms of an initially chosen node and a randomly chosen neighbor, which makes
it difficult to disentangle the effects of the stochastic process itself
relative to the network structure. We introduce a process with two steps, one
that selects a pair of interacting nodes and one that determines the direction
of interaction as a function of the degrees of the two nodes and a parameter
$\alpha$ which sets the likelihood of the higher degree node giving its state.
Traditional voter model behavior can be recovered within the model. We find
that on a complete bipartite network, the traditional voter model is the
fastest process. On a random network with power law degree distribution, we
observe two regimes. For modest values of $\alpha$, exit time is dominated by
diffusive drift of the system state, but as the high nodes become more
influential, the exit time becomes becomes dominated by frustration effects.
For certain selection processes and parameters values, an intermediate regime
occurs where exit occurs after exponential mixing.

###### pacs:

89.75.Fb, 02.50.Ey, 89.75.Hc

##  I Introduction

The voter model has been extensively studied on lattices Liggett (2005) and,
in recent years, on complex networks Suchecki et al. (2005a); Sood and Redner
(2005); Suchecki et al. (2005b); Castellano et al. (2005); Vazquez and Eguiluz
(2008) and is closely related to models of language evolution Baxter et al.
(2006), ecological dynamics Hubbell (2001), opinion dynamics Castellano et al.
(2007), and epidemic spread Pastor-Satorras and Vespignani (2001). The voter
model defines a dynamical process where nodes are each assigned one of two
states, $+1$ or $-1$. Connections are defined on a lattice by nearest
neighbors or on an arbitrary network by edges. Each update step consists of
selecting a pair of nodes and giving the state of one node to the other. In
the most frequently studied version of the model, the first node chosen adopts
the state of the second. The most interesting object of study is the mean exit
time, i.e. the mean time to achieve complete agreement. For the voter model on
complex networks, a node will be chosen second with frequency proportional to
its degree, and so its influence is fixed by the selection process. In this
paper, we introduce a generalized voter model with a single tunable parameter
that allows control of the influence of topology in a manner independent of
the selection process. In this generalized model, the probability of a node
giving its state to its neighbor is proportional to $k^{\alpha}$, where $k$ is
the node’s degree and $\alpha$ is chosen.

Voter model processes based entirely on selection frequency are denoted either
_link update_ or _node update_ Suchecki et al. (2005b). In link update
dynamics, every time step an edge is selected uniformly at random. One of the
two nodes at the ends of the link is then chosen randomly to give its state to
the other. In node update dynamics (described above) a node is selected to
adopt the state of a random neighbor. The neighbor chosen at random is likely
to have high degree, so high degree nodes have more influence.

##  II Generalized voter model

The two processes described above specify both the selection of a pair of
interacting nodes and which node adopts the state of the other. We separate
this process into two distinct steps to better understand the contribution of
each aspect of the process.

Given a network of size $N$, each node $i$ has state $s_{i}=\pm 1$ and degree
$k_{i}$. We define $P_{ij}$ to be the probability of giving to node $i$ the
state of node $j$ during a given time step. There are two independent
components of this event: the probability $S_{ij}$ of selecting an edge
connecting nodes $i$ and $j$ and the probability $W(k_{i},k_{j})$ that a node
with degree $k_{j}$ gives its state to a node with degree $k_{i}$. Thus,

$P_{ij}=S_{ij}W(k_{i},k_{j}).$ (1)

The form of $W(k_{i},k_{j})$ is motivated by comparison with node update
dynamics in uncorrelated networks. In uncorrelated networks, node update
dynamics can be described by considering all nodes of like degree to be
indistinguishable in the ensemble average Sood and Redner (2005). In the
following, $\\{i\\}$ and $\\{j\\}$ refer to the subgroups of all nodes with
degrees $k_{i}$ and $k_{j}$, respectively. Under the node update process, the
probability of giving a node with degree $k_{i}$ the state of a node with
degree $k_{j}$ is

$P_{ij}=n_{i}\frac{k_{j}n_{j}}{\mu_{1}}$ (2)

where $n_{i}$, $n_{j}$ are the fraction of nodes with degrees $k_{i}$ and
$k_{j}$ and $\mu_{1}=\sum_{i}n_{i}k_{i}$ is the average degree. This can be
interpreted as the probability $n_{i}$ of selecting a node with degree $k_{i}$
times the probability of following an edge into a node with degree $k_{j}$. A
node in $\\{j\\}$ has $k_{j}$ edges, so the probability of following an edge
into it is proportional to $k_{j}$. If node pairs are selected in this manner,
a particular pair of nodes in $\\{i\\}$ and $\\{j\\}$ can be chosen by either
picking from $\\{i\\}$ and following an edge to a node in $\\{j\\}$ or by
picking from $\\{j\\}$ and following an edge to a node in $\\{i\\}$. For node
selection, this gives

$S^{n}_{ij}=\frac{n_{i}n_{j}}{2\mu_{1}}(k_{i}+k_{j}).$ (3)

We propose a generalization that includes the standard node update dynamics,
which requires $S_{ij}W(k_{i},k_{j})=P_{ij}$. The only form of
$W(k_{i},k_{j})$ to do this is

$W(k_{i},k_{j})=\frac{k_{j}}{k_{i}+k_{j}}.$ (4)

This form is also consistent with our definition of $W(k_{i},k_{j})$ as a
probability. Notably, it suggests a one parameter generalization:

$W(k_{i},k_{j},\alpha)=\frac{k_{j}^{\alpha}}{k_{i}^{\alpha}+k_{j}^{\alpha}}.$
(5)

Qualitatively, the parameter $\alpha$ determines how much a node asserts its
degree when transmitting its state. For $\alpha>0$, the higher degree node of
a pair is more likely to give its state to the lower degree node, a bias that
increases with $\alpha$. For $\alpha<0$, the opposite is true. The special
case $\alpha=0$ ignores topology in determining the direction of interaction
since $W(k_{i},k_{j},0)=1/2$ always. Edge update dynamics is recovered by
using edge selection to find pairs and setting $\alpha=0$. Node update
dynamics occurs when node selection determines pairs and $\alpha=1$. A
recently investigated “invasion” dynamic, where a node is picked to give its
state to a random neighbor (opposite the traditional model), occurs for node
selection and $\alpha=-1$. If all nodes of the network have the same degree,
as in a mean field or lattice topology, then all values of $\alpha$ are
equivalent to the traditional voter model.

Effectively, this model assumes some connection between the behavior of the
agents and the underlying network on which they live. For example, if this
were to be thought of as a model of opinion dynamics, a value of $\alpha>1$
under node update selection would correspond to a situation where an
individual prefers to behave like those who are more connected than himself.
Celebrities or well-regarded experts are extremely influential, for example,
but the same forces that drive their high visibility also keep them from being
influenced by the non-famous. The forces of influence and accessibility
compete, such that a small value of $\alpha$ makes all nodes able to change
state quickly, but limits the influence of any node. A high value of $\alpha$
makes high degree nodes influential, thus able to order their neighborhood
quickly, but those influential nodes will flip only on rare occasions.

###  II.1 Dynamics

To understand the dynamics, we study the master equation for an arbitrary
network. The probability of a system being in state $\mathbf{s}=\\{s_{i}\\}$
at time $t$ is defined to be $P(\mathbf{s},t)$. Denote by $\mathbf{s}^{i}$ the
state $\mathbf{s}$ where $s_{i}\mapsto-s_{i}$ and let $S_{ij}$ be the
probability of selecting the edge between nodes $i$ and $j$. For brevity, we
write $W(k_{i},k_{j},\alpha)=W_{ij}$. The master equation is

$\frac{d}{dt}{P(\mathbf{s},t)}=\sum_{ij}S_{ij}W_{ij}\left(\frac{1+s_{i}s_{j}}{2}\right)P(\mathbf{s}^{i},t)-S_{ij}W_{ij}\left(\frac{1-s_{i}s_{j}}{2}\right)P(\mathbf{s},t)$
(6)

Let $\rho_{k}$ be the ensemble average probability of a node with degree $k$
being in a $+1$ state. Averaged over all random graphs, nodes of like degree
are indistinguishable and we consider each such like-degree subgraph
separately Sood and Redner (2005). The evolution of arbitrary ensemble average
functions can be found in a straight-forward manner Krapivsky (1992), giving:

$\frac{d\rho_{k}}{dt}=\sum_{l}\frac{S_{kl}}{n_{k}}W_{kl}(\rho_{l}-\rho_{k}).$
(7)

We can find a conserved magnetization, $\rho^{*}$, by choosing coefficients
$C_{i}$ such that

$\frac{d\rho^{*}}{dt}=\sum_{k}C_{k}\frac{d}{dt}\rho_{k}=0.$ (8)

Since $S_{ij}$ is symmetric in $i$ and $j$, this can happen for arbitrary
subgraph densities only if

$C_{k}\frac{W_{kl}}{n_{k}}=C_{l}\frac{W_{lk}}{n_{l}}.$ (9)

This implies that $C_{k}\propto n_{k}k^{\alpha}$. Normalizing,

$\rho^{*}=\frac{\sum_{k}n_{k}k^{\alpha}\rho_{k}}{\sum_{k}n_{k}k^{\alpha}}=\frac{1}{\mu_{\alpha}}\sum_{k}n_{k}k^{\alpha}\rho_{k}$
(10)

where $\mu_{\alpha}$ is the $\alpha^{th}$ moment of the degree distribution.
Note that the ensemble conserved magnetism is independent of the process of
selecting node pairs.

##  III Bipartite Network

The simplest degree homogeneous topology is the fully connected bipartite
network. Such a network is given by two groups of nodes, group $A$ with size
$a$ and group $B$ with size $b$. A node in group $A$ is connected to every
node in group $B$, but none in group $A$. The degree of nodes in $A$, $k_{A}$,
is the size of $B$, giving $k_{A}=b$ and similarly $k_{B}=a$. In this
situation all edges are interchangeable, so there is no difference between the
two selection processes. We will consider only the effect of $\alpha$. Let
$\rho_{a}$ be the concentration of $+1$ opinions in $A$ and $\rho_{b}$ be the
concentration of $+1$ opinions in $B$. In our model, the special value
$\alpha=1$ is equivalent to the case studied in Sood and Redner (2005) on the
same network and we follow a similar procedure, omitting details that can be
found there. From Equation (10), the conserved magnetization is

$\rho^{*}=\frac{1}{a^{\alpha-1}+b^{\alpha-1}}\left(b^{\alpha-1}\rho_{a}+a^{\alpha-1}\rho_{b}\right)$
(11)

For any initial conditions, the ensemble average subgraph densities approaches
$\rho^{*}$. If all nodes in $A$ start as $+1$ and all nodes in $B$ start as
$-1$, then the probability of ending in the $+1$ state is

$P_{+}=\frac{b^{\alpha-1}}{b^{\alpha-1}+a^{\alpha-1}}$ (12)

The mean exit time $T_{\alpha}$ is given by the backward Komologorov equation
Redner (2001). $T_{\alpha}$ solves

$-\frac{a^{\alpha}+b^{\alpha}}{a+b}=(\rho_{a}-\rho_{b})(b^{\alpha-1}\partial_{b}-a^{\alpha-1}\partial_{a})T+\\\
\frac{1}{2}(\rho_{a}+\rho_{b}-2\rho_{a}\rho_{b})(a^{\alpha-2}\partial^{2}_{a}+b^{\alpha-2}\partial^{2}_{b})T.$
(13)

where $\partial_{a}$ and $\partial_{b}$ are partial derivatives with respect
to the initial subgraph densities.The first term describes convection, which
brings the subgraph densities to some equal value, and the second term
describes the diffusion of the network-wide state Sood and Redner (2005). The
convective dynamics can be shown to be fast for all $\alpha$.

The fast step toward equal subgraph densities has a negligible impact on
extinction time and we can consider only the subsequent one dimensional
problem. We define $\rho=\rho_{a}=\rho_{b}$ and apply a change of variables
using Equation (11). After integrating,

$T=-(a^{1-\alpha}+b^{1-\alpha})(a^{\alpha-1}+b^{\alpha-1})\frac{ab}{a+b}(\rho\log(\rho)+(1-\rho)\log(1-\rho)).$
(14)

This has a similar form to the standard voter model, but with a factor that is
symmetric about $\alpha=1$ and non-vanishing. If we take $a=\lambda b$, then

$T\propto(2+\lambda^{1-\alpha}+\lambda^{\alpha-1})\frac{\lambda}{1+\lambda}b.$
(15)

If $\lambda\gg 1$, corresponding to a star-like graph,

$T\sim\lambda^{|\alpha-1|}b.$ (16)

This scaling is confirmed in simulations (see Figure 1). Notably, the standard
voter model, $\alpha=1$, is the fastest process for any complete bipartite
network.

##  IV Arbitrary networks

Similar analysis extends naturally to networks in which a node’s degree
determines the network structure. Many random network models fall into this
category, notably any random network generated by the configuration model,
including those with scale-free distributions, and Erdos-Renyi networks Newman
(2003). Small world networks are not included, however, as certain nodes have
exceptional topological characteristics that are independent of their degree
Watts and Strogatz (1998).

As in Equation (13), we can write the equation satisfied by the mean exit time
on an arbitrary network. In the following, $\rho_{i}$ refers to the density of
$+1$ states on the subgraph of nodes with degree $k_{i}=i$.

$\delta_{t}=\sum_{ij}S_{ij}W^{\alpha}_{ij}(\rho_{j}-\rho_{i})\delta_{i}\partial_{i}T+\frac{1}{2}\sum_{ij}S_{ij}W^{\alpha}_{ij}(\rho_{i}+\rho_{j}-2\rho_{i}\rho_{j})\delta_{i}^{2}\partial_{i}^{2}T$
(17)

The system is again split into a convective term and a diffusive term. The
assumption of fast approach to well-mixed state must be treated more carefully
in our generalized model, but there do exist cases where diffusion dominates.
Namely, it was observed in Sood and Redner (2005) that node update dynamics on
a scale free network has fast convection compared to its diffusive exit time.
There is, in general, a prefactor to the term $\rho_{i}(1-\rho_{j})$ when
evaluating the probability that two nodes have differing states Vazquez and
Eguiluz (2008). Numerical simulations show that this is of order one and not
generally a constant over time when $\alpha\neq 1$, so it does not affect the
final scaling results.

We expect this to not be true for all $\alpha$. In the case of $\alpha\gg 1$,
a node with degree higher than all its neighbors will act to dictate its
neighbors’ states, but only rarely be changed itself. The network in this case
may not be able to quickly approach the global equilibrium given by
$\rho^{*}$, since these locally highest degree nodes will be pinned for a time
dependent on $\alpha$. If this duration is longer than the time for the rest
of the system to become ordered via mixing and diffusion, a quasi-frustrated
state occurs. The exact local topology, rather than just degree distributions,
can dominate the dynamics.

Let us suppose that the exit time is diffusion dominated and will return to
discuss the validity of this assumption. The system can be approximated by a
one dimensional equation in
$\rho=\frac{1}{\mu_{\alpha}}\sum_{i}n_{i}k^{\alpha}_{i}\rho_{i}$:

$-N=\frac{1}{\mu^{2}_{\alpha}}\left(\sum_{ij}S_{ij}W^{\alpha}_{ij}k_{i}^{2\alpha}\right)\rho(1-\rho)\partial_{\rho}^{2}T.$
(18)

And thus

$T\propto-\frac{N\mu^{2}_{\alpha}}{\sum_{ij}S_{ij}W_{ij}^{\alpha}k_{i}^{2\alpha}}.$
(19)

The denominator can be simplified by noting that

$\sum_{ij}S_{ij}W^{\alpha}_{ij}k_{i}^{2\alpha}=\frac{1}{2}\left(\sum_{ij}S_{ij}W^{\alpha}_{ij}k_{i}^{2\alpha}+\sum_{ij}S_{ij}W^{\alpha}_{ji}k_{j}^{2\alpha}\right)=\frac{1}{2}\sum_{ij}S_{ij}k_{i}^{\alpha}k_{j}^{\alpha}.$
(20)

This gives:

$T\propto
N\frac{\mu_{\alpha}^{2}}{\sum_{ij}S_{ij}k_{i}^{\alpha}k_{j}^{\alpha}}.$ (21)

Since $S_{ij}$ is a probability, the sum can be thought of as a weighted
average over selection probabilities. Interestingly, for $\alpha=0$, neither
the form of interaction selection nor the network topology matter. In that
case, $\mu_{0}=1$ and $W_{ij}=1/2$, so

$T_{\alpha=0}\propto N.$ (22)

This agrees with the observation in Suchecki et al. (2005a) that exit times
scale with $N$ in situations where the unweighted magnetization is conserved,
which corresponds exactly with $\alpha=0$.

To go farther, we need to specify the selection scheme and the network. We
focus our consideration on random uncorrelated scale-free networks with degree
distribution $n_{k}\sim k^{-\nu}$. Networks with power law distributions
appear in a variety of social and biological contexts and exhibit a range of
interesting behaviors Newman (2003). Let us first consider node update, for
which $S_{ij}=n_{i}n_{j}\frac{k_{i}+k_{j}}{2\mu_{1}}$. Then Equation (21)
becomes:

$T_{\mathcal{N}}\propto N\frac{\mu_{1}\mu_{\alpha}}{\mu_{\alpha+1}}.$ (23)

The $\alpha^{th}$ moment can be approximated by an integral:

$\mu_{\alpha}\sim\int^{k_{max}}k^{\alpha}n(k)dk$ (24)

up to an effective maximum degree $k_{max}$, defined by
$\int_{k_{max}}^{\infty}n(k)dk=1/N$ Krapivsky and Redner (2002). It is easily
seen that $k_{max}\sim N^{1/(\nu-1)}$.

$T_{\mathcal{N}}\propto\left\\{\begin{array}[]{ll}N^{\frac{\nu-2}{\nu-1}}&\alpha>\nu-1\\\
N^{\frac{2\nu-\alpha-3}{\nu-1}}&\nu-2<\alpha<\nu-1\\\ N&\alpha<\nu-2\\\
\end{array}\right.$ (25)

For $\nu>2$ and any $\alpha$, the exit time increases without bound as system
size increases. We simulated the process on random network generated by the
configuration model Molloy and Reed (1995) and found good agreement with our
predictions (see Figure 2).

For edge update dynamics, $S_{ij}=n_{i}n_{j}\frac{k_{i}k_{j}}{\mu^{2}}.$ Low
degree nodes are selected less frequently under edge selection than node
selection. The diffusive exit time can be calculated similarly, giving:

$T_{\mathcal{E}}=N\left(\frac{\mu_{1}\mu_{\alpha}}{\mu_{\alpha+1}}\right)^{2}.$
(26)

The approximate scaling for edge update is

$T_{\mathcal{E}}\propto\left\\{\begin{array}[]{ll}N^{\frac{\nu-3}{\nu-1}}&\alpha>\nu-1\\\
N^{\frac{3\nu-2\alpha-5}{\nu-1}}&\nu-2<\alpha<\nu-1\\\ N&\alpha<\nu-2\\\
\end{array}\right.$ (27)

This leads to very different scaling behavior. For the parameter regions

$\displaystyle\nu<3$ $\displaystyle\quad\quad\alpha>\nu-1$ (28)
$\displaystyle\nu<\frac{2\alpha+5}{3}$
$\displaystyle\quad\quad\nu-2<\alpha<\nu-1$ (29)

the diffusive exit time vanishes as $N$ increases. Simulations show that there
is a diffusive region which agrees with our predictions for smaller values of
$\alpha$ (see Figure 3).

The convective process involves an exponential decay of each $\rho_{i}$ to its
equilibrium value, with rate determined by the network structure. When this
value is comparable to $1/n_{i}N$ (or $1-1/n_{i}N$) for all $i$, the system
effectively reaches convergence before diffusive time scales matter. We can
read off the dynamics from Equation 7. The decay rate $\tau$ is not strongly
dependent on system size, so the relationship with $N$ comes only from the
exit condition

$\frac{1}{N}\sim e^{-T_{\mathcal{E}}/\tau}.$ (30)

This implies that when diffusive time scales vanish,

$T_{\mathcal{E}}\sim\log(N).$ (31)

Exponential mixing is only observed as a dominant time scale under link update
selection for values of $N$, $\nu$, and $\alpha$ such that frustration does
not occur.

To understand where frustration comes from, consider a node whose degree is
higher than that of all of its neighbors. In the limit
$\alpha\rightarrow\infty$, such a node will take an arbitrarily long time to
flip. For large values of $\alpha$, escape from this frustrated state must
then be the dominant time scale, as diffusive time scales for high $\alpha$
are dependent only on $\nu$. For either selection process, the probability of
a given node being sampled in a time step is a constant, independent of system
size. It is also possible for small connected subgraphs to be frustrated if
every edge out of the subgraph is to a node of lower degree. As a result of
locally highest degree nodes and clusters, a small number of hard-to-flip
nodes can hold back a large network from reaching convergence (Figure 4).

The continuum treatment in Equation 17 averages over network ensembles before
solving for the exit time. In the high $\alpha$ case, considering an averaged
network ignores local leader effects and fails to give an accurate solution.
The approach to an ensemble equilibrium state does not occur, so diffusion
about this state is not a valid assumption (Figure 5).

The degree distribution of locally highest degree nodes can be approximated
for the non-assortative case quite simply. The degree distribution of local
leaders, $p_{ll}(k)$, is the independent product that a node has degree $k$
and that all $k$ neighbors have a degree less than $k$ Blondel et al. (2007):

$p_{ll}(k)=p(k)\left(\sum_{k^{\prime}<k}\frac{k^{\prime}p(k^{\prime})}{\mu_{1}}\right)^{k}.$
(32)

No term is strongly related to system size, so the total number of local
leaders $N_{ll}$ scales linearly with $N$.

The dynamics of these nodes are based on extremely local behavior and thus
very hard to approximate, but we can construct the slowest possible node and
understand its behavior. For a sufficiently large system, one such node will
likely exist and thus dominate the exit time. The probability on any given
time step to flip a local leader with degree $l$ is formally given by:

$P_{flip}(l)=\sum_{i}S_{il}\frac{k_{i}^{\alpha}}{k_{i}^{\alpha}+l^{\alpha}}P(\sigma_{i}\neq\sigma_{l})$
(33)

where $P(\sigma_{i}\neq\sigma_{l})$ is the probability that the neighboring
state differs from the local leader’s state and subscripts index individual
nodes.

For either selection process, when $l$ is large,

$S_{li}\sim\frac{1}{N}.$ (34)

We assume the probability of a neighbor being in the opposite state is finite
and treat the following as a lower bound. Simulations show that this is
reasonable. Assuming that $l^{\alpha}\gg k_{i}^{\alpha}$,

$P(l)\approx\frac{1}{Nl^{\alpha-1}}\left(\frac{1}{l}\sum_{i}k_{i}^{\alpha}\right).$
(35)

The slowest situation occurs when the local leader in question is surrounded
by relatively low degree nodes and $P(l)\sim l^{1-\alpha}/N$. For this case,
the time to flip in units of system size scales as:

$T_{f}\sim l^{\alpha-1}.$ (36)

The highest degree of the local leaders scales linearly with the global
highest degree Blondel et al. (2007), so for $\nu<3$ we have $T_{f}\sim
N^{(\alpha-1)/(\nu-1)}$. This will be true for large $\alpha$, but requires
the system size to be high enough that at least one rare node such as the one
described is likely. Such networks are prohibitively large to sample, but the
qualitative situation described explains why frustration sets in at a smaller
value of $\alpha$ for $\nu=2.4$ than $\nu=2.8$ in Figure 2. Local leaders for
small $\nu$ have a higher degree than for large $\nu$, though numerical
comparison from Equation 32 shows them to be rarer.

##  V Summary

Recent work Baxter et al. (2008) has found an approximate mean exit time for a
duplication process on networks with arbitrary edge weights, assuming that
diffusion is the dominant time scale. In this work, we demonstrate that there
are at least two natural ways for this estimate of exit time to fail. As
observed in Baxter _et al._ , the time for the system to reach a metastable
equilibrium can be at least as large as the diffusive exit time scale. We see
this in the edge selection process for values of $\alpha$ and $\nu$ where the
diffusive exit time vanishes as system size gets large. The frustrated
dynamics in the node selection mode, however, presents a new way in which the
diffusive estimate can fail. System dynamics are driven by a small number of
topologically special nodes, breaking the assumption that a continuum
description applies.

We have defined and analyzed a single parameter voter-model-like stochastic
process that is identical to the original voter model on a homogeneous
network. On a complex network, our generalized voter model has a tunable
dependence on local network topology, allowing us to control the differing
effects of selection and the direction of influence. On complete bipartite
graphs, the traditional voter model is the fastest process to reach an
absorbing state. On scale free networks, the dynamics depend strongly on the
selection process. Node selection has two regimes; a diffusive one
characterized by a well defined average magnetism and diffusion constant based
on global network properties, and one with frustrated dynamics stemming from
the neighborhood around a small number of locally highest degree nodes. Edge
selection, previously considered uninteresting, has three regimes. In addition
to diffusive and frustrated regimes, it also has a middle ground characterized
by exponential mixing. Understanding the dynamics involved in this transition
to frustration would be an interesting avenue for future work.

##  VI Acknowledgements

LMS is supported in part by NSF grant DMS 0554587. CMS would like to thank A.
Stein, S. Cobey, and D. Adams for useful conversations.

Figure 1: Circles are the simulated fit of $d\log{T}/d\log{N}$ for a complete
bipartite network with $M=40$ and $N$ ranging from 100 to 5000. The dotted
line is the scaling of Equation 16 for the same range of $N$.

Figure 2: Simulated values of $d\log{T}/d\log{N}$ under _node_ selection based
on several hundred runs for $N$ from 750–15000. Circles correspond to
$\nu=2.8$, squares to $\nu=2.4$. The dashed line is the scaling based on the
diffusive estimate, calculated by fitting the numerically calculated sum in
Equation 21 for similar values of $N$. Note that frustration begins to
dominate for $\nu=2.4$ at $\alpha>1.6$, causing the deviation from the
diffusive estimate.

Figure 3: Simulated values of $d\log{T}/d\log{N}$ under _edge_ selection from
several hundred runs of values of $N$ from 750–15000. Circles are for
$\nu=2.8$, squares for $\nu=2.4$. The dashed line comes from fitting the
numerically calculated sum in Equation 21 for similar values of $N$. The
horizontal line is the effective slope of $T\sim\log{N}$ for the system sizes
used.

Figure 4: A example of node update dynamics on a small network containing a
global highest degree node ($G$) and a separate local highest degree node
($L$) for $\alpha=1$ and $\alpha=10$. The darkness of a node corresponds to
the average fraction of time spent in a state opposite the final state of the
network over 1000 realizations with identical initial conditions. A darker
node has spent more time in a contrary state than a light node. For
$\alpha=1$, states are well mixed and $T=9.4$. For $\alpha=10$, mixing does
not occur. The local leader and its neighborhood spend most of the time in a
contrary state and $T=4814$.

Figure 5: Typical dynamics of the fraction of active edges (edges connecting
nodes with different states) for diffusion and frustration. Both have
$N=10000$, $\nu=2.8$, and edge selection. The top line, in black, is for
$\alpha=1$ and the bottom line, in gray, is for $\alpha=8$. The $\alpha=1$
case is diffusive and fluctuates to convergence after reaching the ensemble
average value. The $\alpha=8$ case decays exponentially to a value greater
than the ensemble average value, because a locally highest degree node or
cluster is slow to flip. Just before $t=1200$, these nodes flip and the system
reaches the absorbing state. Inset shows early time dynamics on a similar
network for $\alpha=10$ in gray and $\alpha=4$ in black. Both decay
exponentially initially, but the high $\alpha$ case becomes frustrated and the
other continues to convergence.

## References

  * Liggett (2005) T. M. Liggett, _Interacting Particle Systems_ (Springer Verlag, New York, New York, 2005).
  * Suchecki et al. (2005a) K. Suchecki, V. M. Eguiluz, and M. S. Miguel, Physical Review E 72 (2005a).
  * Sood and Redner (2005) V. Sood and S. Redner, Phys Rev Lett (2005).
  * Suchecki et al. (2005b) K. Suchecki, V. M. Eguiluz, and M. S. Miguel, Europhysics Letters (2005b).
  * Castellano et al. (2005) C. Castellano, V. Loreto, A. Barrat, and F. Cecconi, Physical Review E (2005).
  * Vazquez and Eguiluz (2008) F. Vazquez and V. M. Eguiluz (2008), arXiv:0803.1686v1.
  * Baxter et al. (2006) G. J. Baxter, R. A. Blythe, W. Croft, and A. J. McKane, Physical Review E (2006).
  * Hubbell (2001) S. P. Hubbell, _The Unified Neutral Theory of Biodiversity and Biogeography_ (Princeton University Press, Princeton, NJ, 2001).
  * Castellano et al. (2007) C. Castellano, S. Fortunato, and V. Loreto (2007), arXiv:0710.3256v1.
  * Pastor-Satorras and Vespignani (2001) R. Pastor-Satorras and A. Vespignani, Phys Rev Lett (2001).
  * Krapivsky (1992) P. L. Krapivsky, Physical Review A 45 (1992).
  * Redner (2001) S. Redner, _A Guide to First-Passage Processes_ (Cambridge University Press, Cambridge, UK, 2001).
  * Newman (2003) M. E. J. Newman, SIAM Review 45, 167 (2003).
  * Watts and Strogatz (1998) D. Watts and S. Strogatz, Nature (1998).
  * Krapivsky and Redner (2002) P. L. Krapivsky and S. Redner, Journal of Physics A Mathematical and General (2002).
  * Molloy and Reed (1995) M. Molloy and B. Reed, Random Structures and Algorithms (1995).
  * Blondel et al. (2007) V. D. Blondel, J.-L. Guillaume, J. Hendrickx, C. de Kerchove, and R. Lambiotte (2007), arXiv:0707.4064.
  * Baxter et al. (2008) G. J. Baxter, R. A. Blythe, and A. J. McKane (2008), arXiv:0801.3083.