Free Energy Perturbation (FEP) for Drug Discovery

FEP transforms ligand A into ligand B alchemically — atoms are "mutated" via a coupling parameter λ that scales off atoms of A and scales on atoms of B. Simulations at multiple λ windows provide overlap between adjacent states, and the Bennett Acceptance Ratio (BAR) or Multistate BAR (MBAR) estimator combines window ΔG values into ΔΔG​bind with optimal statistical efficiency.

The thermodynamic cycle closure property is fundamental: ΔG​bind(B) − ΔG​bind(A) = ΔG​pert(complex) − ΔG​pert(solvent). This means chemical accuracy requires running both the complex and solvated-ligand legs of each transformation.

Thermodynamic Cycle and RBFE Setup

import numpy as np
from dataclasses import dataclass, field
from typing import Optional


@dataclass
class FEPEdge:
    """Represents one alchemical transformation A → B."""
    ligand_a: str
    ligand_b: str
    n_lambda_windows: int = 12
    simulation_ns_per_window: float = 5.0
    complex_dG: Optional[float] = None     # kcal/mol
    solvent_dG: Optional[float] = None     # kcal/mol
    complex_dG_err: Optional[float] = None
    solvent_dG_err: Optional[float] = None

    @property
    def dG_bind(self) -> Optional[float]:
        """RBFE: ΔΔG_bind = ΔG_complex − ΔG_solvent (thermodynamic cycle)."""
        if self.complex_dG is not None and self.solvent_dG is not None:
            return self.complex_dG - self.solvent_dG
        return None

    @property
    def dG_bind_err(self) -> Optional[float]:
        if self.complex_dG_err is not None and self.solvent_dG_err is not None:
            return np.sqrt(self.complex_dG_err**2 + self.solvent_dG_err**2)
        return None


@dataclass
class FEPNetwork:
    """
    Perturbation network of ligand pairs for relative binding free energy (RBFE).

    Best practice: use a connected graph where each ligand appears in ≥2 edges,
    allowing cycle closure correction and detection of poorly converged edges.
    """
    edges: list[FEPEdge] = field(default_factory=list)
    reference_ligand: Optional[str] = None
    experimental_IC50: dict[str, float] = field(default_factory=dict)

    def add_edge(self, ligand_a: str, ligand_b: str, **kwargs) -> FEPEdge:
        edge = FEPEdge(ligand_a=ligand_a, ligand_b=ligand_b, **kwargs)
        self.edges.append(edge)
        return edge

    def compute_cycle_closure_errors(self) -> list[dict]:
        """
        Find closed cycles in edge graph and compute cycle closure error (CCE).

        CCE = |sum of ΔΔG around a cycle|; should approach 0 for converged sims.
        CCE > 1 kcal/mol indicates at least one edge has convergence issues.
        """
        from itertools import combinations
        closure_errors = []
        edge_dict = {}
        for e in self.edges:
            if e.dG_bind is not None:
                edge_dict[(e.ligand_a, e.ligand_b)] = e.dG_bind
                edge_dict[(e.ligand_b, e.ligand_a)] = -e.dG_bind

        # Find triangular cycles
        ligands = list({l for e in self.edges for l in (e.ligand_a, e.ligand_b)})
        for a, b, c in combinations(ligands, 3):
            if all(k in edge_dict for k in [(a,b),(b,c),(c,a)]):
                cce = edge_dict[(a,b)] + edge_dict[(b,c)] + edge_dict[(c,a)]
                closure_errors.append({"cycle": (a, b, c), "cce_kcal_mol": round(cce, 3)})

        return closure_errors

MBAR Estimator for Multi-State FEP

import numpy as np


def compute_bar_estimate(
    dU_forward: np.ndarray,
    dU_reverse: np.ndarray,
    kT: float = 0.592,
    n_iterations: int = 100,
) -> tuple[float, float]:
    """
    Bennett Acceptance Ratio (BAR) estimate for ΔG between adjacent lambda windows.

    BAR minimises the variance of the free energy estimate by finding the
    optimal overlap between forward (A→B) and reverse (B→A) energy distributions.

    Args:
        dU_forward: Potential energy differences U_B - U_A sampled in state A (kJ/mol or kcal/mol)
        dU_reverse: Potential energy differences U_A - U_B sampled in state B
        kT: Thermal energy (0.592 kcal/mol at 298 K)
        n_iterations: Self-consistency iterations for BAR equation

    Returns:
        (delta_G, delta_G_std) in same units as dU.
    """
    N_f, N_r = len(dU_forward), len(dU_reverse)
    M = np.log(N_f / N_r)
    dG = 0.0

    for _ in range(n_iterations):
        f_fwd = 1.0 / (1.0 + np.exp(+(dU_forward / kT - dG / kT - M)))
        f_rev = 1.0 / (1.0 + np.exp(-(dU_reverse / kT + dG / kT + M)))
        dG_new = kT * (np.log(np.mean(f_fwd)) - np.log(np.mean(f_rev))) + dG
        if abs(dG_new - dG) < 1e-6:
            break
        dG = dG_new

    # BAR statistical uncertainty (Shirts 2008)
    var = kT**2 * (
        1.0 / (N_f * np.mean(f_fwd**2)) + 1.0 / (N_r * np.mean(f_rev**2))
    )
    std = np.sqrt(max(var, 0.0))
    return round(float(dG), 4), round(float(std), 4)


def accumulate_fep_profile(
    window_dGs: list[float],
    window_stds: list[float],
) -> tuple[float, float]:
    """
    Sum lambda-window ΔGs to get total ΔG across the alchemical transformation.

    Args:
        window_dGs: ΔG for each adjacent lambda pair
        window_stds: Statistical uncertainty for each window pair

    Returns:
        (total_dG, total_std) with properly propagated error.
    """
    total_dG = sum(window_dGs)
    total_std = np.sqrt(sum(s**2 for s in window_stds))
    return round(total_dG, 3), round(total_std, 3)

FEP Convergence Diagnostics

def check_fep_convergence(
    fwd_energies: np.ndarray,
    rev_energies: np.ndarray,
    kT: float = 0.592,
) -> dict:
    """
    Diagnose FEP convergence via overlap matrix metrics and hysteresis analysis.

    Three key checks:
    1. Phase-space overlap: Bhattacharyya coefficient of forward/reverse distributions
    2. Forward-reverse hysteresis: |ΔG_fwd + ΔG_rev| < 0.5 kcal/mol for good convergence
    3. Normalised free energy timeseries drift (requires blocking analysis externally)

    Args:
        fwd_energies: U_B - U_A sampled in state A
        rev_energies: U_A - U_B sampled in state B
        kT: Thermal energy factor

    Returns:
        Convergence diagnostics dict.
    """
    # Bhattacharyya overlap between Boltzmann-weighted distributions
    combined = np.concatenate([fwd_energies, -rev_energies])
    bins = np.linspace(combined.min(), combined.max(), 50)
    h_fwd, _ = np.histogram(fwd_energies, bins=bins, density=True)
    h_rev, _ = np.histogram(-rev_energies, bins=bins, density=True)
    overlap = float(np.sum(np.sqrt(h_fwd * h_rev + 1e-12)) * (bins[1] - bins[0]))

    # One-way estimates for hysteresis
    dG_fwd = -kT * np.log(np.mean(np.exp(-fwd_energies / kT)))
    dG_rev = -kT * np.log(np.mean(np.exp(-rev_energies / kT)))
    hysteresis = abs(dG_fwd - dG_rev)

    diagnosis = "converged" if overlap > 0.3 and hysteresis < 1.0 else (
        "marginal" if overlap > 0.15 or hysteresis < 2.0 else "unconverged"
    )

    return {
        "phase_space_overlap": round(overlap, 3),
        "hysteresis_kcal_mol": round(hysteresis, 3),
        "dG_forward_kcal_mol": round(dG_fwd, 3),
        "diagnosis": diagnosis,
        "recommendation": (
            "Increase simulation time or add lambda windows" if diagnosis != "converged"
            else "Acceptable convergence"
        ),
    }

FEP Method Comparison

Key Takeaways

• RBFE via FEP is currently the most accurate computational method for lead optimisation rank-ordering within a congeneric series, achieving ~0.8–1.2 kcal/mol RMSE when the force field, charge model, and binding mode preparation are correct.
• The thermodynamic cycle mandates running both complex and solvent legs of each alchemical transformation. Cycle closure error (CCE) diagnosis — checking that ΔΔG around closed triangles sums near zero — is the best practical convergence test.
• Phase-space overlap between adjacent lambda windows is the key sampling diagnostic. Bhattacharyya coefficient >0.3 or forward-reverse hysteresis <1.0 kcal/mol indicates sufficient overlap. If overlap is poor, add intermediate lambda windows.
• FEP fails for large perturbations (>6 heavy atom changes, ring-size changes, net charge changes). For scaffold hops, use ABFE in separate legs or fall back to ML-based ΔΔG predictors (e.g. δ-ML corrections on top of RBFE).

Try it live

Explore a related model hands-on: open Molecule Screening Engine — Similarity Search Live and change its parameters while it runs. Everything stays in your browser.