Protein Structure and Binding Site Analysis

Structure-based drug design (SBDD) relies on knowing the three-dimensional shape of a protein target's binding site. A drug-like small molecule must fit into the pocket like a key into a lock — making specific interactions with residues that are important for the protein's function. Until 2021, the primary bottleneck was obtaining a protein structure: X-ray crystallography takes months and can fail for membrane proteins or flexible targets. AlphaFold2 has substantially changed this: high-confidence structures are now available within hours for most human protein targets, enabling routine SBDD even for previously intractable targets.

This guide covers obtaining and validating protein structures, detecting binding pockets computationally, scoring their druggability, and preparing structures for docking pipelines.

Protein Structure Sources

Before running AlphaFold2, always check whether an experimental structure already exists:

Retrieving Structures Programmatically

import requests
import json
import os

def fetch_pdb_structure(pdb_id: str, output_dir: str = ".") -> str:
    """
    Download a PDB structure file from the RCSB PDB.
    Returns the path to the downloaded file.
    """
    url = f"https://files.rcsb.org/download/{pdb_id.upper()}.pdb"
    response = requests.get(url, timeout=30)
    response.raise_for_status()

    filepath = os.path.join(output_dir, f"{pdb_id.lower()}.pdb")
    with open(filepath, "w") as f:
        f.write(response.text)
    print(f"Downloaded {pdb_id}: {filepath}")
    return filepath


def fetch_alphafold_structure(uniprot_id: str, output_dir: str = ".") -> str:
    """
    Download an AlphaFold2 predicted structure from the AlphaFold DB.
    Returns the path to the downloaded PDB file.
    """
    url = f"https://alphafold.ebi.ac.uk/files/AF-{uniprot_id}-F1-model_v4.pdb"
    response = requests.get(url, timeout=30)
    if response.status_code == 404:
        raise ValueError(f"No AlphaFold structure found for {uniprot_id}")
    response.raise_for_status()

    filepath = os.path.join(output_dir, f"AF_{uniprot_id}.pdb")
    with open(filepath, "w") as f:
        f.write(response.text)
    print(f"Downloaded AlphaFold structure for {uniprot_id}: {filepath}")
    return filepath


def get_alphafold_confidence(uniprot_id: str) -> dict:
    """
    Retrieve per-residue pLDDT confidence scores from the AlphaFold API.
    High pLDDT (>70) indicates reliable structural regions.
    """
    url = f"https://alphafold.ebi.ac.uk/api/prediction/{uniprot_id}"
    response = requests.get(url, timeout=30)
    response.raise_for_status()
    data = response.json()
    if not data:
        raise ValueError(f"No prediction data for {uniprot_id}")
    return data[0]  # First model entry


# Example
pdb_path = fetch_alphafold_structure("P00533")  # EGFR kinase domain
conf = get_alphafold_confidence("P00533")
print(f"Mean pLDDT: {sum(conf['plddt']) / len(conf['plddt']):.1f}")

Understanding AlphaFold2 Confidence Scores

AlphaFold2 outputs a per-residue confidence score called pLDDT (predicted Local Distance Difference Test), ranging 0–100:

For docking, only use regions with pLDDT > 70. Low-confidence loops often flank binding pockets and can create artefactual docking poses.

import numpy as np
from Bio import PDB

def analyze_pocket_residue_confidence(pdb_path: str, pocket_residues: list[int]) -> dict:
    """
    Check AlphaFold pLDDT confidence for binding site residues.
    In AlphaFold PDB files, pLDDT is stored in the B-factor column.

    Args:
        pdb_path:         Path to AlphaFold PDB file
        pocket_residues:  List of residue sequence numbers in the binding site

    Returns:
        Dict with residue -> pLDDT value
    """
    parser = PDB.PDBParser(QUIET=True)
    structure = parser.get_structure("protein", pdb_path)
    model = next(structure.get_models())
    chain = next(model.get_chains())

    confidence = {}
    for residue in chain.get_residues():
        resnum = residue.get_id()[1]
        if resnum in pocket_residues:
            # B-factor stores pLDDT in AlphaFold PDB outputs
            b_factors = [atom.get_bfactor() for atom in residue.get_atoms()]
            confidence[resnum] = round(np.mean(b_factors), 1)

    return confidence

Binding Pocket Detection with fpocket

fpocket is the leading open-source tool for detecting void spaces on protein surfaces that could accommodate a small molecule. It uses Voronoi tessellation to identify alpha spheres (empty spheres that fit inside protein cavities) and clusters them into pockets, then scores each pocket for druggability.

import subprocess
import os
import pandas as pd
from pathlib import Path

def run_fpocket(pdb_path: str, output_dir: str = None) -> dict:
    """
    Run fpocket on a protein PDB file to detect binding pockets.

    Returns:
        Dict with list of detected pockets and their properties.
    """
    pdb_path = Path(pdb_path)
    if output_dir is None:
        output_dir = pdb_path.parent / (pdb_path.stem + "_out")

    result = subprocess.run(
        ["fpocket", "-f", str(pdb_path)],
        capture_output=True,
        text=True,
        cwd=str(pdb_path.parent),
    )
    if result.returncode != 0:
        raise RuntimeError(f"fpocket failed: {result.stderr}")

    # Parse the summary info file
    info_path = pdb_path.parent / (pdb_path.stem + "_out") / (pdb_path.stem + "_info.txt")
    pockets = parse_fpocket_info(str(info_path))
    return pockets


def parse_fpocket_info(info_path: str) -> list[dict]:
    """Parse fpocket output info file into a list of pocket dicts."""
    pockets = []
    current_pocket = {}

    with open(info_path) as f:
        for line in f:
            line = line.strip()
            if line.startswith("Pocket"):
                if current_pocket:
                    pockets.append(current_pocket)
                pocket_num = int(line.split()[1].rstrip(":"))
                current_pocket = {"pocket_id": pocket_num}
            elif ":" in line and current_pocket:
                key, _, value = line.partition(":")
                key = key.strip().lower().replace(" ", "_")
                try:
                    current_pocket[key] = float(value.strip())
                except ValueError:
                    current_pocket[key] = value.strip()

    if current_pocket:
        pockets.append(current_pocket)
    return pockets

Druggability Assessment

Not every detected pocket is druggable. A pocket is druggable if a drug-like small molecule can bind with sufficient affinity (< 10 μM) to modulate the protein's function. The key pocket properties that predict druggability are:

def score_pocket_druggability(pocket: dict) -> str:
    """
    Apply a simple rule-based druggability filter on fpocket output.
    Returns 'druggable', 'borderline', or 'undruggable'.
    """
    volume = pocket.get("volume_score", 0)
    drug_score = pocket.get("drug_score", 0)
    hydro = pocket.get("hydrophobicity_score", 0)

    if drug_score > 0.5 and volume > 300 and hydro > 0.3:
        return "druggable"
    elif drug_score > 0.2:
        return "borderline"
    else:
        return "undruggable"


# Example: rank pockets by druggability
pockets = run_fpocket("egfr_kinase.pdb")
pocket_df = pd.DataFrame(pockets)
pocket_df["druggability_class"] = pocket_df.apply(score_pocket_druggability, axis=1)
pocket_df = pocket_df.sort_values("drug_score", ascending=False)

print(pocket_df[["pocket_id", "drug_score", "volume_score", "druggability_class"]].head(10).to_string())

Protein Preparation for Docking

Raw PDB files require preparation before use in docking. Key steps:

• Remove non-protein molecules: delete crystallographic waters, detergents, and co-purification artifacts (keep structural waters if they mediate binding).
• Add hydrogens: PDB files rarely contain hydrogen atoms. Add them at physiological pH (7.4) using Reduce, OpenBabel, or rdkit/protonate3D.
• Assign protonation states: use Epik (commercial) or Dimorphite-DL (free) to predict ionisation states of histidine, aspartate, glutamate, lysine at target pH.
• Resolve missing loops: if the binding site loop is missing (common in crystal structures), model it with MODELLER or PDBFixer.
• Check for alternate conformations: many crystal structures have two side-chain conformers. Select the one with higher occupancy or the one compatible with known Inhibitors.

from pdbfixer import PDBFixer
from openmm.app import PDBFile

def prepare_receptor(input_pdb: str, output_pdb: str, ph: float = 7.4):
    """
    Prepare a protein structure for docking using PDBFixer:
    1. Add missing residues
    2. Add missing heavy atoms
    3. Add hydrogens at specified pH

    Requires: pip install pdbfixer openmm
    """
    fixer = PDBFixer(filename=input_pdb)

    # Find and add missing residues (loops)
    fixer.findMissingResidues()
    fixer.findNonstandardResidues()
    fixer.replaceNonstandardResidues()

    # Remove heterogens (waters, ligands, co-factors)
    fixer.removeHeterogens(keepWater=False)

    # Add missing heavy atoms
    fixer.findMissingAtoms()
    fixer.addMissingAtoms()

    # Add hydrogens at specified pH
    fixer.addMissingHydrogens(ph)

    # Save
    with open(output_pdb, "w") as f:
        PDBFile.writeFile(fixer.topology, fixer.positions, f)

    print(f"Prepared receptor saved to {output_pdb}")


prepare_receptor("raw_egfr.pdb", "egfr_prepared.pdb", ph=7.4)

Cryptic Pocket Detection

Cryptic pockets are binding sites that are not visible in the apo (unliganded) structure but open up upon ligand binding or during conformational dynamics. They are particularly important for allosteric drug discovery. Detection methods include:

• Molecular dynamics simulations: run MD on the apo protein and apply fpocket to each frame. Pockets that appear transiently are cryptic binding sites.
• Probe-mapping (SiteMap, FTMap): place small organic probes on the protein surface and energy-minimize. Clusters of probes identify hotspots.
• CryptoSite, AlloSigMA: ML-based tools that predict cryptic pockets from sequence and structure without MD.

import subprocess
import pandas as pd
from pathlib import Path
import re

def detect_cryptic_pockets_from_trajectory(
    trajectory_frames: list[str],  # List of PDB file paths from MD
    min_frames_present: int = 3,   # Must appear in at least N frames
) -> pd.DataFrame:
    """
    Run fpocket on each MD frame and identify pockets that appear
    in multiple frames but may not exist in the static structure.
    """
    pocket_counts = {}  # pocket_id -> list of drug_scores across frames

    for frame_pdb in trajectory_frames:
        try:
            pockets = run_fpocket(frame_pdb)
        except Exception:
            continue
        for pocket in pockets:
            pocket_id = f"pocket_{pocket['pocket_id']}"
            if pocket_id not in pocket_counts:
                pocket_counts[pocket_id] = []
            pocket_counts[pocket_id].append(pocket.get("drug_score", 0))

    results = []
    for pocket_id, scores in pocket_counts.items():
        if len(scores) >= min_frames_present:
            results.append({
                "pocket_id": pocket_id,
                "n_frames": len(scores),
                "mean_drug_score": round(sum(scores) / len(scores), 3),
                "max_drug_score": round(max(scores), 3),
            })

    return pd.DataFrame(results).sort_values("mean_drug_score", ascending=False)

Key Takeaways

• AlphaFold2 provides high-quality predicted structures for most human proteins. Only use regions with pLDDT > 70 for structure-based design; low-confidence loops around binding sites must be validated.
• fpocket is the standard open-source pocket detector; its drug_score output (SVM-based) is a reliable druggability signal — target pockets with drug_score > 0.5 and volume 300–1500 ų.
• Always prepare the receptor before docking: add hydrogens at pH 7.4, remove non-relevant heterogens, and repair missing loops flanking the binding site.
• Cryptic pockets accessible only in dynamic states can be uncovered by applying pocket detection to MD trajectory frames — this is key for allosteric drug design.
• When using an AlphaFold2 structure for docking, compare it to the nearest PDB homolog. If the predicted binding site geometry matches a known drugged homolog, confidence in the docking results is much higher.

Try it live

Explore a related model hands-on: open Multi-Chain Docking (AlphaFold 3 Logic) and change its parameters while it runs. Everything stays in your browser.