Filtering de novo binder designs before you order them

A generative pipeline will hand you thousands of designs. A contract research organisation will make and measure a few dozen. Everything between those two numbers is this skill: which designs go into the order, in what rank, and on what evidence.

The filters here are enrichment, not proof. The published numbers are worth stating plainly before anything else:

campaign	what was ordered	hit rate
Rosetta-era yeast-display screens, deep-learning filtered (Bennett 2023)	libraries of ~10⁴ designs per target	< 1% — and only ~2.3% of designs pass the filter at all
Adaptyv Bio EGFR competition round 2, 2024	402 designs, already selected from open submissions	14.0% (53 of 378 with usable labels)
Autonomous multi-target campaign, 2026	1,320 designs, 30 per target, heavily filtered and ensembled	26.8% overall, 49% for the top-ranked design of each target

A filtered design is a better bet, not a binder. Nothing below changes that, and a skill that implies otherwise is selling you something.

What this skill assumes you already have

This is an analysis skill: it computes over predictions you have, and does not run a predictor. For each design you need

input	what it is	where it comes from
design model	the complex your generator proposed — binder plus target	your backbone/sequence design pipeline
re-predicted complex	the same sequence pair folded again, independently	a co-folding model
PAE matrix	predicted aligned error, per residue pair, for the re-prediction	the predictor's confidence output
per-residue pLDDT	the re-prediction's local confidence	the same file
binder sequence	one-letter, binder chain only	anywhere

The registry has skills for producing the predictions: esm runs ESMFold2 and the ESM family, boltz2-nim runs Boltz2 for complexes with a confidence output, and alphafold covers reading pLDDT and PAE out of AlphaFold DB records. This skill starts after that.

Nothing here needs a GPU, an account, or a licence key. Everything runs on CPU with permissively licensed packages. Two filters that the literature leans on — Rosetta interface ΔG/ΔSASA and the Lawrence–Colman shape complementarity Sc — are not implemented here, because both ship through Rosetta or PyRosetta, which are not OSI-licensed and need a paid licence for commercial use. See What needs a licence near the end for what to substitute and what you lose.

Set up

python3 -m venv .venv
./.venv/bin/pip install --quiet "biopython>=1.85" "DockQ>=2.1"
./.venv/bin/python -c "import numpy, Bio, DockQ; print(numpy.__version__, Bio.__version__)"

Run every Python block below with ./.venv/bin/python. numpy is BSD-3-Clause, Biopython ships under the permissive Biopython Licence Agreement, and DockQ is MIT.

Do not pin numpy here. DockQ 2.1.3 requires numpy below 2, so asking for numpy>=2.0 alongside it fails resolution outright rather than warning; the line above resolves to numpy 1.26.4 and Biopython 1.88.

Get the inputs used below

Four public sources, no account on any of them. About 49 MB to download, and roughly 330 MB on disk once the structure archive is unpacked.

# 1. A labelled benchmark: 402 designs against EGFR, expressed and measured by one lab.
#    Adaptyv Bio EGFR design competition round 2. Data ODbL, code Apache-2.0.
curl -sSL -o egfr_round2.csv \
  https://raw.githubusercontent.com/adaptyvbio/egfr_competition_2/main/results/result_summary.csv
curl -sSL -o egfr_structures.zip \
  https://api.adaptyvbio.com/storage/v1/object/public/egfr_design_competition_2/structure_predictions.zip
unzip -q -o egfr_structures.zip -x "__MACOSX/*"
rm -rf egfr_structures && mv structure_predictions egfr_structures

# 2. A real two-chain prediction with its PAE matrix, plus the reference scores for it.
#    From the IPSAE repository (MIT).
base=https://raw.githubusercontent.com/DunbrackLab/IPSAE/main/Example
curl -sSL -O $base/fold_aurka_0_tpx2_0_full_data_0.json
curl -sSL -O $base/fold_aurka_0_tpx2_0_model_0.cif
curl -sSL -O $base/fold_aurka_0_tpx2_0_model_0_10_10.txt
curl -sSL -O $base/RAF1_KSR1_MEK1_9f755_scores_alphafold2_multimer_v3_model_1_seed_000.json.gz
curl -sSL -O $base/RAF1_KSR1_MEK1_9f755_unrelaxed_alphafold2_multimer_v3_model_1_seed_000.pdb

# 3. Experimental structures, from the PDB. 1OL5 is the crystal structure of the complex
#    predicted in (2); 1BRS is barnase-barstar; 1UBQ and 1BJ1 are the awkward cases.
for id in 1OL5 1BRS 1UBQ 1BJ1; do
  curl -sSL -o $id.pdb https://files.rcsb.org/download/$id.pdb
done
ls -1 *.csv *.json *.cif *.pdb | head -20

The Adaptyv round 2 release is the labelled ground truth used throughout: 402 designs submitted by 100+ entrants, all expressed by cell-free synthesis and measured on the same SPR instrument, with a binder/non-binder call, a K_D where one could be fitted, and the ColabFold AlphaFold2-multimer scores the organisers used for their own ranking. Confirmed reachable 2026-08-27.

Before any filter — check the field means what its name says

Most of the numbers below are read out of a file rather than computed, and the reading is where the mistakes are. Run this once against your own predictor's output before you write a single threshold.

import collections, statistics

def bfactors_pdb(path):
    """(chain, resnum) -> list of B-factor column values, one per atom."""
    g = collections.defaultdict(list)
    for line in open(path):
        if line.startswith("ATOM"):
            g[(line[21], line[22:27])].append(float(line[60:66]))
    return g

def bfactors_cif(path):
    """Same, for an AlphaFold3-style mmCIF (label_asym_id, label_seq_id, B_iso)."""
    g = collections.defaultdict(list)
    for line in open(path):
        if line.startswith("ATOM"):
            f = line.split()
            g[(f[6], f[8])].append(float(f[14]))
    return g

def audit(path, reader):
    g = reader(path)
    flat = [v for vals in g.values() for v in vals]
    spread = max(max(v) - min(v) for v in g.values())
    plausible = max(flat) <= 100.0 and min(flat) >= 0.0 and statistics.median(flat) > 50.0
    print(f"{path:22s} residues {len(g):5d}  range {min(flat):6.2f}-{max(flat):6.2f}  "
          f"median {statistics.median(flat):6.2f}  max within-residue spread {spread:5.2f}  "
          f"looks like pLDDT: {plausible}")

audit("egfr_structures/aureliabustos.bce_var433.pdb", bfactors_pdb)   # AlphaFold2 PDB
audit("fold_aurka_0_tpx2_0_model_0.cif", bfactors_cif)                # AlphaFold3 mmCIF
audit("1UBQ.pdb", bfactors_pdb)                                       # X-ray monomer
audit("1BJ1.pdb", bfactors_pdb)                                       # X-ray complex

# and what "the residue's pLDDT" means when the values are genuinely per-atom
cif = bfactors_cif("fold_aurka_0_tpx2_0_model_0.cif")
ca = {}
for line in open("fold_aurka_0_tpx2_0_model_0.cif"):
    if line.startswith("ATOM"):
        f = line.split()
        if f[3] == "CA":
            ca[(f[6], f[8])] = float(f[14])
chain_a = [k for k in ca if k[0] == "A"]
print(f"AF3 chain A mean pLDDT: {statistics.mean(ca[k] for k in chain_a):.2f} from CA, "
      f"{statistics.mean(statistics.mean(cif[k]) for k in chain_a):.2f} from per-residue "
      f"atom means")

Four things this prints, and each of them is a bug someone has shipped:

• The B-factor column of a deposited crystal structure is a crystallographic B-factor, not pLDDT. 1UBQ runs 2.00–42.75 and 1BJ1 reaches 142.29. Both are perfectly ordinary; neither is a confidence score, and for a B-factor low is good, the opposite direction. A value above 100 proves the column is not pLDDT. A median below 50 is a strong hint. Barnase–barstar, the tightest protein–protein complex known, scores a mean "pLDDT" of 26 if you make this mistake.
• AlphaFold2 writes one pLDDT per residue, repeated on every atom; AlphaFold3 writes a genuine per-atom value — its own documentation calls pLDDT "a per-atom confidence estimate". The audit prints a within-residue spread of 0.00 for the AF2 files and 47.10 for the AF3 mmCIF. So CA is the residue's pLDDT in AF2 and only one of its atoms in AF3 — the last line of the block prints the consequence, the same chain reading 94.68 from CA and 92.18 from per-residue atom means.
• Averaging over atoms instead of residues silently weights tryptophans over glycines. Across the 400 EGFR design structures the per-atom binder mean differs from the per-residue binder mean by 0.92 pLDDT on average.
• A modified residue is a HETATM record even in a predicted structure. The audit above reports 317 residues for a chain pair that has 319, because the two phosphothreonines of this AURKA model are HETATM and an ATOM-only parse drops them. Harmless in a confidence audit; not harmless if one of them is in your interface.

Filter 1 — interface confidence, from the PAE matrix

ipTM is a global number, and it is not the interface

ipTM is one scalar for the whole prediction. It does not tell you which interface in a multi-chain model is real, and per-chain-pair variants do not fix the ranking.

import gzip, json, itertools
import numpy as np

scores = json.load(gzip.open(
    "RAF1_KSR1_MEK1_9f755_scores_alphafold2_multimer_v3_model_1_seed_000.json.gz"))
pae = np.asarray(scores["pae"], dtype=float)
chains = np.array([l[21] for l in
                   open("RAF1_KSR1_MEK1_9f755_unrelaxed_alphafold2_multimer_v3_model_1_seed_000.pdb")
                   if l.startswith("ATOM") and l[12:16].strip() == "CA"])
xyz = np.array([[float(l[30:38]), float(l[38:46]), float(l[46:54])] for l in
                open("RAF1_KSR1_MEK1_9f755_unrelaxed_alphafold2_multimer_v3_model_1_seed_000.pdb")
                if l.startswith("ATOM")])
atom_chains = np.array([l[21] for l in
                        open("RAF1_KSR1_MEK1_9f755_unrelaxed_alphafold2_multimer_v3_model_1_seed_000.pdb")
                        if l.startswith("ATOM")])

print(f"one global ipTM for the whole model: {scores['iptm']}")
print(f"pairwise ipTM reported by the predictor: {scores['pairwise_iptm']}")
for a, b in itertools.combinations(sorted(set(chains)), 2):
    A, B = xyz[atom_chains == a], xyz[atom_chains == b]
    n = sum(int((np.linalg.norm(A[i:i+500, None, :] - B[None, :, :], axis=-1) < 5.0).sum())
            for i in range(0, len(A), 500))
    print(f"  chains {a}-{b}: {n:5d} heavy-atom contacts under 5 A")

The global ipTM of 0.53 would fail a "confident interface" threshold of 0.6 even though this model contains two well-packed interfaces. The pairwise values invert the ordering: B-C, with 170 contacts, scores 0.604 — above A-B, which has 1254.

Use ipTM to reject nothing. Compute a per-pair score from the PAE matrix instead.

ipAE, ipSAE and the direction that matters

Three scores are computed from the same cross-chain block of the PAE matrix, and they differ in how they weight it.

• ipAE (often pae_interaction): the mean PAE over residue pairs with one residue in each chain. Simple, and dominated by the many pairs far from the interface.
• pDockQ / pDockQ2: logistic fits over interface pLDDT and contact count, or over PAE at contacting pairs. Calibrated to DockQ, so they answer "is this pose right", not "does it bind".
• ipSAE: for each residue i in the aligned chain, take only the partner residues with PAE below a cutoff, set d0 from how many those are, and average 1/(1+(PAE/d0)²) over them; the chain-pair score is the maximum over i. Restricting to confident pairs is what stops a big, badly-predicted chain from washing the interface out.

import json
import numpy as np

def load_af3_pae(json_path, protein_chains):
    """AlphaFold3 PAE is per TOKEN, and a protein chain has more tokens than residues.
    Collapse to one row/column per residue before doing anything else."""
    d = json.load(open(json_path))
    pae = np.asarray(d["pae"], dtype=float)
    ch, rid = d["token_chain_ids"], d["token_res_ids"]
    keep, seen = [], set()
    for i, (c, r) in enumerate(zip(ch, rid)):
        if c in protein_chains and (c, r) not in seen:
            seen.add((c, r))
            keep.append(i)
    keep = np.array(keep)
    return pae[np.ix_(keep, keep)], np.array([ch[i] for i in keep])

def d0(n):
    """Yang and Skolnick length normalisation, floored at 1.0 as ipSAE does."""
    n = np.asarray(n, dtype=float)
    return np.maximum(1.0, 1.24 * np.sign(n - 15) * np.abs(n - 15) ** (1 / 3) - 1.8)

def ipsae_asym(pae, chains, aligned, scored, cutoff=10.0):
    """ipSAE for aligned -> scored. Asymmetric on purpose."""
    rows, cols = chains == aligned, chains == scored
    ok = np.zeros_like(pae, dtype=bool)
    ok[np.ix_(rows, cols)] = pae[np.ix_(rows, cols)] < cutoff
    n0 = ok.sum(axis=1)
    dd = d0(n0)
    best = 0.0
    for i in np.where(rows)[0]:
        m = ok[i]
        if m.any():
            best = max(best, float((1.0 / (1.0 + (pae[i, m] / dd[i]) ** 2)).mean()))
    return best

def ipae(pae, chains, a, b):
    return float(pae[np.ix_(chains == a, chains == b)].mean())

pae, chains = load_af3_pae("fold_aurka_0_tpx2_0_full_data_0.json", {"A", "B"})
ab = ipsae_asym(pae, chains, "A", "B")
ba = ipsae_asym(pae, chains, "B", "A")
print(f"residues after collapsing tokens: {len(chains)}  "
      f"(A {(chains=='A').sum()}, B {(chains=='B').sum()})")
print(f"ipSAE A->B (align on the 276-residue chain, score the 43-residue one): {ab:.6f}")
print(f"ipSAE B->A (the other direction):                                     {ba:.6f}")
print(f"ipSAE_max {max(ab, ba):.6f}   ipSAE_min {min(ab, ba):.6f}")
print(f"mean interface PAE (ipAE): {ipae(pae, chains, 'A', 'B'):.3f} A")
print(open("fold_aurka_0_tpx2_0_model_0_10_10.txt").read().rstrip())

Two things to take from that output.

The implementation agrees with the reference. The published values for this model are 0.448952 and 0.866498; the code above returns 0.448952 and 0.866531. The last digits differ only because a residue tokenised more than once has more than one PAE row and the two implementations pick different ones.

That collapse is not cosmetic. This file carries 296 tokens for chain A's 276 residues — AlphaFold3 tokenises modified residues one token per atom, and this AURKA model is phosphorylated on Thr160 and Thr161, giving 11 tokens each. Index the PAE matrix by token and B→A reads 0.874849 instead of 0.866531. The error grows with how much of your interface is modified, and nothing in the file warns you.

The two directions differ by a factor of nearly two on the same interface, and this is mechanical rather than mysterious. d0 is set from the number of confident partner residues, so the direction that scores the small chain gets a small d0 and a harsh score. In binder design the binder is always the small chain, so align on the target and score the binder — which is the min of the two directions for any realistic target/binder size ratio. That is the conservative convention, and the one benchmarked in 2026 across 3,532 designs.

It is also, on this particular complex, wrong: DockQ against the crystal structure (below) says the model is right, and the pessimistic direction reads 0.449. Conservative filters throw away good designs. That is what they are for; you just have to know the rate.

Where the numbers come from

filter	threshold	source	what sits either side
interface PAE (ipAE)	< 10 Å	Bennett et al. 2023, Nat Commun 14:2625 — the AF2 initial-guess filter, paired with af2_complex_rmsd < 5 Å	~2.3% of raw designs pass; on the EGFR benchmark it keeps 105/378 and 24/53 binders
interface PAE, normalised	i_pAE ≤ 0.35	BindCraft settings_filters/default_filters.json (MIT). ColabDesign divides PAE by 31, so this is ≈10.9 Å raw	the same filter as the row above, in different units
ipTM	> 0.8 confident, < 0.6 failed, between is a grey zone	AlphaFold 3 docs/output.md, verbatim	global, so it cannot rank interfaces within a model
i_pTM	≥ 0.5	BindCraft default	looser than the row above and, on the EGFR benchmark, almost non-selective — it keeps 280 of 378
binder pLDDT	≥ 80 (0.8 on BindCraft's 0–1 scale)	BindCraft default; AlphaFold DB bands are >90 very high, 70–90 confident, 50–70 low	keeps 222/378 and 39/53 binders
ipSAE	no published cutoff; rank, do not threshold	Dunbrack 2025 (IPSAE, MIT); Overath et al. 2025 computed >200 features per design and found the AF3-derived ipSAE beat both ipAE and ipTM, at 1.4× the average precision of ipAE	the score scales with interface size, so a cutoff transferred between targets is not the same filter
design vs re-prediction RMSD	≤ 3.5 Å binder, < 5 Å complex	BindCraft default; Bennett et al. 2023	see Filter 3
DockQ	≥ 0.23 acceptable, ≥ 0.49 medium, ≥ 0.80 high	printed by DockQ v2 itself	a quality band for a pose, not a binding probability
shape complementarity Sc	≥ 0.55 per model, ≥ 0.60 averaged	BindCraft default; Lawrence and Colman 1993, J Mol Biol 234:946	needs Rosetta — see What needs a licence
surface hydrophobicity	≤ 0.35	BindCraft default (exposed apolar residues over exposed residues)	a developability filter, not a binding one
interface residues / H-bonds	≥ 7 / ≥ 3, unsatisfied ≤ 4	BindCraft defaults	see Filter 4

Filter 2 — pLDDT, scoped to the right chain

The target is usually large and always well predicted, so a complex-wide mean pLDDT is a statement about the target. Score the binder, and score the binder's interface.

import collections
import numpy as np

def per_residue_plddt(path):
    """AF2 PDB: one value per residue, repeated across its atoms. Take the CA."""
    out = collections.defaultdict(dict)
    for line in open(path):
        if line.startswith("ATOM") and line[12:16].strip() == "CA":
            out[line[21]][line[22:27]] = float(line[60:66])
    return out

def interface_residues(path, binder="A", target="B", cutoff=5.0):
    coords = collections.defaultdict(list)
    for line in open(path):
        if line.startswith("ATOM"):
            coords[(line[21], line[22:27])].append(
                (float(line[30:38]), float(line[38:46]), float(line[46:54])))
    tgt = np.array([c for (ch, _), v in coords.items() if ch == target for c in v])
    hits = []
    for (ch, res), v in coords.items():
        if ch != binder:
            continue
        d = np.linalg.norm(np.array(v)[:, None, :] - tgt[None, :, :], axis=-1)
        if (d < cutoff).any():
            hits.append(res)
    return hits

for name in ["round1zeroshot.K5Q_N70S_K71R_N73T_S87T_N88D_R179K_K183R_E213D_S214P",
             "chrisxushaoyong.hu_nano2_4_85252b",
             "elian.elian2"]:
    p = f"egfr_structures/{name}.pdb"
    pl = per_residue_plddt(p)
    iface = interface_residues(p)
    binder = np.array(list(pl["A"].values()))
    both = np.array(list(pl["A"].values()) + list(pl["B"].values()))
    print(f"{name[:34]:36s} binder {binder.mean():6.2f}  complex {both.mean():6.2f}  "
          f"interface {np.mean([pl['A'][r] for r in iface]):6.2f}  "
          f"binder min {binder.min():6.2f}")

On the 378 labelled EGFR designs, ranked by area under the ROC curve against the binder/non-binder call:

score	AUROC	average precision
interface pLDDT, binder side	0.684	0.226
binder mean pLDDT	0.656	0.217
ipTM	0.636	0.207
interface PAE	0.612	0.210
complex mean pLDDT	0.609	0.220
ESM2 pseudo-log-likelihood	0.547	0.212
binder minimum pLDDT	0.439	0.129

Base rate 0.140. Two readings, both useful. Scoping pLDDT to the interface beats scoping it to the binder, which beats the complex mean — as much as any of these beats anything. And the minimum pLDDT over the binder is worse than chance: a design with one flexible terminus is not a worse binder, and filtering on a minimum removes long designs rather than bad ones.

Filter 3 — self-consistency

Bennett et al. split design failure in two. Type 1: the sequence does not fold to the intended monomer. Type 2: it folds and does not form the intended interface. Confidence scores speak to neither directly — a model can be confidently wrong. Self-consistency does: fold the sequence again, independently, and ask whether you get the structure you designed.

DockQ is the standard measure, and it separates the two failures — fnat is interface recovery, LRMSD is the pose, iRMSD is the interface geometry.

./.venv/bin/DockQ fold_aurka_0_tpx2_0_model_0.cif 1OL5.pdb 2>&1 | grep -vE '^\*|^