Free Energy Perturbation Calculations in Drug Discovery: Alchemical Transformations, Thermodynamic Cycles, and Binding Affinity Prediction

Introduction

Free energy perturbation (FEP) calculations represent a cornerstone of computational structure-based drug design, enabling the rigorous estimation of relative and absolute binding free energies for ligand-protein complexes [1, 2]. The theoretical foundation of FEP rests on statistical mechanics and the Zwanzig equation, which relates the free energy difference between two thermodynamic states to an ensemble average of the Boltzmann factor of their energy difference [3, 4]. In drug discovery contexts, FEP is employed to rank-order congeneric series of compounds, guide lead optimization, and predict resistance mutations in target proteins [5, 6, 7]. The method has been applied across diverse therapeutic areas including oncology, infectious disease, and metabolic disorders [5, 8, 9]. For veterinary applications, FEP can inform the design of species-selective therapeutics targeting pathogens such as Mycobacterium tuberculosis thymidylate kinase or viral proteases [1, 10].

Theoretical Foundations of Free Energy Perturbation

The Zwanzig Equation

The central relationship in FEP is the Zwanzig equation, also known as the exponential averaging formula [3, 4]. For a system evolving from state A to state B, the free energy difference is given by:

ΔG(A→B) = -kBT ln ⟨exp(-(UB - UA)/kBT)⟩A

where UA and UB are the potential energies of states A and B, kB is the Boltzmann constant, T is the temperature, and the angle brackets denote an ensemble average over configurations sampled from state A [3, 4]. This equation is exact in the limit of infinite sampling but suffers from convergence issues when the energy difference between states is large [3, 11]. To mitigate this, the transformation is typically divided into multiple intermediate windows or λ-states, each representing a linear combination of the Hamiltonians of the end states [3, 12, 13].

Alchemical Transformations

Alchemical transformations refer to the computational mutation of one chemical species into another within a simulation box, without requiring the physical synthesis of intermediate compounds [3, 2]. In relative binding free energy (RBFE) calculations, a ligand A is alchemically transformed into ligand B both in the solvated state and in the protein-bound state [14, 15, 6]. The difference between these two free energy changes yields the relative binding affinity:

ΔΔGbind = ΔGbound(A→B) - ΔGfree(A→B)

This approach cancels systematic errors arising from force field inaccuracies and incomplete sampling, provided the perturbations are chemically similar [15, 16, 17]. The alchemical pathway is parameterized by a coupling parameter λ, which smoothly scales the interactions of the perturbed atoms from those of the initial state to those of the final state [3, 12, 13].

Thermodynamic Cycles

Thermodynamic cycles provide the conceptual framework for decomposing complex free energy calculations into tractable sub-calculations [18, 19, 20]. In the context of absolute binding free energy (ABFE) calculations, a thermodynamic cycle connects the bound and unbound states of a ligand with a protein through alchemical decoupling steps [19, 21, 22]. The cycle typically involves:

1. Decoupling the ligand from the protein binding site in the bound state.
2. Decoupling the ligand from bulk solvent in the unbound state.
3. Applying restraints to the ligand in the binding site to prevent translational and rotational drift [21, 22].

The net free energy change around the closed cycle is zero, allowing the binding free energy to be computed as the difference between the decoupling free energies in the bound and unbound environments [18, 19, 22]. For RBFE, the thermodynamic cycle connects two ligands through their respective bound and unbound states, as described above [14, 15, 6].

Computational Methods and Algorithms

Bennett Acceptance Ratio (BAR)

The Bennett Acceptance Ratio (BAR) is an estimator for free energy differences that achieves lower variance than the Zwanzig equation by using information from both forward and reverse perturbations [3, 11]. BAR solves for the free energy difference ΔG that satisfies:

⟨f(UB - UA - C)⟩A = ⟨f(UA - UB + C)⟩B

where f(x) = 1/(1 + exp(x/kBT)) is the Fermi function and C = ΔG - kBT ln(nB/nA) with nA and nB representing the number of samples from each state [3, 11]. BAR is widely implemented in modern FEP workflows and provides robust estimates when the phase space overlap between adjacent λ-windows is adequate [3, 12, 11]. The multistate Bennett Acceptance Ratio (MBAR) extends this formalism to multiple states simultaneously, improving statistical efficiency [3, 11].

λ-Dynamics and Alchemical Intermediate Methods

λ-Dynamics is an alternative to the traditional stepwise λ-window approach, in which λ is treated as a dynamic variable that evolves during the simulation [16, 7, 23]. This method allows the system to spontaneously sample multiple alchemical states within a single trajectory, enabling efficient exploration of chemical space [16, 7]. Multiple molecule λ-dynamics extends this concept to concurrently perturb both the ligand and protein residues, facilitating the study of drug resistance mutations [7]. Recent advances include the use of chemical intermediates to bridge distant transformations, as implemented in methods such as PairMap and IMERGE-FEP [11, 24]. These approaches insert intermediate chemical structures to improve phase space overlap and convergence [11, 24].

Dual-LAO and State Function Corrections

The Dual-LAO (Dual Lambda Adaptive Optimization) method accelerates RBFE calculations by adaptively optimizing the λ-schedule and sampling protocol [3]. This approach reduces the number of required λ-windows while maintaining accuracy for both simple and complex transformations [3]. State function-based corrections provide a post-processing framework to correct systematic errors in large-scale FEP calculations, improving agreement with experimental data without requiring additional simulations [25]. These corrections leverage the fact that free energy is a state function, allowing errors to be decomposed and addressed systematically [25].

Workflow for FEP-Guided Drug Design

The typical workflow for applying FEP calculations in drug discovery involves several sequential steps, as illustrated in Figure 1.

flowchart TD
    A[Target Protein Structure], > B[Ligand Series Design]
    B, > C[Pose Prediction and Docking]
    C, > D[System Preparation and Parameterization]
    D, > E[Alchemical Transformation Setup]
    E, > F[λ-Window Simulation and Sampling]
    F, > G[Free Energy Estimation via BAR/MBAR]
    G, > H[Rank Ordering and Experimental Validation]
    H, > I{Agreement with Experiment?}
    I, >|Yes| J[Lead Optimization]
    I, >|No| K[Re-evaluate Poses or Force Fields]
    K, > C

Figure 1. Workflow for FEP-guided drug design. The process begins with a target protein structure and a series of proposed ligands. After pose prediction and system preparation, alchemical transformations are simulated across multiple λ-windows. Free energy differences are estimated using BAR or MBAR, and the resulting rank order is validated experimentally.

Input Pose Quality and System Preparation

The accuracy of FEP calculations is highly sensitive to the quality of the input protein-ligand complex structure [26, 17]. Benchmarking studies have demonstrated that using experimentally determined or high-quality predicted holo structures significantly improves FEP performance compared to apo structures or poorly docked poses [26, 17]. Machine learning-based structure prediction methods, such as AlphaFold-derived models, can be used when experimental structures are unavailable, but their accuracy for FEP applications depends on the conformational state of the binding site [26]. System preparation involves assigning force field parameters, solvating the complex in a water box, adding counterions, and equilibrating the system [14, 18, 27].

Sampling and Convergence

Adequate sampling of conformational space is critical for converged free energy estimates [28, 21, 29]. Enhanced sampling techniques, including selective integrated tempering sampling and replica exchange, can accelerate convergence by overcoming energy barriers [30, 29]. The convergence-adaptive roundtrip method dynamically adjusts simulation lengths based on the observed convergence of free energy estimates, reducing computational cost while maintaining accuracy [29]. Water sampling is particularly important for binding sites with buried water molecules, and divide-and-conquer approaches that enhance water sampling have been shown to improve ABFE calculations [21].

Resource Allocation and Automation

Automated workflows for FEP calculations have been developed to increase throughput and reduce user intervention [12, 28]. The ALCHEMD framework provides an accessible and automated pipeline for RBFE calculations, integrating system preparation, simulation execution, and analysis [12]. On-the-fly optimization of resource allocation distributes computational effort across λ-windows based on their individual convergence rates, improving overall efficiency [28]. GPU acceleration of molecular dynamics engines, such as GROMACS, has dramatically reduced the wall-clock time required for FEP simulations, enabling routine application in industrial settings [27].

Applications in Drug Discovery

Relative Binding Free Energy Calculations

RBFE calculations are the most widely used FEP application in drug discovery, enabling the rank ordering of congeneric ligand series [15, 6, 8]. The FEP+ protocol has been benchmarked extensively for predicting binding affinities of congeneric ligands toward soluble proteins, demonstrating good correlation with experimental data for well-behaved systems [15]. Applications include the design of dipeptidyl peptidase 8/9 inhibitors, phosphodiesterase-1 inhibitors, and integrin αvβ6 inhibitors [6, 8, 9]. Core flipping, where the central scaffold of a ligand is modified, presents a particular challenge for RBFE calculations, but λ-dynamics approaches have been developed to handle such transformations [16].

Absolute Binding Free Energy Calculations

ABFE calculations estimate the absolute binding affinity of a ligand to a protein, which is useful for evaluating novel chemotypes or fragments [19, 21, 22]. These calculations require careful handling of restraints to prevent the ligand from drifting out of the binding site during decoupling [21, 22]. The RED-E-function-based equilibrium parameter finder automates the selection of optimal restraint parameters, improving the reproducibility and accuracy of ABFE calculations [22]. ABFE has been applied to discover ROR1 inhibitors using microsecond-scale molecular dynamics simulations and metadynamics [19].

Membrane Proteins and Challenging Targets

FEP calculations for membrane proteins present additional challenges due to the anisotropic environment of the lipid bilayer and the slow dynamics of membrane-embedded domains [18]. Specialized protocols have been developed to handle ligand binding to G protein-coupled receptors, ion channels, and transporters [18]. These protocols require careful treatment of the membrane environment, including the use of appropriate lipid models and the application of restraints to maintain the protein orientation [18].

Integration with Machine Learning

Machine learning methods are increasingly integrated into FEP workflows to improve efficiency and accuracy [31, 32, 33]. Machine learning scoring functions can be used to rapidly screen large compound libraries, with FEP calculations reserved for the most promising candidates [31, 32]. Force field parameters can be fine-tuned to experimental free energy measurements using machine learning optimization, improving the accuracy of subsequent FEP predictions [33]. Augmented data approaches combine machine learning predictions with FEP results to narrow the gap between scoring functions and rigorous free energy methods [32].

Limitations and Best Practices

Force Field Accuracy

The accuracy of FEP calculations is fundamentally limited by the quality of the underlying molecular mechanics force field [4, 33]. Errors in partial charges, van der Waals parameters, and torsional potentials propagate into free energy estimates [4, 33]. Machine learning force fields trained on quantum mechanical data offer the potential for improved accuracy, as demonstrated for hydration free energy calculations [4]. However, these force fields require careful validation for protein-ligand binding applications [4, 33].

Sampling Limitations

Incomplete sampling of protein and ligand conformational space remains a major source of error in FEP calculations [28, 21, 29]. Slow conformational changes, such as side chain rearrangements or loop motions, may not be adequately sampled within typical simulation timescales [28, 21]. Enhanced sampling methods and adaptive resource allocation strategies can mitigate these issues but do not eliminate them entirely [28, 30, 29].

System-Specific Considerations

The performance of FEP calculations varies significantly across different protein targets and ligand series [15, 17]. Factors such as binding site flexibility, water-mediated interactions, and the presence of metal ions or cofactors can affect accuracy [15, 17]. Benchmarking against experimental data for the specific system of interest is essential to establish confidence in FEP predictions [15, 17].

Future Directions

Ongoing developments in FEP methodology focus on improving accuracy, reducing computational cost, and expanding applicability to challenging targets [2, 29]. Advances in alchemical intermediate methods, such as PairMap and IMERGE-FEP, are enabling accurate calculations for chemically distant transformations [11, 24]. The integration of FEP with machine learning and automated workflows is making the technology more accessible to non-specialist users [12, 31, 32]. For veterinary drug discovery, FEP holds promise for designing species-selective therapeutics that target pathogens while minimizing off-target effects in host animals [1, 10].

References

[1] Chikhale RV, Islam MA, Suryawanshi VS et al. Free energy perturbation-derived identification of natural compounds targeting Mycobacterium tuberculosis thymidylate kinase. Sci Rep. 2026. https://pubmed.ncbi.nlm.nih.gov/42192177/

[2] Qian R, Xue J, Xu Y et al. Alchemical Transformations and Beyond: Recent Advances and Real-World Applications of Free Energy Calculations in Drug Discovery. J Chem Inf Model. 2024. https://pubmed.ncbi.nlm.nih.gov/39360948/

[3] Ansari N, Aviat F, Hénin J et al. Dual-LAO for calculating fast and robust relative binding free energies of simple and complex transformations. Commun Chem. 2026. https://pubmed.ncbi.nlm.nih.gov/42098254/

[4] Xie X, Weber JL, Svensson M et al. Accurate Hydration Free Energy Calculations for Diverse Organic Molecules With a Machine Learning Force Field. J Chem Theory Comput. 2026. https://pubmed.ncbi.nlm.nih.gov/41851056/

[5] Li Y, Liu S, Hu S et al. Structure-Based Design and Free Energy Perturbation-Guided Synthesis of Novel VEGFR2 and HDAC Dual Inhibitors with Potent Antitumor Efficacy. J Med Chem. 2026. https://pubmed.ncbi.nlm.nih.gov/41704127/

[6] Nozal V, Beyens O, Peeters S et al. Highly potent dipeptidyl peptidase 8/9 (DPP8/9) inhibitors designed via relative binding free energy calculations. Eur J Med Chem. 2025. https://pubmed.ncbi.nlm.nih.gov/40609224/

[7] Liesen MP, Hayes RL, Brooks Iii CL et al. Multiple Molecule λ-Dynamics: Probing Drug Resistance with Concurrent Protein and Ligand Perturbations. J Phys Chem Lett. 2025. https://pubmed.ncbi.nlm.nih.gov/40506019/

[8] Li Z, Jiang MY, Liu R et al. Discovery of highly potent phosphodiesterase-1 inhibitors by a combined-structure free energy perturbation approach. Acta Pharm Sin B. 2024. https://pubmed.ncbi.nlm.nih.gov/39807312/

[9] Harrison BA, Dowling JE, Bursavich MG et al. The Discovery of MORF-627, a Highly Selective Conformationally-Biased Zwitterionic Integrin αvβ6 Inhibitor for Fibrosis. J Med Chem. 2024. https://pubmed.ncbi.nlm.nih.gov/39446989/

[10] Islam SM, Rahman MS, Singla S et al. Computational Discovery of MERS-CoV Main Protease Inhibitors Through Screening and Molecular Dynamics Simulations. Proteins. 2026. https://pubmed.ncbi.nlm.nih.gov/41887806/

[11] Schoenmaker L, Jiskoot DA, Scheen J et al. IMERGE-FEP: Improving Relative Free Energy Calculation Convergence with Chemical Intermediates. J Phys Chem B. 2025. https://pubmed.ncbi.nlm.nih.gov/39976528/

[12] Liu R, Zhong Y, Yao Y et al. ALCHEMD: Bridging Accessibility and Accuracy in Automated Relative Binding Free Energy Workflows. J Chem Theory Comput. 2026. https://pubmed.ncbi.nlm.nih.gov/41474838/

[13] Zhou M, Shao X, Cai W et al. Zooming across the Alchemical Space. J Phys Chem Lett. 2025. https://pubmed.ncbi.nlm.nih.gov/40273944/

[14] Alencar Araripe D, Díaz-Holguín A, Poso A et al. Doing More with Less: Accurate and Scalable Ligand Free Energy Calculations by Focusing on the Binding Site. J Chem Inf Model. 2026. https://pubmed.ncbi.nlm.nih.gov/41685578/

[15] Sampangin Venkatesh SK, Das A, Lakkaniga NR. Evaluating the FEP+ protocol for predicting binding affinity of congeneric ligands towards various soluble proteins. RSC Med Chem. 2025. https://pubmed.ncbi.nlm.nih.gov/41355862/

[16] Gartan P, Brooks CL 3rd, Reuter N. Core Flipping in Lead Optimization: Rank Ordering Using λ-Dynamics. J Chem Inf Model. 2025. https://pubmed.ncbi.nlm.nih.gov/40522341/

[17] Ohadi D, Kumar K, Ravula S et al. Input Pose is Key to Performance of Free Energy Perturbation: Benchmarking with Monoacylglycerol Lipase. J Chem Inf Model. 2024. https://pubmed.ncbi.nlm.nih.gov/39560439/

[18] Zhang H, Im W. Ligand-Binding Free Energy Simulations of Membrane Proteins. Adv Exp Med Biol. 2026. https://pubmed.ncbi.nlm.nih.gov/41652165/

[19] Gupta SRR, Rameshwari R, Singh IK. Discovery of a potent ROR1 inhibitor using μs-scale MD simulations, wt-metadynamics, and absolute binding free energy calculations. J Mol Model. 2025. https://pubmed.ncbi.nlm.nih.gov/40996455/

[20] Dudas B, Athanasiou C, Mobarec JC et al. Quantifying Cooperativity through Binding Free Energies in Molecular Glue Degraders. J Chem Theory Comput. 2025. https://pubmed.ncbi.nlm.nih.gov/40326883/

[21] Liu R, Yao Y, Huang W et al. Divide-and-Conquer ABFE: Improving Free Energy Calculations by Enhancing Water Sampling. J Chem Theory Comput. 2025. https://pubmed.ncbi.nlm.nih.gov/40127297/

[22] Huang W, Liu R, Yao Y et al. RED-E-Function-Based Equilibrium Parameter Finder: Finding the Best Restraint Parameters in Absolute Binding Free Energy Calculations. J Phys Chem Lett. 2025. https://pubmed.ncbi.nlm.nih.gov/39718976/

[23] Liesen MP, Vilseck JZ. Superimposing Ligands with a Ligand Overlay as an Alternate Topology Model for λ-Dynamics-Based Calculations. J Phys Chem B. 2024. https://pubmed.ncbi.nlm.nih.gov/39515788/

[24] Furui K, Shimizu T, Akiyama Y et al. PairMap: An Intermediate Insertion Approach for Improving the Accuracy of Relative Free Energy Perturbation Calculations for Distant Compound Transformations. J Chem Inf Model. 2025. https://pubmed.ncbi.nlm.nih.gov/39800967/

[25] Liu R, Lai Y, Yao Y et al. State Function-Based Correction: A Simple and Efficient Free-Energy Correction Algorithm for Large-Scale Relative Binding Free-Energy Calculations. J Phys Chem Lett. 2025. https://pubmed.ncbi.nlm.nih.gov/40458911/

[26] Furui K, Ohue M. Benchmarking HelixFold3-Predicted Holo Structures for Relative Free Energy Perturbation Calculations. ACS Omega. 2025. https://pubmed.ncbi.nlm.nih.gov/40160764/

[27] Chen Y, Yang J. Acceleration of the GROMACS Free-Energy Perturbation Calculations on GPUs. ACS Omega. 2025. https://pubmed.ncbi.nlm.nih.gov/40521454/

[28] Koby SB, Gutkin E, Patel S et al. Automated On-the-Fly Optimization of Resource Allocation for Efficient Free Energy Simulations. J Chem Inf Model. 2025. https://pubmed.ncbi.nlm.nih.gov/40328725/

[29] Yao Y, Liu R, Li W et al. Convergence-Adaptive Roundtrip Method Enables Rapid and Accurate FEP Calculations. J Chem Theory Comput. 2024. https://pubmed.ncbi.nlm.nih.gov/39236257/ *** Disclaimer: This article is for educational and informational purposes only. It is not intended to substitute for professional veterinary advice, diagnosis, treatment, or regulatory guidance. Always consult a licensed veterinarian or qualified specialist regarding animal health, disease diagnosis, and therapeutic decisions.

[30] Xia Y, Lin X, Hu J et al. SPONGE-FEP: An Automated Relative Binding Free Energy Calculation Accelerated by Selective Integrated Tempering Sampling. J Chem Theory Comput. 2025. https://pubmed.ncbi.nlm.nih.gov/39868875/

[31] van Pinxteren DJM, Jespers W. Integrating Machine Learning into Free Energy Perturbation Workflows. J Chem Inf Model. 2025. https://pubmed.ncbi.nlm.nih.gov/40958764/

[32] Valsson Í, Warren MT, Deane CM et al. Narrowing the gap between machine learning scoring functions and free energy perturbation using augmented data. Commun Chem. 2025. https://pubmed.ncbi.nlm.nih.gov/39922899/

[33] Rufa D, Fass J, Chodera JD. Fine-tuning molecular mechanics force fields to experimental free energy measurements. bioRxiv. 2025. https://pubmed.ncbi.nlm.nih.gov/39829785/

[34] Kongtaworn N, Sinsulpsiri S, Hanpaibool C et al. Molecular Dynamics and Solvated Interaction Energy Prioritize Cannabidiol and Cannabinol as Variant-Spanning SARS-CoV-2 RBD-ACE2 Interface Blockers. Molecules. 2026. https://pubmed.ncbi.nlm.nih.gov/42075932/

[35] Friesner RA, Murphy RB, Zhang Y et al. Glide WS: Methodology and Initial Assessment of Performance for Docking Accuracy and Virtual Screening. J Chem Theory Comput. 2025. https://pubmed.ncbi.nlm.nih.gov/41339251/