Bio Structural Biology Structure Preparation

Version Compatibility

Reference examples tested with: pdbfixer 1.9+, openmm 8.1+

reduce/Reduce2, propka3, and pdb2pqr30 are external command-line tools installed separately (conda install -c conda-forge reduce propka pdb2pqr or their own pip/binary packages); the code calls them via subprocess.

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Structure Preparation

"Make this receptor simulation-ready - add hydrogens and fix the protonation" -> add the atoms the experiment never resolved, assign pH- and environment-dependent protonation/tautomer states, and fill modeled-absent atoms and short loops.

• Python: pdbfixer.PDBFixer for missing atoms/residues + hydrogens at a stated pH; openmm.app.Modeller.addHydrogens for force-field-ready H with explicit variants
• CLI: reduce/mmtbx.reduce2 for H-bond-network flip and His tautomer optimization; propka3 for pKa prediction; pdb2pqr30 --with-ph --ff=AMBER --titration-state-method propka for pKa + PQR

Governing Principle

A deposited or predicted structure is NOT simulation-ready, and "preparing" it stacks a second layer of inference on top of a model that was already inferred from data. Every atom added, every proton placed, and every loop built is a hypothesis, not a measurement - so the prepared file is a distinct object from the deposited one, and the assumptions (pH, protonation model, what was built) must travel with it or downstream H-bond, salt-bridge, interface, and docking results become unreproducible.

Hydrogens are the first trap. X-ray crystallography rarely resolves them: a hydrogen carries ~1 electron and only becomes visible in density at roughly sub-1.2 Angstrom resolution, so almost every crystal structure arrives with zero (or only polar) hydrogens. They must be ADDED, and where a hydrogen goes is decided by the PROTONATION and TAUTOMER state of its residue, which crystallographic density alone frequently cannot distinguish. His has three relevant states (HID/HIE/HIP: proton on ND1, on NE2, or both/charged); Asn, Gln, and His side-chain amide/ring groups each have a 180-degree FLIP that swaps look-alike atoms (O for N, N for C) that density at typical resolution cannot tell apart (Word et al. 1999 J Mol Biol 285:1735). Placing a proton "at standard pKa 7" is wrong for any residue whose microenvironment shifts its pKa: a buried Asp/Glu can stay protonated well above pH 7, a Cys in a catalytic or metal site can be a deprotonated thiolate, a Lys buried near acidic residues can lose its charge. Use a pKa predictor (PROPKA, H++) for titratable residues and an all-atom H-bond-network optimizer (reduce) for flips and His tautomers - never trust the standard state for anything not freely solvated.

Missing atoms and short missing LOOPS are the second trap. A residue truncated to Cbeta, or a chain that jumps 45 -> 58, almost always means the region was DISORDERED (too mobile to model into density), not deleted - the atoms exist in reality (see structure-navigation, structure-validation). Building them back is legitimate for making a system topologically complete, but a built side chain or loop is a GUESS among many possible conformations, must be flagged as such, and must never be reported as experimental. Long gaps, terminal extensions, and anything spanning a domain are beyond what a preparation tool should invent - hand those to a modeling method, not PDBFixer.

Decision: tool by preparation task

Decision: standard protonation state vs pKa predictor

The one-line rule: standard states are defensible only for residues in a bulk-solvent-like environment; any titratable residue that is buried, charged-clustered, metal-adjacent, or in a pocket needs a predictor (PROPKA/H++) and an H-bond-network pass (reduce), with the chosen pH stated.

PDBFixer preparation pipeline

Goal: Turn a raw PDB/mmCIF into a hydrogen-complete, gap-filled structure at a chosen pH, recording what was added.

Approach: Call the PDBFixer finders in their required order (missing residues, then nonstandard, then missing atoms) so addMissingAtoms sees both sets; strip crystallization heterogens while keeping (or dropping) water deliberately; then add hydrogens at an explicitly chosen pH. The order matters - hydrogens are added last, after heavy atoms exist.

from pdbfixer import PDBFixer
from openmm.app import PDBFile

fixer = PDBFixer(filename='receptor.pdb')  # or PDBFixer(pdbid='1VII') to fetch from RCSB

fixer.findMissingResidues()          # short internal gaps + terminals, from SEQRES vs modeled
fixer.findNonstandardResidues()      # e.g. MSE (selenomethionine), modified residues
fixer.replaceNonstandardResidues()   # map them back to standard parents
fixer.removeHeterogens(keepWater=False)  # drop buffer ions/cryoprotectants; keepWater=True to retain
fixer.findMissingAtoms()             # truncated side chains + the residues found above
fixer.addMissingAtoms()              # build heavy atoms; built loops are HYPOTHESES, log them

# pH 7.0 is a CHOICE, not a safe default - state it and match the experimental/biological condition.
# addMissingHydrogens detects existing disulfides and leaves those Cys as CYX (no SG hydrogen).
fixer.addMissingHydrogens(pH=7.0)

with open('receptor_prepared.pdb', 'w') as out:
    PDBFile.writeFile(fixer.topology, fixer.positions, out, keepIds=True)

# missingResidues is a dict {(chain_index, residue_index): [resname, ...]} - flatten it for provenance.
built = [(ci, pos, name) for (ci, pos), names in fixer.missingResidues.items() for name in names]
print(f'built {len(fixer.missingResidues)} missing-residue segment(s), {len(built)} residue(s); pH=7.0')

Optimize flips and His tautomers with reduce

Goal: Resolve Asn/Gln/His amide and ring orientations and His protonation by all-atom H-bond-network scoring, which density at typical resolution cannot settle.

Approach: Run reduce with building enabled so it adds hydrogens AND evaluates the 180-degree flip of each Asn/Gln/His plus His NH placement, choosing the orientation that optimizes the local hydrogen-bond network and minimizes clashes. reduce optimizes GEOMETRY, not titration - pair it with a pKa predictor for charge states.

import subprocess

# -build runs -OH -ROTEXOH -HIS -FLIP: adds H and optimizes OH/His rotation plus Asn/Gln/His flips
# (a superset of -FLIP, not an alias). reduce scores orientation by small-probe all-atom contacts (Word 1999).
with open('receptor_reduced.pdb', 'w') as out:
    subprocess.run(['reduce', '-build', 'receptor_prepared.pdb'], stdout=out, check=True)

# Reduce2 (CCTBX/Phenix) is the maintained successor: mmtbx.reduce2 receptor_prepared.pdb

Predict pKa and write a PQR for electrostatics

Goal: Decide which titratable residues deviate from standard states at a target pH, and emit charges+radii for Poisson-Boltzmann electrostatics.

Approach: PROPKA predicts per-residue pKa from the 3D environment; PDB2PQR wraps PROPKA to assign protonation at --with-ph, add hydrogens for the chosen force field, and write a PQR. Feed the PQR to APBS for the electrostatic potential (a separate downstream step; see the electrostatics note below).

import subprocess

# propka3 writes receptor_prepared.pka; the SUMMARY lists predicted pKa vs model (standard) pKa.
subprocess.run(['propka3', 'receptor_prepared.pdb'], check=True)

# PDB2PQR assigns states at pH via PROPKA and writes charges/radii for the named force field.
# --ff must MATCH the downstream force field; --with-ph 7.0 is the stated titration condition.
subprocess.run([
    'pdb2pqr30', '--ff=AMBER', '--with-ph', '7.0',
    '--titration-state-method', 'propka', '--keep-chain',
    'receptor_prepared.pdb', 'receptor.pqr',
], check=True)

Electrostatics note: the PQR is the input to APBS (Jurrus et al. 2018 Protein Sci 27:112) for the Poisson-Boltzmann potential/surface; PDB2PQR can emit an APBS input file. Keep force field, pH, and radii set identical between preparation and the APBS run.