††thanks: Email<EMAIL_ADDRESS>Email:
<EMAIL_ADDRESS>
# Robustness and Stability of Spin Glass Ground States to Perturbed
Interactions

Vaibhav Mohanty  Rudolf Peierls Centre for Theoretical Physics, University of
Oxford, Oxford, OX1 3NP, UK  MD-PhD Program and Program in Health Sciences and
Technology, Harvard Medical School, Boston, MA 02125, USA
and Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Ard A.
Louis  Rudolf Peierls Centre for Theoretical Physics, University of Oxford,
Oxford, OX1 3NP, UK

###### Abstract

Across many scientific and engineering disciplines, it is important to
consider how much the output of a given system changes due to perturbations of
the input. Here, we investigate the glassy phase of $\pm J$ spin glasses at
zero temperature by calculating the robustness of the ground states to flips
in the sign of single interactions. For random graphs and the Sherrington-
Kirkpatrick model, we find relatively large sets of bond configurations that
generate the same ground state. These sets can themselves be analyzed as
subgraphs of the interaction domain, and we compute many of their topological
properties. In particular, we find that the robustness, equivalent to the
average degree, of these subgraphs is much higher than one would expect from a
random model. Most notably, it scales in the same logarithmic way with the
size of the subgraph as has been found in genotype-phenotype maps for RNA
secondary structure folding, protein quaternary structure, gene regulatory
networks, as well as for models for genetic programming. The similarity
between these disparate systems suggests that this scaling may have a more
universal origin.

##  I Introduction

Systems in which the input can be represented as a sequence of characters
appear across science and engineering almost universally and especially
commonly in biology and computer science. In these fields, notions of system
robustness or sensitivity can be defined to quantify the outputs’ average
resistance to small changes in the input sequences.

In biology, the mapping from genotypes (which store the information) to
phenotypes (which describe biological properties) can be abstracted as
genotype-phenotype maps (GP maps). Examples include 4 letter RNA sequences and
20 letter protein sequences that can be mapped to their physical folded
states, and gene-regulatory networks, which can, for example, be described by
Boolean networks [1] where a set of weights represent the gene interaction
strengths. The responses of such systems to changes in the input sequences
have been extensively studied computationally and analytically [2, 3, 4, 5, 6,
7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22].

For GP maps, an important concept is the set of genotypes (sequences) that map
to a particular phenotype, often called the neutral set. These present a
number of commonalities across GP maps [12, 13, 18]: the neutral sets are
typically highly connected so that they can be viewed as networks which can be
traversed by single mutational steps, leading to enhanced evolvability [4,
19]. Neutral sets are also called neutral networks. In some cases, the neutral
network can split into smaller component networks which are disconnected due
to biophysical constraints [5, 20]. The neutral set size (NSS) can vary over
many orders of magnitude, and is typically strongly biased, with a small
fraction of the phenotypes taking up the majority of genotypes. Such phenotype
bias can strongly affect evolutionary outcomes [8, 11, 21, 23].

A key property of the neutral set for GP maps is the mutational robustness,
typically defined in the literature as the fraction $\rho_{p}$ of single-
character mutations in a genotype that produce the same phenotype $p$,
averaged over the neutral set of all genotypes that produce phenotype $p$.
That is to say, for a biological input-output map $f(g)$ which takes in a
genotype $g$ of $d$ characters chosen from an alphabet
$K=\\{K_{0},\dots,K_{k-1}\\}$, the robustness of phenotype $p$ is defined as

$\displaystyle\rho_{p}(f)=\frac{1}{\absolutevalue{\mathcal{G}_{p}}d(k-1)}\sum_{g\in\mathcal{G}_{p}}n_{p,g}$
(1)

where $\mathcal{G}_{p}$ is the neutral set of all genotypes $g$ whose output
is the phenotype $p$, and $n_{p,g}$ is the number of nearest-neighbors of $g$
mapping to $p$ defined as the number of genotypes $g^{\prime}$ satisfying
$f(g^{\prime})=f(g)=p$ that differ from $g$ by a Hamming distance of 1. Thus,
$\rho_{p}(f)\in[0,1]$ measures the mean probability that a mutation from
$g\in\mathcal{G}_{p}$ to a neighboring genotype $g^{\prime}\in\mathcal{G}_{p}$
results in the same phenotype $p$.

What all the biological models mentioned above hold in common is that they
each have high levels of mean robustness to changes in their inputs, excepting
pathological or adversarial examples. To quantify this statement, consider the
naive uncorrelated expectation that
$\rho_{p}\approx\absolutevalue{\mathcal{G}_{p}}/k^{d}$ because, in a large
randomly-assigned GP map, the probability that a nearest-neighboring genotype
yields phenotype $p$ is approximately equal to the probability that any
genotype drawn at random from the entire input space yields phenotype $p$.
This probability is $\absolutevalue{\mathcal{G}_{p}}/k^{d}$, so
$n_{p,g}\approx d(k-1)\absolutevalue{\mathcal{G}_{p}}/k^{d}$, and
$\rho_{p}\approx\absolutevalue{\mathcal{G}_{p}}/k^{d}$. And indeed this is
true for completely uncorrelated and randomly-assigned GP maps [12]. However,
biology and computer science-inspired GP maps such as the RNA secondary
structure maps, protein folding maps, gene regulatory networks, linear genetic
programs, and digital logic gate maps exhibit a remarkably similar observed
scaling law for robustness
$\rho_{p}\sim\mathcal{O}(\log\absolutevalue{\mathcal{G}_{p}})$ rather than
$\rho_{p}\sim\mathcal{O}(\absolutevalue{\mathcal{G}_{p}})$ [18]. Such enhanced
mutational robustness is necessary for the neutral sets to percolate [12]
which, in turn, is critical for the evolutionary process because it allows for
neutral exploration of the neutral set, greatly enhancing the ability of a
population to find novel phenotypic variation [4]. This enhanced robustness
also facilitates the exploration of fitness landscapes [24].

Robustness has a direct counterpart in computer science, the sensitivity,
which measures the likelihood that flipping a single input bit will alter the
output bit of Boolean functions $f:\\{0,1\\}^{d}\rightarrow\\{0,1\\}$ that map
binary sequences of length $d$ onto a single binary output. In other words,
low sensitivity corresponds to high robustness. Huang’s [25] famously short
solution to the decades-old Sensitivity Conjecture [26] concerning induced
subgraphs of the $n$-dimensional hypercube graph has recently brought a great
deal of attention to the sensitivity analysis of Boolean functions. While
scaling laws for sensitivity are typically not measured in the manner done for
biological robustness, there are many mathematical equivalencies between the
two. In the literature (see e.g., [26, 25, 27, 28]), quantities are defined
such as local sensitivity, function sensitivity, and average sensitivity,
which are similar to biological robustness. In this paper, we do not
explicitly use these definitions from sensitivity analysis, instead opting for
the GP map inspired ones. But, we point out that there is a seemingly
understudied connection between robustness in biological systems and
sensitivity in Boolean functions.

Given the wide range of systems for which high robustness is observed, we ask
here whether a similar phenomenon can be found in spin glasses, which have a
rich history in statistical and condensed matter physics. They have been
intensely studied since the 1970s [29, 30], and have led to many important
insights in physics and other related disciplines, including computer science
[31, 32, 33]. More recently, the spin glass Hamiltonian has been used as a
phenomenological model for epistatic genotype-to-fitness landscapes in which
different sites (e.g. DNA, genes expressions, or amino acids) may couple to
each other [34, 35, 36, 37, 38, 39, 40]. An important application has been to
viruses [35, 36, 37, 38, 39]. By taking sequence data over time, inverse
statistical physics methods can be employed to “learn” the interactions
(couplings/bond configurations) between different sites. In the context of
evolution, therefore, one can interpret the ground state of the spin glass
energy landscape as the global fitness peak on an evolutionary fitness
landscape. The interactions between spins in such a system can depend on a
number of biological or environmental factors [38, 34]. The interesting
question in the context of fitness landscapes is, what is the robustness of
the sequences of spins to mutations in the bond configurations (the epistatic
couplings)?

Figure 1: $G(V,E)$ for (a and b) two representative examples of random graphs
and (c and d) two representative examples of complete graphs (Sherrington-
Kirkpatrick model). Each Ising spin $s_{i}$ is placed on a vertex, and each
interaction $J_{ij}$ is placed on an edge. Each spin also experiences an
external random field $h_{i}$.

In this paper, we consider this question by investigating the spin glass phase
at $T=0$ by calculating the robustness of $\pm J$ spin glass ground state
configurations to a sign flip perturbation of a single bond. Our investigation
is related to the concept of spin glass bond chaos or disorder chaos [41],
which considers how the ground state of a spin glass changes when all the
couplings are perturbed by a small amount; notable investigations (see refs.
[41, 42, 43, 44, 45]) have been conducted on continuous $J$ spin glasses,
whereas we consider the $\pm J$ spin glasses. Another similar concept is the
recently-introduced $\sigma$-criticality [46], which considers the effects on
the ground state of only a single-bond perturbation, as we do here, but again
this is a continuous $J$ spin glass. Our use of $\pm J$ spin glasses is an
important distinction between our work and these previous studies because the
discreteness of the bond configuration domain of the $\pm J$ spin glass,
modeled as a hypercube graph, allows us to investigate the topological
features of the subgraphs formed by bond configurations mapping to the same
ground state. The bond chaos and $\sigma$-criticality investigations are
mainly interested in how the spin glass ground states change, and we consider
these aspects here as well, but in this work we are mainly interested in the
universality of the subgraph network features in the bond configuration
domain.

In addition to spin glasses on random graphs, we examine special cases, namely
the Sherrington-Kirkpatrick model and 1D Edwards-Anderson model. Various
network topological properties are computed for interaction domain subgraphs
comprised of bond configurations that all map to the same ground state spin
configuration. We find that these subgraphs obey the same logarithmic scaling
law between robustness and neutral set size as the analogous biological and
computer science GP maps described above, suggesting that this high robustness
may hold for a much wider set of physical systems.

##  II Model and Methods

###  II.1 Spin Glass Model

Consider an undirected, unweighted random graph $G(V,E)$ with vertex set $V$
and edge set $E$. We place Ising spins on each vertex, and each edge
represents a nonzero interaction between spins. A spin configuration
$s\in\\{\pm 1\\}^{\absolutevalue{V}}$ can be written as a sequence of $+1$ and
$-1$ values, so it is essentially a binary sequence of length
$\absolutevalue{V}$. A set of interactions (bond configuration) $J\in\\{\pm
1\\}^{\absolutevalue{E}}$ similarly is a sequence of $+1$ and $-1$ values of
length $\absolutevalue{E}$. The spin glass Hamiltonian

$\displaystyle\mathcal{H}_{G}(s;J)=-\sum_{\\{i,j\\}\in
E}J_{ij}s_{i}s_{j}-\sum_{i\in V}h_{i}s_{i}$ (2)

contains couplings between all spins which are connected by an edge in $G$.
The single-spin, external magnetic field interactions $h_{i}$ are
independently and identically distributed uniformly on the interval
$[-10^{-4},+10^{-4}]$, noting that the external field $h_{i}\ll
J_{ij}\in\\{-1,+1\\}$; we incorporate an external field in order to break the
possible degeneracies of the spin glass ground state. In our numerical
simulations described below, we find that this choice of distribution for
$h_{i}$ is sufficient to ensure a unique ground state for all of our
simulations.

The input-output map considered in our study is the spin glass ground state
optimization function

$\displaystyle\Gamma_{G}:\\{\pm 1\\}^{\absolutevalue{E}}\rightarrow\\{\pm
1\\}^{\absolutevalue{V}}$ (3)

defined for the graph $G$. For a bond configuration $J$, $\Gamma_{G}(J)$
outputs the ground state configuration $s$ that minimizes the Hamiltonian eq.
2. The most common task in spin glass theory is to find $s$ given a particular
bond configuration $J$. In this paper, we study an inverse problem, namely the
relationship between the set of all bond configurations $\\{J\\}$ that
generate a particular output $s$.

To efficiently represent this system, we note that the collection of all
binary sequences of length $n$ can be represented by a $n$-dimensional
undirected hypercube graph $Q_{n}(U,F)$. This is accomplished by mapping each
binary sequence to a vertex in $Q_{n}(U,F)$ and placing edges between two
vertices if the corresponding sequences have a Hamming distance of 1 between
them. The hypercube graph has vertex set $U$ with $\absolutevalue{U}=2^{n}$
and edge set $F$ with $\absolutevalue{F}=2^{n-1}n$.

The domain of $\Gamma_{G}$ accordingly has a mapping to the
$\absolutevalue{E}$-dimensional hypercube graph $Q_{\absolutevalue{E}}(U,F)$.
In general, for graphs $G$ that can produce geometrical frustration in the
spin glass, $\Gamma_{G}(J)$ follows no pattern and is difficult to calculate
[47], even more so because of the degeneracy-breaking external random field
interactions $\\{h_{i}\\}$. But, because of frustration, two sets of bond
configurations $J^{(i)}$ and $J^{(j)}$ corresponding to adjacent vertices in
$Q_{n}(U,F)$ often have $\Gamma_{G}(J^{(i)})=\Gamma_{G}(J^{(j)})$. The
vertices corresponding to all $J$ such that $\Gamma_{G}(J)=s$ for some fixed
$s$ induce a subgraph $H_{s}(U_{s},F_{s})$ of $Q_{\absolutevalue{E}}$. It
follows that $\bigcup_{s}U_{s}=U$. In the GP map literature, a neutral network
(or neutral set) for phenotype (here, the ground state) $s$ is the graph
$H_{s}(U_{s},F_{s})$.

In this paper, we numerically compute topological properties of each neutral
network of spin glass bond configurations which all map to the same ground
state spin configuration. We consider multiple random graphs for $G$ as well
as the fully-connected graph; the latter case is known as the Sherrington-
Kirkpatrick (SK) model of a spin glass [29]. In our simulations, we also
impose the condition that every node in $G$ has at least one neighbor.
Additionally, in this paper we also consider the 1-dimensional Edwards-
Anderson (EA) model [30], for which the relationship between robustness
$\absolutevalue{F_{s}}$ and vertex count $\absolutevalue{U_{s}}$ (equivalent
to robustness) becomes analytically solvable. Our simulated spin glasses are
relatively small because (1) we need to calculate the ground states exactly
for each bond configuration, (2) every possible bond configuration is
considered for each spin glass in order to accurately determine neutral
network properties. In other words, the number of times eq. 2 is computed in
order to find all ground states for a random graph $G(V,E)$ scales as
$O\left(2^{\absolutevalue{V}}\times 2^{\absolutevalue{E}}\right)$, which
forces us to use relatively small systems. Our largest model
($\absolutevalue{V}=8$, $\absolutevalue{E}=28$) involves calculation of an
exact ground state for every single bond configuration (over 268 million bond
configurations); this is a similar size to systems for which exact ground
state calculations like ours were carried out in the context of the spin glass
literature (see e.g., ref. [48]).

###  II.2 Definitions of Topological Quantities

The following parameters are computed for the neutral networks in bond
configuration space:

#### Robustness.

The neighbor count $n_{p,g}$ in eq. 1 is equivalent to the degree of a vertex
$g\in U_{s}$ within a neutral network $H_{s}(U_{s},F_{s})$. The mean degree is
related to the number of edges by

$\displaystyle\sum_{v\in U_{s}}\deg(v)=2\absolutevalue{F_{s}}.$ (4)

Here, we compute the robustness $\rho_{s}$, simply dividing the above quantity
by the size of the subgraph $\absolutevalue{U_{s}}$ and by the length of the
input sequence $\absolutevalue{E}$:

$\displaystyle\rho_{s}=\phi_{ss}\equiv\frac{2\absolutevalue{F_{s}}}{\absolutevalue{U_{s}}\absolutevalue{E}}\in[0,1].$
(5)

The notation $\phi_{ss}$ will become clear below when we also treat the
transition probability $\phi_{rs}$ of a bond perturbation leading to a
different ground state $r$. It is important to notice that the robustness is
nothing other than the edge-to-vertex ratio
$\absolutevalue{F_{s}}/\absolutevalue{U_{s}}$ normalized to the range $[0,1]$
and also is equivalently the normalized mean degree within the neutral network
$H_{s}$. In many “real-world” GP maps including RNA and protein folding,
Boolean threshold networks, and genetic algorithms, it has been observed that
$\rho_{s}\sim\log\absolutevalue{U_{s}}$ or equivalently
$\absolutevalue{F_{s}}\sim\absolutevalue{U_{s}}\log\absolutevalue{U_{s}}$. We
will test this scaling for spin glass systems.

#### Transition probability.

We define the transition probability $\phi_{rs}$ to be the probability that,
given a bond configuration which maps to ground state $s$, a single-bond
perturbation changes the ground state from $s$ to $r$. Graph theoretically, we
can think of this in the following way: consider two neutral networks of
$Q_{\absolutevalue{E}}(U,F)$ called $H_{s}(U_{s},F_{s})$ and
$H_{r}(U_{r},F_{r})$ so that $s$ and $r$ are two different spin configurations
($s\neq r$). Let $T_{rs}=T_{sr}$ be the set of edges connecting $H_{s}$ and
$H_{r}$. The transition probability of $s\rightarrow r$ due to a single-bond
perturbation is

$\displaystyle\phi_{rs}\equiv\frac{\absolutevalue{T_{rs}}}{\absolutevalue{U_{s}}\absolutevalue{E}}\in[0,1],\quad
s\neq r.$ (6)

In many real-world GP maps, it has been observed that often
$\phi_{rs}\sim\absolutevalue{U_{r}}$, or equivalently
$\absolutevalue{T_{rs}}\sim\absolutevalue{U_{r}}\absolutevalue{U_{s}}$ [8,
12]. We will also test this scaling for spin glass systems.

We now see that, for the case where $s=r$, the transition probability
$T_{ss}$, based on the definition above, simply counts the number of edges
within the induced subgraph $H_{s}(U_{s},F_{s})$, and therefore
$T_{ss}=F_{s}$. In the definition of $\phi_{ss}$ in eq. 5, an additional
factor of 2 comes from the double counting of edges; recall, the robustness is
equivalent to the (normalized) mean degree, which means that the number of
edges must be double-counted. From the definitions, it is easy to check that
the normalization condition $\sum_{r}\phi_{rs}=1$ holds, where the sum
includes $r=s$.

#### Rank-size plot.

A rank-size plot is useful for visualizing how many orders of magnitude the
neutral set sizes span. It may be used to deduce if there is a power law
(i.e., generalized Zipf’s Law) relationship between rank and neutral set size,
as is seen for some GP maps.

#### Degree distribution, clustering coefficients, assortativity, and
betweenness centrality.

We also compute the degree distribution, clustering coefficients,
assortativity, and betweenness centrality for the neutral networks studied
here. Degree distributions give us information about the modality and skew of
the degrees of the vertices in the neutral networks. Clustering coefficients,
which turn out to be trivially zero for these systems, are informative about
the number of neighboring nodes which share a neighbor. Assortativity measures
correlation between the degrees of neighboring vertices. Betweenness
centrality of a node measures the number of shortest paths in the network
passing through that node; it is often compared with, and is expected to be
positively correlated with, the degree of the node [49]. Because our systems
are small, generalizable trends regarding these quantities are more difficult
to resolve, but calculating them assists our physical intuition, so we present
them in the Supplementary Material along with their mathematical definitions.

Figure 2: Linear-linear plot of robustness $\rho_{s}$ versus normalized
subgraph size $\absolutevalue{U_{s}}/2^{\absolutevalue{E}}$ (equivalent to
NSS) for representative random graphs (a and b) and representative SK models
(c and d). The dashed curve arises from the linear least squares fit
calculated for Figure 3, with the abscissa transformed to linear scale. A
logarithmic relationship between robustness and NSS is established.  Figure 3:
Linear-log plot of robustness $\rho_{s}$ versus normalized subgraph size
$\absolutevalue{U_{s}}/2^{\absolutevalue{E}}$ (equivalent to NSS) for
representative random graphs (a and b) and representative SK models (c and d).
The dashed line is the line of best of fit in the linear-log scale. A linear
relationship between robustness and the log of NSS is established. Regression
means and 95% confidence intervals are given for the slope of the best fit
lines. The dotted line shows the random null expectation for robustness.
Figure 4: Transition probability $\phi_{rs}$ (for $r\neq s$ versus normalized
subgraph size $\absolutevalue{U_{s}}/2^{\absolutevalue{E}}$ (equivalent to
NSS) for representative random graphs (a and b) and representative SK models
(c and d). The dashed line is the line of best of fit in the linear-log scale;
the dotted line is the identity line (in which abscissa = ordinate). In panel
(d), for the SK model, the line of best fit deviates most from the identity
line (it is outside the viewing window), but the data are consistent with a
linear fit, as for the other models. Regression means and 95% confidence
intervals are given for the slope of the best fit lines.

##  III Results

We now present results for the properties described in the previous section
for the neutral networks $H_{s}(U_{s},F_{s})$ for spin glasses defined on
random graphs as well as for the SK model. Two instances of random graphs and
two SK models of differing sizes are used as the representative simulation
examples for this main text; these graphs $G$ are shown in Figure 1. To
demonstrate the consistency of our results across many instances of the random
graphs, several additional random graph instances and their network
topological properties are found in the plots in the Supplementary Material.
We also present the 1D EA model as a special case where we can calculate the
exact relationship between robustness and neutral network size.

###  III.1 Robustness

For all topologies for $G$, we find that each induced subgraph
$H_{s}(U_{s},F_{s})$ has exactly one connected component, regardless of size.
Whether this will also hold for larger systems is unclear, but this lack of
multiple component networks is different from, for example, the RNA GP map,
where biophysical constraints (mainly that $\text{GC}\rightarrow\text{UA}$ and
$\text{CG}\rightarrow\text{AU}$ are not possible by single point mutations)
lead to fragmentation of the neutral networks [5, 20]. In Figure 2, we plot
the robustness of neutral networks versus the neutral set size. Spin glasses
on random and complete graphs all show behavior consistent with the
$\rho_{s}\sim\log\absolutevalue{U_{s}}$ relationship, which is also seen for
the closely related scaling of robustness with neutral set size found for many
GP maps. For these small models, the scale of the variation is too small to
fully confirm the expected log-scaling, but the robustness is significantly
enhanced when compared to an uncorrelated random null model (see Figure 3).
This null model is calculated by taking each spin glass ground state mapping
$\Gamma_{G}(J)$ and randomizing the input-output pairings while keeping the
subgraph size the same. A single-bond perturbation $J\mapsto J^{\prime}$ will
result in a spin configuration $s=\Gamma_{G}(J^{\prime})$ being selected with
probability $\absolutevalue{U_{s}}/2^{\absolutevalue{E}}$, regardless of
$\Gamma_{G}(J)$. Thus,
$\phi_{rs}\approx\absolutevalue{U_{r}}/2^{\absolutevalue{E}}$, even for $r=s$,
in the null model.

Figure 5: Plot of normalized induced subgraph size
$\log_{10}(\absolutevalue{U_{s}}/2^{\absolutevalue{E}})$ (equivalent to the
NSS or frequency) versus the rank of the size for (a and b) random graphs and
(c and d) SK models. All the models exhibit a similar rapid decay of the
neutral set size.

###  III.2 Transition Probabilities

Transition probabilities are plotted in Figure 4 for the largest induced
subgraph of each spin glass model. As found for biologically inspired GP maps
[8, 12], $\phi_{rs}$ is typically much smaller than the $\rho_{s}$; this is
apparent from comparing the vertical axes in Figure 4 and Figure 3. As a null
model, we again use

$\displaystyle\phi_{rs}\approx\absolutevalue{U_{r}}/2^{\absolutevalue{E}},\quad
r\neq s.$ (7)

Overall, for the random graphs, as can be seen in Figure 4, this null model
curve matches the linear least squares fit, suggesting that vertices of
subgraphs $H_{r}(U_{r},F_{r})$ (for $r\neq s$) with nonzero transition
probability are approximately randomly distributed with frequency
$\approx\absolutevalue{U_{r}}/2^{\absolutevalue{E}}$ in the neighborhoods of
all vertices $v\in\absolutevalue{U_{s}}$. To interpret these findings about
robustness and transition probabilities, consider a vertex $v\in
H_{s}(U_{s},F_{s})$. In its set of nearest neighbors, this vertex is expected
to see an overrepresentation (relative to the null model) of other vertices
belonging to $H_{s}(U_{s},F_{s})$, and it sees a random assortment of other
vertices $u\in H_{r}(U_{r},F_{r})$ for various $r\neq s$ with probability
proportional to $\absolutevalue{U_{r}}/2^{\absolutevalue{E}}$. These findings
are in accordance with other GP map studies.

For the SK models, there is less agreement between the null expectation and
the observed lines of best fit. We speculate that the increased bond density
in $G$ is responsible for this effect on transition probabilities, but our
systems are currently too small to investigate this effect.

###  III.3 Size-Rank Distributions

In the GP map literature, there has been a lot of interest in phenotype bias,
the observation that the neutral set sizes can vary over many orders of
magnitude, which can determine evolutionary outcomes [11, 21] even when
natural selection is also at play. In Figure 5 we show that such phenotype
bias also exists for spin glass systems. The rank plots show a consistent
behavior independent of spin glass graph topology $G$. Recent work on GP maps
has suggested that there are two main classes of rank plots [9, 50, 14]. For
the first class, the distribution of neutral set sizes obeys a Zipf-like power
law, which arises from models in which input site ordering is strongly
constrained (including Boolean neural networks [51, 52]). The other class is a
log-normal distribution, which appears, for example, in RNA secondary
structure GP maps [50, 11]. An open question is whether or not these spin
glass systems also fall into one of these two classes. The current systems are
still too small to conclusively answer this question, but the log-log rank
versus size plots in Supplementary Material suggest a deviation from Zipf’s
law.

###  III.4 Effects of Bond Perturbation on Ground State Spin Configuration

Our single-bond perturbation numerical experiments are informative regarding
the nature of the ground state in the $T=0$ glassy phase. Our robustness
calculations have already shown that a single-bond perturbation often leads to
no change in the ground state, much more often than what the random null model
predicts. Given that the robustness is equivalent to the normalized mean
degree within a neutral network, we can understand robustness at the
individual bond level by looking at the relationship between the degree of an
individual vertex within the neutral network and various physical parameters
such as ground state energy. We can ask, how is the degree within a neutral
network related to the ground state energy? How many spin flips occur in the
ground state energy when a single bond is perturbed? Is the unperturbed ground
state a local minimum in the energy landscape generated by the new, perturbed
bond configuration?

From the degree distributions within neutral networks with
$\absolutevalue{E}=15$ (shown in the Supplementary Material), we see that the
range of degrees for the vertices in most neutral networks tends to span from
4 up to 15. In the Supplementary Material, we empirically confirm the expected
result that the degree of a vertex is positively correlated with its
betweenness centrality; that is, there tend to be a larger number of shortest
paths passing through vertices with higher degree. Intuitively, one should
then expect that vertices with higher degree have ground states which are more
energetically stable to bond perturbation. In other words, the ground state
energy is sufficiently low, and the minimum is deep such that a single-bond
perturbation tends not to change the ground state. This trend is what we see
in the numerical results presented in the Supplementary Material, where we
also discuss the different scenarios when ground state change due to bond
perturbation does occur.

Figure 6: Results for 1D Edwards-Anderson model: (a) graph representation of
the 1D Edwards-Anderson model with $\absolutevalue{V}=11$, (b) neutral set
size versus rank plot on log-log scale, (c) degree distribution of neutral set
vertices, and (d) robustness $\rho_{s}$ versus log of the normalized neutral
set size $\absolutevalue{U_{s}}/2^{\absolutevalue{E}}$. The dashed line is the
analytical result from eq. 8.

###  III.5 Analytically Tractable Special Case: 1D Edwards-Anderson Model

The Edwards-Anderson (EA) model is a special case which deserves individual
treatment. The EA model, the original theory of spin glasses [30], is simply a
spin glass on a lattice with nearest neighbor interactions only.

For the 1D Edwards-Anderson (EA) model with periodic boundary conditions, the
topology of which is shown in panel (a) of Figure 6, the behavior of
$\rho_{s}$ is analytically tractable. Let us call a bond configuration with an
even number of anti-ferromagnetic interactions ($J_{ij}=-1$) an even bond
configuration and a bond configuration with an odd number of anti-
ferromagnetic interactions an odd bond configuration.

Every even bond configuration has an exactly determinable ground state with no
spin frustration, which can be found using the following algorithm:

  1. 1.

Choose an arbitrary site $i$, and set the $i$-th spin to $s_{i}=+1$ or
$s_{i}=-1$ arbitrarily.

  2. 2.

Determine $s_{i+1}$ by setting $s_{i+1}=J_{i,i+1}s_{i}$, so $s_{i+1}=s_{i}$ if
$J_{i,i+1}=+1$ or $s_{i+1}=-s_{i}$ if $J_{i,i+1}=-1$.

  3. 3.

Continue to $s_{i+2}$ by similarly setting $s_{i+2}=J_{i+1,i+2}s_{i+1}$, and
continue around the entire spin chain until we have set $s_{i-1}$. At this
point one can check and confirm that $-J_{i-1,i}s_{i-1}s_{i}=-1$, so all bonds
are satisfied, and each pairwise interaction decreases the energy by the same
amount ($=-1$).

  4. 4.

If $\sum_{i=1}^{\absolutevalue{V}}h_{i}s_{i}<0$, then all spins need to be
flipped $s_{i}\rightarrow-s_{i}$. If
$\sum_{i=1}^{\absolutevalue{V}}h_{i}s_{i}>0$, then all the spins are kept as
is. The ground state energy is exactly
$E_{gs}=-\absolutevalue{V}-\sum_{i=1}^{\absolutevalue{V}}h_{i}s_{i}$, and we
have now ensured $\sum_{i=1}^{\absolutevalue{V}}h_{i}s_{i}>0$.

Having an even number of anti-ferromagnetic interactions will ensure that the
last spin $s_{i-1}$ does not experience frustration, and the bond between the
last spin $s_{i-1}$ and the first spin $s_{i}$ contributes a negative amount
to the energy: $-J_{i-1,i}s_{i-1}s_{i}=-1$. On the other hand, odd bond
configurations will always have one bond at which there is frustration, and
there will be a positive contribution to the energy
$-J_{i-1,i}s_{i-1}s_{i}=+1$. The ground state for odd bond configurations will
be identical to the ground state of one of the neighboring even bond
configurations, but the energy will, of course, be higher. Exactly which
neighboring bond configuration that is depends on the random fields $h_{i}$.

Because an odd bond configuration has the same ground state as one of its
neighboring even bond configurations, the neutral network $H_{s}(U_{s},F_{s})$
for every ground state must be a star graph $S_{\absolutevalue{U_{s}}-1}$ that
has $\absolutevalue{U_{s}}$ nodes, with one central node corresponding to an
even bond configuration, and $\absolutevalue{U_{s}}-1$ peripheral nodes
corresponding to odd bond configurations. The number of nodes of each star
graph neutral network depends on the random fields $h_{i}$. Fig. 6(c) shows
that the degree distribution matches our theory that all neutral networks are
star graphs; for each neutral network, there is one node with
$\absolutevalue{U_{s}}-1$ neighbors, and $\absolutevalue{U_{s}}-1$ nodes with
one neighbor.

To calculate robustness, we note that the star graph
$S_{\absolutevalue{U_{s}}-1}$ has $\absolutevalue{U_{s}}$ nodes and
$\absolutevalue{F_{s}}=\absolutevalue{U_{s}}-1$ edges. Plugging
$\absolutevalue{F_{s}}$ into eq. 5 immediately gives us that the robustness is

$\displaystyle\rho_{s}=\frac{1}{\absolutevalue{U_{s}}\absolutevalue{E}}\sum_{v\in
U_{s}}\deg(v)=\frac{2}{\absolutevalue{E}}\left(1-\frac{1}{\absolutevalue{U_{s}}}\right).$
(8)

The points in Figure 6(d) fall exactly along this curve. The maximum number of
nodes in the star graph is
$\absolutevalue{U_{s}}=\absolutevalue{E}+1=\absolutevalue{V}+1$ because each
even bond configuration has at most $\absolutevalue{E}=\absolutevalue{V}$
neighboring odd bond configurations. A size-rank plot is also shown for the 1D
EA model in panel (b).

The 1D EA model serves as an example in which linear-log scaling is not seen
for $\rho_{s}$. Nevertheless, the relationship between robustness and subgraph
size found in the 1D EA model exhibits a $\rho_{s}$ that is significantly
higher than would be expected from a random mapping $\Gamma_{G}$ of inputs to
outputs, so this 1D model still exhibits enhanced robustness. For the 2D EA
model we could not find an analytically tractable $\rho_{s}$.

##  IV Discussion

Our main result is that, by studying the mapping from spin glass bond
configurations to the $T=0$ ground state, we observe several properties that
are also observed for other input-output maps, such as the GP maps found in
biology and computer science. These include (1) redundancy, in that many
genotypes in the context of GP maps, or bond configurations in this case, map
the the same output (a phenotype for the GP maps, or a ground state for the
spin glasses), (2) phenotype frequency bias, in that the number of bond
configurations (genotypes) mapping to the ground state (phenotype) vary
significantly, and (3) enhanced mutational robustness, in that the probability
of a single-bond perturbation changing the ground state is larger than
expected from a random null model, and moreover that the robustness of a
ground state to single bond perturbations often scales in the same way with
the logarithm of the NSS (the number of bond configurations mapping to a
ground state) as found for biological GP maps [9, 10, 12, 16, 17]. Another
interesting result, also seen for GP maps [8, 12] is that, in contrast to the
robustness, the transition probabilities, defined as the likelihood of a flip
of the spin yielding a different ground state, do scale proportionally to the
NSS, as one would expect from a random model. The similarity to the GP map
behavior suggests that there may be a more universal argument (based, for
example, on algorithmic information theory [53, 54]) for these scaling
properties.

Our spin glass models are relatively small because finding the ground state of
a spin glass is typically computationally expensive and scales badly with
system size. Depending on graph topology, finding the ground state of a spin
glass can be NP hard [47]. Knowledge that the robustness of a ground state is
large may potentially offer improvements to ground state-finding algorithms by
providing a measure of stability of certain ground states as a function of
parameter space. It would also be interesting to check some of our results on
significantly larger graphs. We find, for example, that all our subgraphs that
map to the same ground state form only a single component; that is, they are
connected by single-bond perturbations. Will this percolation property hold
for larger systems, or will these subgraphs start to fragment?

Mapping epistatic interactions onto spin glass Hamiltonians to derive
sequence-to-fitness maps has been especially important for the study of viral
evolution [35, 36, 37, 38, 39]. These models typically have continuous $J$, so
it will be interesting to see if the kind of scaling properties of robustness
that we find here for the $\pm J$ spin glasses carry over to these more
complex systems. If, as we expect, a concept akin to high robustness persists,
then this may have implications for the stability of fitness peaks to
environmental changes, and may affect how easy it is to find continuously
increasing paths to fitness maxima [24].

Another potentially interesting future direction of research is to explore
these results about robustness in the complementary setting of the sensitivity
of Boolean functions [26, 25, 27, 28]. It would be interesting to see whether
similar high robustness/low sensitivity results can be found in this related
context.

##  V Acknowledgements

The authors would like to thank Wilber Lim, Mehrana Nejad, Sloan Nietert, and
Shuofeng Zhang for helpful discussions. VM has been supported by a Marshall
Scholarship and award numbers T32GM007753 and T32GM144273 from the National
Institute of General Medical Sciences. The content is solely the
responsibility of the authors and does not necessarily represent the official
views of the National Institute of General Medical Sciences, the National
Institutes of Health, or the Marshall Aid Commemoration Commission.

## References

  * Kauffman [1969] S. Kauffman, Homeostasis and differentiation in random genetic control networks, Nature 224, 177 (1969).
  * Wagner [2005] A. Wagner, Distributed robustness versus redundancy as causes of mutational robustness, BioEssays 27, 176 (2005).
  * Wagner [2007] A. Wagner,  _Robustness and evolvability in living systems_ , 3rd ed., Princeton studies in complexity (Princeton Univ. Press, Princeton, NJ, 2007).
  * Wagner [2008] A. Wagner, Robustness and evolvability: a paradox resolved, Proceedings of the Royal Society B: Biological Sciences 275, 91 (2008).
  * Aguirre  _et al._ [2011] J. Aguirre, J. M. Buldú, M. Stich, and S. C. Manrubia, Topological structure of the space of phenotypes: the case of RNA neutral networks, PLoS ONE 6, e26324 (2011).
  * Payne and Wagner [2013] J. L. Payne and A. Wagner, Constraint and contingency in multifunctional gene regulatory circuits, PLoS Computational Biology 9, e1003071 (2013).
  * Payne  _et al._ [2014] J. L. Payne, J. H. Moore, and A. Wagner, Robustness, evolvability, and the logic of genetic regulation, Artificial Life 20, 111 (2014).
  * Schaper and Louis [2014] S. Schaper and A. A. Louis, The arrival of the frequent: How bias in genotype-phenotype maps can steer populations to local optima, PLoS ONE 9, e86635 (2014).
  * Greenbury  _et al._ [2014] S. F. Greenbury, I. G. Johnston, A. A. Louis, and S. E. Ahnert, A tractable genotype–phenotype map modelling the self-assembly of protein quaternary structure, Journal of The Royal Society Interface 11, 20140249 (2014).
  * Greenbury and Ahnert [2015] S. F. Greenbury and S. E. Ahnert, The organization of biological sequences into constrained and unconstrained parts determines fundamental properties of genotype–phenotype maps, Journal of The Royal Society Interface 12, 20150724 (2015).
  * Dingle  _et al._ [2015] K. Dingle, S. Schaper, and A. A. Louis, The structure of the genotype–phenotype map strongly constrains the evolution of non-coding RNA, Interface Focus 5, 20150053 (2015).
  * Greenbury  _et al._ [2016] S. F. Greenbury, S. Schaper, S. E. Ahnert, and A. A. Louis, Genetic correlations greatly increase mutational robustness and can both reduce and enhance evolvability, PLOS Computational Biology 12, e1004773 (2016).
  * Ahnert [2017] S. E. Ahnert, Structural properties of genotype–phenotype maps, Journal of The Royal Society Interface 14, 20170275 (2017).
  * Weiß and Ahnert [2018] M. Weiß and S. E. Ahnert, Phenotypes can be robust and evolvable if mutations have non-local effects on sequence constraints, Journal of The Royal Society Interface 15, 20170618 (2018).
  * Nichol  _et al._ [2019] D. Nichol, M. Robertson-Tessi, A. R. A. Anderson, and P. Jeavons, Model genotype–phenotype mappings and the algorithmic structure of evolution, Journal of The Royal Society Interface 16, 20190332 (2019).
  * Camargo and Louis [2020] C. Q. Camargo and A. A. Louis, Boolean threshold networks as models of genotype-phenotype maps, Complex Networks XI , 143 (2020).
  * Hu  _et al._ [2020] T. Hu, M. Tomassini, and W. Banzhaf, A network perspective on genotype–phenotype mapping in genetic programming, Genetic Programming and Evolvable Machines 15 (2020).
  * Manrubia  _et al._ [2021] S. Manrubia, J. A. Cuesta, J. Aguirre, S. E. Ahnert, L. Altenberg, A. V. Cano, P. Catalán, R. Diaz-Uriarte, S. F. Elena, J. A. García-Martín, P. Hogeweg, B. S. Khatri, J. Krug, A. A. Louis, N. S. Martin, J. L. Payne, M. J. Tarnowski, and M. Weiß, From genotypes to organisms: State-of-the-art and perspectives of a cornerstone in evolutionary dynamics, Physics of Life Reviews 38, 55 (2021).
  * Payne and Wagner [2019] J. L. Payne and A. Wagner, The causes of evolvability and their evolution, Nature Reviews Genetics 20, 24 (2019).
  * Schaper  _et al._ [2012] S. Schaper, I. G. Johnston, and A. A. Louis, Epistasis can lead to fragmented neutral spaces and contingency in evolution, Proceedings of the Royal Society B: Biological Sciences 279, 1777 (2012).
  * Dingle  _et al._ [2022] K. Dingle, F. Ghaddar, P. Šulc, and A. A. Louis, Phenotype bias determines how natural rna structures occupy the morphospace of all possible shapes, Molecular Biology and Evolution 39 (2022).
  * Wright and Laue [2022] A. H. Wright and C. L. Laue, Evolving Complexity is Hard (2022), arXiv:2209.13013 [cs].
  * Johnston  _et al._ [2022] I. G. Johnston, K. Dingle, S. F. Greenbury, C. Q. Camargo, J. P. K. Doye, S. E. Ahnert, and A. A. Louis, Symmetry and simplicity spontaneously emerge from the algorithmic nature of evolution, Proceedings of the National Academy of Sciences 119, e2113883119 (2022).
  * Greenbury  _et al._ [2022] S. F. Greenbury, A. A. Louis, and S. E. Ahnert, The structure of genotype-phenotype maps makes fitness landscapes navigable, Nature Ecology & Evolution (2022).
  * Huang [2019] Huang, Induced subgraphs of hypercubes and a proof of the Sensitivity Conjecture, Annals of Mathematics 190, 949 (2019).
  * Nisan [1994] N. Nisan, On the degree of boolean functions as real polynomials, computational complexity 4, 301 (1994).
  * Bernasconi [1996] A. Bernasconi, Sensitivity vs. block sensitivity (an average-case study), Information Processing Letters 59, 151 (1996).
  * Chakraborty [2011] S. Chakraborty, On the sensitivity of cyclically-invariant Boolean functions, Discrete Mathematics & Theoretical Computer Science 13, 51 (2011).
  * Sherrington and Kirkpatrick [1975] D. Sherrington and S. Kirkpatrick, Solvable Model of a Spin-Glass, Physical Review Letters 35, 1792 (1975).
  * Edwards and Anderson [1975] S. F. Edwards and P. W. Anderson, Theory of spin glasses, Journal of Physics F: Metal Physics 5, 965 (1975).
  * Mezard  _et al._ [1987] M. Mezard, G. Parisi, and M. A. Virasoro,  _Spin glass theory and beyond_ , World Scientific lecture notes in physics No. v. 9 (World Scientific, Singapore ; New Jersey, 1987).
  * Mezard and Montanari [2009] M. Mezard and A. Montanari,  _Information, physics, and computation_ , Oxford graduate texts (Oxford University Press, Oxford ; New York, 2009).
  * Nishimori [2001] H. Nishimori,  _Statistical physics of spin glasses and information processing: an introduction_ , International series of monographs on physics No. 111 (Oxford University Press, Oxford ; New York, 2001).
  * Guo  _et al._ [2019] Y. Guo, M. Vucelja, and A. Amir, Stochastic tunneling across fitness valleys can give rise to a logarithmic long-term fitness trajectory, Science Advances 5, eaav3842 (2019).
  * Louie  _et al._ [2018] R. H. Y. Louie, K. J. Kaczorowski, J. P. Barton, A. K. Chakraborty, and M. R. McKay, Fitness landscape of the human immunodeficiency virus envelope protein that is targeted by antibodies, Proceedings of the National Academy of Sciences 115, E564 (2018).
  * Butler  _et al._ [2016] T. C. Butler, J. P. Barton, M. Kardar, and A. K. Chakraborty, Identification of drug resistance mutations in HIV from constraints on natural evolution, Physical Review E 93, 022412 (2016).
  * Barton  _et al._ [2016] J. P. Barton, N. Goonetilleke, T. C. Butler, B. D. Walker, A. J. McMichael, and A. K. Chakraborty, Relative rate and location of intra-host HIV evolution to evade cellular immunity are predictable, Nature Communications 7, 11660 (2016).
  * Shekhar  _et al._ [2013] K. Shekhar, C. F. Ruberman, A. L. Ferguson, J. P. Barton, M. Kardar, and A. K. Chakraborty, Spin models inferred from patient-derived viral sequence data faithfully describe HIV fitness landscapes, Physical Review E 88, 062705 (2013).
  * Hopf  _et al._ [2017] T. A. Hopf, J. B. Ingraham, F. J. Poelwijk, C. P. I. Schärfe, M. Springer, C. Sander, and D. S. Marks, Mutation effects predicted from sequence co-variation, Nature Biotechnology 35, 128 (2017).
  * Cocco  _et al._ [2018] S. Cocco, C. Feinauer, M. Figliuzzi, R. Monasson, and M. Weigt, Inverse statistical physics of protein sequences: a key issues review, Reports on Progress in Physics 81, 032601 (2018).
  * Bray and Moore [1987] A. J. Bray and M. A. Moore, Chaotic nature of the spin-glass phase, Physical Review Letters 58, 57 (1987).
  * Aspelmeier [2008] T. Aspelmeier, Free-energy fluctuations and chaos in the sherrington-kirkpatrick model, Physical Review Letters 100, 117205 (2008).
  * Ney-Nifle and Young [1997] M. Ney-Nifle and A. P. Young, Chaos in a two-dimensional Ising spin glass, Journal of Physics A: Mathematical and General 30, 5311 (1997).
  * Sasaki  _et al._ [2005] M. Sasaki, K. Hukushima, H. Yoshino, and H. Takayama, Temperature chaos and bond chaos in edwards-anderson ising spin glasses: domain-wall free-energy measurements, Physical Review Letters 95, 267203 (2005).
  * Wang  _et al._ [2016] W. Wang, J. Machta, and H. G. Katzgraber, Bond chaos in spin glasses revealed through thermal boundary conditions, Physical Review B 93, 224414 (2016).
  * Newman and Stein [2022] C. M. Newman and D. L. Stein, Ground-state stability and the nature of the spin glass phase, Physical Review E 105, 044132 (2022).
  * Barahona [1982] F. Barahona, On the computational complexity of Ising spin glass models, Journal of Physics A: Mathematical and General 15, 3241 (1982).
  * Boettcher and Kott [2005] S. Boettcher and T. M. Kott, Exact enumeration of ground states in the Sherrington-Kirkpatrick spin glass, Physical Review B 72, 212408 (2005).
  * Oldham  _et al._ [2019] S. Oldham, B. Fulcher, L. Parkes, A. Arnatkevic̆iūtė, C. Suo, and A. Fornito, Consistency and differences between centrality measures across distinct classes of networks, PLOS ONE 14, e0220061 (2019).
  * Manrubia and Cuesta [2017] S. Manrubia and J. A. Cuesta, Distribution of genotype network sizes in sequence-to-structure genotype–phenotype maps, Journal of The Royal Society Interface 14, 20160976 (2017).
  * Van den Broeck and Kawai [1990] C. Van den Broeck and R. Kawai, Learning in feedforward Boolean networks, Physical Review A 42, 6210 (1990).
  * Valle-Pérez  _et al._ [2019] G. Valle-Pérez, C. Q. Camargo, and A. A. Louis, Deep learning generalizes because the parameter-function map is biased towards simple functions, arXiv:1805.08522 [cs, stat] (2019).
  * Dingle  _et al._ [2018] K. Dingle, C. Q. Camargo, and A. A. Louis, Input–output maps are strongly biased towards simple outputs, Nature Communications 9, 761 (2018).
  * Dingle  _et al._ [2020] K. Dingle, G. V. Pérez, and A. A. Louis, Generic predictions of output probability based on complexities of inputs and outputs, Scientific Reports 10, 4415 (2020).