The Unreasonable Redundancy of Nature’s Protein Folds

The graph-theoretic split

We need a way to split proteins based on how the residues are connected to each other.
For that, a protein is naturally a graph: each residue is a node, and edges connect
residues that are close in space. We use C-alpha atoms from the confident part of the
prediction, connect each residue to its spatial nearest neighbors, and give close
neighbors stronger weights than distant ones. In the current version, each residue
sees its 15 nearest spatial neighbors.

This turns the fragmentation problem into a connectivity problem. A compact
domain becomes a dense local graph. A high-confidence linker becomes a narrow
bridge between two dense regions. Graph theory gives us tools for asking
whether that bridge is really part of one unit, or whether it is holding two
independent pieces together.

Spectral bisection asks for the weakest global connection in this graph (why the quantity
we compute finds this connection is a little bit black magic to me, ask the graph theorists).
We found that the spectral bisection points of a protein correlate very well with the points
we’d cut at if we were manually identifying different protein regions. A standalone version
of the splitter is included in Supplementary info.

Clustering the fragments

Once we split proteins into their “interacting units”, the unit of clustering is no longer
one predicted protein. It is one compact fragment. We can cluster those fragments by
structural similarity and ask how much independent fold signal the distillation sets
actually contain.

We cluster MGnify at roughly 30% sequence identity with MMseqs2, which gives about 40
million sequence clusters. From there, we discard sequence singletons, then use the
OpenFold3-predicted structures
released through the OpenFold datasets portal for the remaining MGnify sequences.



This is the handoff point from sequence-space cleaning to structure-space cleaning:
sequence singletons are removed first, then the OpenFold3 predictions are fragmented
and filtered.

We fragment those predicted structures and filter the fragments with quality metrics
meant to keep examples amenable to training a generative model (we will write more
about this in a later post). The structural clustering below starts from the resulting
set: about 2 million MGnify fragments.

We use Foldseek for the
first pass of clustering. Foldseek’s fast mode uses both
structural and sequence-derived signals, which is what makes
it practical, but also means it can split fragments that are
structurally very similar and sequence-divergent.

Common pitfall: Foldseek singletons are not always singletons

A Foldseek singleton is not necessarily a new fold. It only means no
other fragment crossed the thresholds in that particular Foldseek run.
To check this failure mode, we took 1,000 fragments that Foldseek had
labeled as singletons and compared them against each other with
pairwise TM-score. At TM ≥ 0.8, 373 of those fragments fell into 69
connected components. The largest hidden cluster had 35 members.
These were supposed to be singletons.



Foldseek’s fast pass mixes 3Di and sequence-derived signals. The TM-align audit is slower,
but it asks the question we actually care about here: whether the backbones superpose.

The practical lesson is to be skeptical of the singletons. If we treat
every Foldseek singleton as an independent structural mode, we
overestimate novelty and give the sampler a distorted picture of fold
space. TM-score is much slower, but it is the right ground-truth audit
pass when the question is whether two fragments really share a fold.

So we added a second pass over cluster representatives. Instead of
comparing every fragment to every other fragment, we compared
representatives against representatives with a more structure-centered
alignment, then confirmed candidate merges with our TM-align
implementation. Merge criterion: min(tm_norm_a, tm_norm_b) ≥ 0.7.
Both directions had to independently confirm fold-level similarity.
If two representatives passed that test, we merged their clusters with
union-find.

These are the results that surprised us. After sequence clustering,
fragmentation, quality filtering, and dropping singleton clusters from
this summary, MGnify is not two million independent structural examples.
The repeated part of the dataset is closer to twenty-five thousand
structural neighborhoods, with most of the mass concentrated in a
small head of the distribution. The top 1,000 MGnify clusters are only
about 4.0% of the multi-member cluster list, but they contain 71.5% of
the fragments in those clusters.

That changes what “sampling from MGnify” means. Sampling fragments
uniformly mostly revisits the same common neighborhoods. Sampling
clusters uniformly goes too far the other way and gives the smallest
repeated neighborhoods too much influence. The sampler needs to live somewhere between those two
extremes.

Sampling from a redundant world

If we are training a generative model on MGnify, how should we sample from these clusters?
The standard recipe is: pick a cluster uniformly, then pick a member uniformly from inside
it. For data this skewed, that overshoots in the opposite direction.
With uniform cluster sampling, the top 1,000 MGnify clusters
— which hold 71.5% of multi-member fragments — would be sampled only about 4%
of the time.
We instead use a balancing exponent γ (implemented as cluster_size_exponent),
where the aggregate sampling probability of a cluster scales like N^γ with N
the cluster size. γ = 1 recovers the dataset distribution; γ = 0 weights every
cluster equally; values in between trade off natural abundance against fold diversity.

The right γ depends on what the downstream model is supposed to learn. A generative
model trying to cover fold space wants lower γ, since it benefits from seeing
rare folds more than once per epoch. A folding model trying to match natural
conditional distributions wants higher γ. We don’t have a settled answer, but
γ around 0.5 has been a reasonable starting point in our setups — it preserves the
head’s dominance while flattening the long tail. The effective-cluster count in the figure
is a useful sanity check: it’s the number of clusters you would need if they were all
equally weighted to give the same diversity as your reweighted distribution.



At γ = 0.5, aggregate cluster weight grows like the square root of cluster size:
a cluster with 100 times more fragments gets 10 times the sampling mass.

Conclusion

What surprised us is how redundant natural fold space looks once you pick the right
unit of clustering. After cleaning up predicted structures, cutting away obvious
noise, splitting multi-domain chains, auditing Foldseek singletons, and clustering
the resulting fragments, most of the mass sits in a small number of structural
neighborhoods. Natural proteins do not appear to be exploring backbone space
uniformly. They seem to reuse a relatively small set of fold solutions over
and over.

Natural enzymes often evolve by modifying existing proteins: duplication,
divergence, loop changes, active-site mutations, cofactors, metals, and changes
in the local environment around a substrate. What surprised us was not that
nature reuses folds, but how strongly that reuse shows up once we process
predicted structures into training units and cluster them at scale.

For enzyme design, this leaves two different possibilities. One is the nature-like route:
choose a familiar scaffold and learn how to engineer the active-site neighborhood with much
higher precision. In that view, the loops, first-shell residues, ligand pose, cofactors, and
interaction geometry matter more than making the backbone globally novel. If that is the
right regime, then simply adding more natural sequence-derived structures may not help much
by itself; it may mostly give us more examples of the same scaffold families.

The other possibility is more speculative. Evolution is historically constrained; new
enzymes do not appear from nowhere, and natural fold space may be shaped by what was easy
to reach through duplication and divergence. If design models become good enough, there may
be useful backbone space that nature never explored. But that raises a harder question: can
models trained mostly on natural folds learn to explore outside the natural fold manifold,
or do they inherit the same redundancy we are measuring here?

We will find out in the lab as we try to design enzymes, seeing which designs actually
express, fold, and catalyze. More on this later.

Supplementary info

Where the MGnify distillation advantage actually shows up

One reason I am suspicious of treating MGnify as just more sequence data: the
performance advantage does not appear uniformly across every interface benchmark.
It shows up strongly in antibody-antigen prediction, while most of the broader
cofolding benchmarks remain closer together. This is not a clean causal experiment
— architecture, compute, and training details all move at once — but it
fits the intuition that protein-protein interaction accuracy improves when the model
has seen much more of natural protein space.

Spectral linker split

This is a standalone version of the graph split used above. It assumes
you have already filtered out low-confidence residues, and that
ca_coords contains only the C-alpha coordinates you still
want to consider. The function returns the C-alpha indices that sit
closest to Fiedler sign-change boundaries; in our pipeline those
residues become linker/cut residues before fragment IDs are assigned.

import numpy as np
from scipy.sparse import csr_matrix
from scipy.spatial import KDTree


def spectral_linker_indices(
    ca_coords: np.ndarray,
    *,
    k: int = 15,
    sigma: float = 8.0,
    connectivity_threshold: float = 0.05,
    boundary_size: int = 4,
    min_partition_size: int = 50,
) -> np.ndarray:
    """
    Find linker-like C-alpha positions with recursive spectral bisection.

    Parameters
    ----------
    ca_coords:
        Array with shape (N, 3), containing only the C-alpha coordinates that
        survived whatever filtering you want to apply first, usually pLDDT.
    k:
        Number of spatial nearest neighbors used to build the residue graph.
    sigma:
        Gaussian bandwidth for edge weights: exp(-d^2 / (2 sigma^2)).
    connectivity_threshold:
        Stop splitting when the Fiedler value is above this threshold.
    boundary_size:
        Number of residues closest to each sign-change boundary to mark.
    min_partition_size:
        Do not try to split partitions smaller than this.

    Returns
    -------
    np.ndarray
        Sorted indices into ca_coords. These are the residues to treat as
        linkers/cuts before assigning final fragments.
    """
    coords = np.asarray(ca_coords, dtype=float)
    if coords.ndim != 2 or coords.shape[1] != 3:
        raise ValueError("ca_coords must have shape (N, 3)")

    n = len(coords)
    if n < min_partition_size:
        return np.array([])

    k_eff = min(k, n - 1)
    if k_eff < 1:
        return np.array([])

    tree = KDTree(coords)
    distances, neighbors = tree.query(coords, k=k_eff + 1)

    # Skip the self-neighbor in column 0.
    row = np.repeat(np.arange(n), k_eff)
    col = neighbors[:, 1:].ravel()
    weights = np.exp(-(distances[:, 1:] ** 2).ravel() / (2 * sigma**2))

    W = csr_matrix((weights, (row, col)), shape=(n, n))
    W = W.maximum(W.T)

    linker_indices: set[int] = set()

    def split(node_indices: np.ndarray) -> None:
        if len(node_indices) < min_partition_size:
            return

        idx = np.sort(node_indices)
        sub_W = W[idx][:, idx].toarray()

        degree = sub_W.sum(axis=1)
        if np.count_nonzero(degree) < 2:
            return

        d_inv_sqrt = np.zeros_like(degree)
        nonzero = degree > 0
        d_inv_sqrt[nonzero] = 1.0 / np.sqrt(degree[nonzero])

        # Normalized graph Laplacian: L_sym = I - D^-1/2 W D^-1/2.
        L_sym = np.eye(len(idx)) - d_inv_sqrt[:, None] * sub_W * d_inv_sqrt[None, :]

        evals, evecs = np.linalg.eigh(L_sym)
        if len(evals) < 2:
            return

        fiedler_value = evals[1]
        fiedler_vector = evecs[:, 1]
        if fiedler_value > connectivity_threshold:
            return

        boundary_order = np.argsort(np.abs(fiedler_vector))[:boundary_size]
        linker_indices.update(idx[boundary_order].tolist())

        keep = np.ones(len(idx), dtype=bool)
        keep[boundary_order] = False

        partition_a = idx[keep & (fiedler_vector < 0)]
        partition_b = idx[keep & (fiedler_vector >= 0)]

        split(partition_a)
        split(partition_b)

    split(np.arange(n))
    return np.array(sorted(linker_indices))

Foldseek command

For the fragment clustering pass, this is the Foldseek command we
used before the representative-level structural audit.

# Foldseek command for reproducibility
$FOLDSEEK cluster "$DB_DIR/fragDB" "$DB_DIR/fragCluDB" "$TMP_DIR" \
  --threads "$THREADS" \
  -c 0.8 --cov-mode 0 \
  --tmscore-threshold 0.7 \
  --tmscore-threshold-mode 1 \
  --lddt-threshold 0.6 \
  --max-seqs 2000 \
  -e 0.001 \
  --alignment-type 2 \
  --cluster-reassign 1 \
  -v 2

Arda Goreci, "The Unreasonable Redundancy of Nature's Protein Folds",
Ligo Biosciences Blog, May 20, 2026.

@article{goreci2026naturesproteinfolds,
  author = {Arda Goreci},
  title = {The Unreasonable Redundancy of Nature's Protein Folds},
  journal = {Ligo Biosciences Blog},
  year = {2026},
  month = may,
  day = {20},
  publisher = {Ligo Biosciences}
}