Molecular Dynamics Simulations for Binding Affinity Refinement: From Fundamental Principles to Advanced Applications in Drug Discovery

Abstract

This article provides a comprehensive overview of molecular dynamics (MD) simulations for predicting and refining binding affinities in biomolecular complexes. It covers foundational principles, including force fields and thermodynamic ensembles, and details key methodological approaches such as alchemical free energy calculations (FEP, TI, BAR) and end-point methods (MM/PBSA, MM/GBSA). The content addresses critical challenges like sampling limitations and computational cost, offering optimization strategies and advanced sampling techniques. Furthermore, it explores the validation of MD-derived affinities against experimental data and presents comparative analyses of different methods. Designed for researchers, scientists, and drug development professionals, this review serves as a strategic guide for leveraging MD simulations to accelerate rational drug design and improve the accuracy of binding affinity predictions.

The Foundation of Binding Affinity: Core Principles and the Role of MD Simulations

Binding affinity is a fundamental parameter in drug discovery, quantifying the strength of interaction between a molecule and its target protein. Accurately predicting and measuring this interaction is crucial for the rapid development and optimization of novel therapeutics. The binding affinity is determined by the mechanism of recognition between proteins and ligands, with models including lock and key, induced fit, and conformational selection helping to explain this process. From a computational perspective, accurately predicting binding affinity remains a significant challenge, as current strategies often produce results that diverge by orders of magnitude from experimental data. This application note examines the key thermodynamic concepts of Kd, Ki, and IC50, explores their interrelationships, details experimental and computational protocols for their determination, and frames these methodologies within the context of molecular dynamics simulations for binding affinity refinement.

Key Parameters: Definitions and Theoretical Foundations

Thermodynamic Basis of Binding Affinity

The formation of a ligand-protein complex is a two-state process governed by the rates of association (on rate) and dissociation (off rate). For the interaction where L represents the ligand, P the protein, and LP the ligand-protein complex:

[ \text{L} + \text{P} \underset{\text{k}{\text{off}}}{\overset{\text{k}{\text{on}}}{\rightleftharpoons}} \text{LP} ]

The binding affinity is mathematically defined by the dissociation constant Kd, which represents the concentration of ligand required to occupy 50% of the available protein targets at equilibrium. The relationship between the rate constants and Kd is expressed as:

[ \text{K}d = \frac{\text{k}{\text{off}}}{\text{k}_{\text{on}}} = \frac{[\text{L}][\text{P}]}{[\text{LP}]} ]

where kon is the association rate constant (M⁻¹s⁻¹), koff is the dissociation rate constant (s⁻¹), and [L], [P], and [LP] represent molar concentrations. The affinity constant Ka is the reciprocal of Kd. The inhibition constant Ki represents the dissociation constant specifically for enzyme-inhibitor complexes measured through inhibition kinetics, whereas Kd is preferred when determined more directly through methods like isothermal titration calorimetry or surface plasmon resonance.

Table 1: Fundamental Parameters in Binding Affinity Characterization

Parameter	Definition	Units	Key Interpretation
Kd	Dissociation constant; concentration at which half binding sites are occupied	Molar (M)	Lower Kd indicates tighter binding
Ki	Inhibition constant; dissociation constant for enzyme-inhibitor binding	Molar (M)	Measured through inhibition kinetics
IC50	Half-maximal inhibitory concentration	Molar (M)	Functional measure of inhibition potency
kon	Association rate constant	M⁻¹s⁻¹	Speed of complex formation
koff	Dissociation rate constant	s⁻¹	Stability of the complex; related to residence time (1/koff)

Distinguishing Kd, Ki, and IC50

While often discussed interchangeably, Kd, Ki, and IC50 represent distinct parameters with important contextual differences. Kd provides a direct measure of binding affinity between two molecules, reflecting how tightly a compound binds to its target. It is a thermodynamic parameter intrinsic to the compound-target interaction. In contrast, IC50 is an operational measure of functional potency that quantifies the concentration necessary to reduce a specific biological process by half its maximum value under specific experimental conditions. While Kd offers insights into binding strength, IC50 reflects the composite effect of various interactions present in the assay system. The relationship between these parameters is critical: two compounds with the same Kd can have very different IC50 values in different experiments due to variations in experimental conditions such as target concentration and substrate levels.

Diagram 1: Parameter relationships in binding studies.

Experimental Methodologies and Protocols

Direct Binding Affinity Determination

Surface Plasmon Resonance (SPR) Protocol for Kd Determination: SPR enables real-time monitoring of molecular interactions without labeling requirements. The methodology involves immobilizing one binding partner (ligand) on a sensor chip surface and flowing the other partner (analyte) over it in solution.

• Sensor Chip Preparation: Select appropriate sensor chip (CM5 for covalent immobilization). Activate carboxyl groups with EDC/NHS mixture. Dilute ligand to 1-10 µg/mL in appropriate immobilization buffer. Inject ligand solution until desired immobilization level is reached. Deactivate excess reactive groups with ethanolamine.
• Equilibration: Condition sensor chip with running buffer until stable baseline achieved.
• Kinetic Titration: Prepare analyte serial dilutions in running buffer. Inject analyte concentrations in random order with adequate dissociation time between injections. Include blank runs for double-referencing.
• Data Analysis: Subtract reference cell and blank injection responses. Fit sensoryrams globally to 1:1 binding model to extract kon and koff. Calculate Kd from ratio koff/kon.

Isothermal Titration Calorimetry (ITC) Protocol: ITC directly measures heat changes during binding interactions, providing Kd, stoichiometry (n), and thermodynamic parameters (ΔH, ΔS).

• Sample Preparation: Thoroughly dialyze both protein and ligand against identical buffer. Precisely degas all solutions to prevent bubble formation.
• Instrument Setup: Load protein solution into sample cell. Fill syringe with ligand solution. Set appropriate temperature, reference power, and stirring speed.
• Titration Experiment: Program automated injections with sufficient spacing between injections. Typically use 15-25 injections of 2-10 µL each.
• Data Analysis: Integrate heat peaks per injection. Subtract dilution heats. Fit binding isotherm to appropriate model to extract Kd, n, and ΔH.

Table 2: Comparison of Key Experimental Methods for Binding Affinity Determination

Method	Key Measurements	Sample Consumption	Throughput	Key Advantages
SPR	kon, koff, Kd	Low (µg)	Medium	Real-time kinetics, label-free
ITC	Kd, n, ΔH, ΔS	High (mg)	Low	Direct thermodynamic profile
Fluorescence Spectroscopy	Kd, Ki, IC50	Low	Medium-High	Sensitive, adaptable to HTS
Radioligand Binding	Kd, Ki, Bmax	Very Low	Medium	Highly sensitive, established protocols
NanoBRET Target Engagement	Cellular Kd-apparent	Medium	High	Live-cell binding assessment

Competition Binding and IC50 Determination

Radioligand Binding Assay Protocol for Ki Determination: This method determines the inhibitory constant (Ki) of unlabeled compounds by measuring their ability to compete with a radiolabeled ligand for the target.

• Membrane Preparation: Homogenize tissue or cells in ice-cold buffer. Centrifuge at high speed to pellet membranes. Wash pellets and resuspend in assay buffer.
• Saturation Binding: Incubate varying concentrations of radioligand with fixed membrane protein concentration. Determine specific binding by subtracting nonspecific binding measured in presence of excess unlabeled ligand.
• Competition Binding: Incubate fixed concentration of radioligand (approximately Kd concentration) with varying concentrations of test compound. Use 8-12 concentrations in triplicate.
• Filtration and Detection: Terminate incubations by rapid filtration through GF/B filters. Wash filters and measure bound radioactivity by scintillation counting.
• Data Analysis: Fit competition curve to determine IC50. Convert IC50 to Ki using Cheng-Prusoff equation: Ki = IC50 / (1 + [L]/Kd), where [L] is radioligand concentration and Kd is its dissociation constant.

Conversion of IC50 to Kd Using the Cheng-Prusoff Equation: The Cheng-Prusoff equation provides a mathematical relationship to convert IC50 values to Ki under specific conditions:

[ Ki = \frac{IC{50}}{1 + \frac{[L]}{K_d}} ]

where [L] is the concentration of the tracer ligand and Kd is its dissociation constant. This conversion requires several assumptions: the binding must be at equilibrium, follow the law of mass action, involve competitive inhibition, and the radioligand concentration should not significantly deplete the free radioligand concentration.

The Scientist's Toolkit: Essential Research Reagents and Materials

Table 3: Key Research Reagent Solutions for Binding Affinity Studies

Reagent/Material	Function	Application Examples
Biacore Sensor Chips	Immobilization surface for binding partners	SPR-based kinetic studies
Radiolabeled Ligands	High-sensitivity tracer molecules	Competition binding assays
NanoBRET Tracers	Cell-permeable fluorescent probes	Live-cell target engagement assays
ITC Reference Buffer	Matched buffer for baseline correction	Accurate thermodynamic measurements
Protease Inhibitor Cocktails	Prevent protein degradation during assays	Maintain target integrity in long incubations
GF/B Filter Plates	Separate bound from free ligand	Filtration-based binding assays
Scintillation Cocktails	Detect bound radioligand	Radioligand binding quantification
Expression Systems	Recombinant protein production	Ensure adequate supply of purified target

Computational Refinement Through Molecular Dynamics Simulations

Molecular Dynamics for Binding Affinity Prediction

Molecular dynamics (MD) simulations provide a powerful approach for investigating protein-ligand binding mechanisms and predicting binding affinities. Unlike static crystal structures, MD captures the dynamic nature of proteins and ligands, which constantly interact with each other and their aqueous environment, following specific conformation distributions. A thermodynamic ensemble—obtained through MD simulation—represents the distribution of complex conformations in equilibrium, with all conformations collectively contributing to the affinity of protein-ligand binding.

Recent advances include Dynaformer, a graph-based deep learning model that leverages MD trajectories to predict binding affinities by learning geometric characteristics of protein-ligand interactions. This approach has demonstrated state-of-the-art scoring and ranking power, outperforming methods that rely solely on static crystal structures. The model was trained on a large-scale MD trajectory dataset containing 3,218 different protein-ligand complexes, with 10-nanosecond simulations performed for each complex and 100 snapshots sampled from each simulation to characterize the conformational space.

Diagram 2: MD workflow for affinity prediction.

Advanced Sampling Methods for Free Energy Calculations

Various computational methods have been developed to estimate binding free energies, each with different levels of rigor and computational cost:

Umbrella Sampling (US) with Restraints Protocol: This enhanced sampling method uses biased MD simulations along carefully chosen reaction coordinates.

• System Setup: Obtain initial protein-ligand structure. Solvate in explicit water box with appropriate ions. Energy minimize and equilibrate with position restraints on protein-ligand complex.
• Reaction Coordinate Definition: Select distance between protein and ligand centers of mass as primary reaction coordinate. Define additional restraints on orientation and conformation as needed.
• Steered MD: Pull ligand away from binding site along reaction coordinate to generate initial configurations for US windows.
• Umbrella Sampling: Run multiple parallel simulations (windows) with harmonic biases placed at different values along reaction coordinate. Ensure sufficient overlap between windows.
• PMF Construction: Use Weighted Histogram Analysis Method to combine data from all windows and obtain potential of mean force. Calculate binding affinity from PMF depth.

Binding Free Energy Estimator 2 (BFEE2) Methodology: BFEE2 addresses the substantial shift in configurational enthalpy and entropy following ligand-protein binding, which is challenging to represent in brute-force simulations. The software uses a stratification scheme with geometrical restraints to reduce conformational entropy of the biomolecular system, resulting in improved convergence of the potential of mean force calculation.

Table 4: Computational Methods for Binding Affinity Prediction

Method	Theoretical Basis	Accuracy (RMSE)	Computational Cost	Best Use Cases
Molecular Docking	Empirical/Knowledge-based scoring	2-4 kcal/mol	Low (CPU minutes)	Initial virtual screening
MM/PBSA/GBSA	Molecular mechanics with implicit solvent	1.5-3 kcal/mol	Low-Medium	Post-docking refinement
Umbrella Sampling	Potential of mean force along RC	1-2 kcal/mol	Medium-High	Detailed binding mechanism
Free Energy Perturbation	Alchemical transformation	0.5-1.5 kcal/mol	High	Lead optimization
Dynaformer (ML-MD)	Graph neural networks on MD trajectories	~1 kcal/mol	High	High-accuracy ranking

Case Study: GPCR Binding Affinity Prediction

G protein-coupled receptors represent an important class of drug targets where MD simulations have provided valuable insights. Recent studies have performed binding free energy calculations for β1 adrenergic receptor (β1AR) in both active and inactive states using a modified Bennett Acceptance Ratio method adapted for membrane proteins. The results showed significant correlation with experimental pKD values, demonstrating the potential of computational approaches to capture state-dependent binding affinities. Spatial distribution function analysis revealed that agonists bound to inactive states of β1AR showed higher mobility than those bound to active states, with orthosteric binding pockets in active states being more compact, resulting in tighter ligand-receptor contacts.

The accurate determination and prediction of binding affinity remains a cornerstone of drug discovery, with Kd, Ki, and IC50 providing complementary information about compound-target interactions. While experimental methods directly measure these parameters, computational approaches—particularly molecular dynamics simulations—offer powerful tools for understanding binding mechanisms and predicting affinities. The integration of machine learning with MD trajectories, as exemplified by Dynaformer, represents a promising avenue for improving prediction accuracy. As these methods continue to evolve, they will play an increasingly important role in accelerating the drug discovery process by enabling more reliable in silico screening and optimization of candidate compounds.

Molecular dynamics (MD) is a powerful computational technique that simulates the physical movements of atoms and molecules over time. By applying the principles of classical mechanics, MD provides insights into the time-dependent behavior of molecular systems, offering a dynamic view that is often inaccessible through experimental methods alone. Within drug discovery and development, MD simulations have become an indispensable tool for binding affinity refinement, allowing researchers to probe the interactions between proteins and ligands at an atomic level. This protocol outlines the theoretical foundations of MD, with a specific focus on its application in predicting and optimizing the binding affinity of small molecules to their target receptors, a critical step in rational drug design.

Theoretical Foundations

Newton's Laws of Motion: The Governing Principles

The entire framework of molecular dynamics is built upon Newton's laws of motion. MD simulations calculate the time evolution of a molecular system by solving Newton's equations of motion for each atom. The core equation is derived from Newton's second law, which states that the force F acting on an object is equal to its mass m multiplied by its acceleration a:

F = m × a [1]

In a molecular system, the force on an atom can also be expressed as the negative gradient of the potential energy V with respect to the atom's position r:

F = -∇V(r) [1]

This fundamental relationship connects the potential energy of a system, which arises from the interactions between atoms, to the forces that drive their motion. By integrating these equations over time, MD simulations generate a trajectory that describes how the positions and velocities of all atoms in the system change over time.

Force Fields: The Mathematical Heart of MD

The concept of a force field is central to molecular dynamics. A force field is a mathematical model that describes the potential energy of a system as a function of the nuclear coordinates. The accuracy of an MD simulation is critically dependent on the quality of the force field used. Most modern force fields express the total potential energy as a sum of several terms that account for different types of intra- and intermolecular interactions:

V_total = V_bonded + V_non-bonded = V_bonds + V_angles + V_dihedrals + V_{electrostatic} + V_{van der Waals} [1]

These force fields, such as CHARMM, AMBER, and OPLS, are parameterized for specific types of molecules and interactions [2] [3]. For example, CHARMM is widely used for biological systems, while AMBER excels in modeling nucleic acids and proteins. The choice of force field can dramatically influence the outcome of a simulation, making careful selection essential for meaningful results, particularly in binding affinity studies where the accurate representation of intermolecular forces is paramount.

Integration Algorithms: Propagating the System Through Time

To generate a dynamical trajectory, the equations of motion must be solved numerically. This is achieved through integration algorithms that calculate atomic positions and velocities at a future time t + Δt based on the current state at time t. The time step Δt is a critical parameter, typically on the order of femtoseconds (10^-15 seconds), to ensure numerical stability and accurately capture the fastest vibrations in the system (e.g., bond stretching).

Several efficient numerical algorithms are commonly used, including the Verlet method, the leap-frog method, and the velocity Verlet method [1]. These algorithms strike a balance between computational efficiency, numerical accuracy, and stability, allowing for the simulation of systems over nanoseconds to microseconds, and sometimes even beyond, which is necessary to observe biologically relevant processes like ligand binding and unbinding.

Table 1: Key Components of a Classical Force Field for Biomolecular Simulation

Energy Component	Mathematical Form	Physical Description
Bond Stretching	V = kb(r - r0)2	Energy required to stretch or compress a bond from its equilibrium length.
Angle Bending	V = kθ(θ - θ0)2	Energy required to bend an angle from its equilibrium value.
Torsional Dihedral	V = kφ[1 + cos(nφ - δ)]	Energy associated with rotation around a central bond.
van der Waals	V = 4ε[(σ/r)12 - (σ/r)6]	Non-bonded interactions due to transient dipoles (Lennard-Jones potential).
Electrostatic	V = (q1q2)/(4πε0r)	Non-bonded interactions between permanent partial charges.

Table 2: Common Numerical Integrators in Molecular Dynamics

Algorithm	Key Features	Advantages	Common Time Step (fs)
Verlet	Calculates new positions r(t+Δt) using positions and forces from previous step.	Time-reversible, good energy conservation.	1-2
Leap-Frog	Velocities and positions are calculated at staggered times.	Computationally efficient.	1-2
Velocity Verlet	Positions, velocities, and accelerations are computed at the same time.	Numerically stable, calculates velocities explicitly.	1-2

Application Note: Protocol for Binding Affinity Refinement

The following protocol describes a methodology for using MD simulations to refine binding affinity predictions for protein-ligand complexes, with a specific application to G-Protein Coupled Receptors (GPCRs). The protocol leverages the Bennett Acceptance Ratio (BAR) method for binding free energy calculation, which has demonstrated a strong correlation (R² = 0.7893) with experimental binding data [2].

Detailed Protocol: BAR-Based Binding Free Energy Calculation for GPCRs

This protocol is adapted from a study that successfully predicted binding affinities for agonists bound to the β1 adrenergic receptor (β1AR) [2].

Stage 1: Initial System Setup

• Structure Preparation: Obtain the initial coordinates for the GPCR-ligand complex from a crystal structure or a high-quality homology model. For the β1AR case, structures were sourced from the Protein Data Bank (PDB).
• Force Field Selection: Parameterize the system using an appropriate all-atom force field (e.g., CHARMM or AMBER). Assign partial charges and other parameters to the ligand using tools compatible with the chosen force field.
• System Solvation and Ionization:
  • Embed the protein-ligand complex in a pre-equilibrated membrane patch (e.g., POPC lipid bilayer) to mimic the native GPCR environment.
  • Solvate the entire system, including membrane, protein, and ligand, in a water model such as TIP3P.
  • Add ions (e.g., Na⁺ and Cl⁻) to neutralize the system's net charge and to achieve a physiologically relevant salt concentration (e.g., 0.15 M NaCl).

Stage 2: Simulation Equilibration

• Energy Minimization: Perform steepest descent or conjugate gradient minimization to remove any bad steric contacts introduced during system setup. Run until the maximum force is below a chosen threshold (e.g., 1000 kJ/mol/nm).
• Position-Restrained MD (NVT Ensemble): Run a short simulation (e.g., 100 ps) with heavy atoms of the protein and ligand harmonically restrained. This allows the solvent and ions to relax around the fixed solute. Use a temperature coupling algorithm (e.g., Berendsen or Nosé-Hoover) to maintain the system at the target temperature (e.g., 310 K).
• Position-Restrained MD (NPT Ensemble): Run a second short simulation (e.g., 100 ps) with the same restraints, but now use a pressure coupling algorithm (e.g., Berendsen or Parrinello-Rahman) to maintain the system at the target pressure (e.g., 1 bar). This ensures the correct density of the system.
• Unrestrained MD (NPT Ensemble): Run a final equilibration phase (e.g., 1-5 ns) without any positional restraints to allow the entire system to fully equilibrate. Monitor properties like potential energy, temperature, pressure, and root-mean-square deviation (RMSD) of the protein backbone to confirm stability.

Stage 3: Production Simulation and Free Energy Calculation

• Alchemical Setup: Define the perturbation pathway that "annihilates" the ligand in the bound state (protein-ligand complex) and in the unbound state (ligand in solution). Divide this pathway into a series of intermediate states, defined by a scaling parameter λ (where λ=0 is the fully interacting state and λ=1 is the non-interacting state). Typically, 10-20 λ windows are used for sufficient sampling.
• Sampling at Lambda Windows: For each λ window, run an independent MD simulation to sample the configuration space. The simulation length per window should be sufficient for convergence; studies have used anywhere from 1-10 ns per window. This can be run using the GROMACS simulation package [2].
• Binding Free Energy Calculation with BAR: Use the Bennett Acceptance Ratio (BAR) method to calculate the free energy difference between adjacent λ windows. The BAR method optimally uses the work values from both the forward (λ→λ+Δλ) and backward (λ+Δλ→λ) transitions to provide a statistically robust estimate of the free energy change. The total binding free energy (ΔG_bind) is the sum of these increments over the entire alchemical path.

The Scientist's Toolkit: Essential Research Reagents and Software

Table 3: Key Computational Tools for MD-Based Binding Affinity Studies

Tool / Resource	Type	Function in Research	Application Context
GROMACS [2]	MD Simulation Engine	High-performance software package used to run energy minimization, equilibration, and production MD simulations.	Core simulation engine for the BAR protocol on GPCR targets.
CHARMM/AMBER [2] [3]	Force Field	Provides the mathematical parameters for bonded and non-bonded interactions, defining the physics of the simulation.	Essential for accurate parameterization of proteins, membranes, and ligands.
OpenMM Toolkit [4]	MD Simulation Library	A versatile, high-performance toolkit for molecular simulation. It forms the core of automated pipelines like drMD.	Used in drMD for running publication-quality simulations with reduced expertise.
BAR Method [2]	Free Energy Algorithm	An alchemical method for calculating binding free energies from simulations at multiple lambda states.	Core analysis method for achieving high correlation with experimental binding data.
OMol25 Dataset [5] [6]	Training Dataset	A massive dataset of over 100 million DFT-level molecular snapshots used to train machine-learning force fields.	Foundation for next-generation, highly accurate neural network potentials (NNPs).
Neural Network Potentials (NNPs) [3]	Machine-Learning Force Field	Force fields that learn interatomic potentials from quantum mechanical data, offering DFT accuracy at lower cost.	Used to overcome limitations of classical force fields, e.g., for charged fluids and reactive systems.
drMD [4]	Automated Pipeline	Simplifies running MD simulations using OpenMM, reducing the barrier to entry for non-experts.	Provides an automated, user-friendly interface for running simulations.

Emerging Trends and Future Perspectives

The field of molecular dynamics is rapidly evolving with the integration of machine learning (ML) and the development of novel, highly accurate datasets. The recent release of the Open Molecules 2025 (OMol25) dataset, containing over 100 million molecular snapshots calculated with density functional theory (DFT), is poised to revolutionize the field [5] [6]. This dataset enables the training of neural network potentials (NNPs) that can predict molecular energies and forces with near-DFT accuracy but at a fraction of the computational cost, unlocking the simulation of large, chemically diverse systems that were previously out of reach [3].

Furthermore, foundation models like Boltz-2 are beginning to bridge the gap between structural prediction and binding affinity estimation. These AI models can approach the accuracy of rigorous methods like free energy perturbation (FEP) while being orders of magnitude faster, potentially transforming virtual screening and lead optimization workflows in drug discovery [7]. As these technologies mature, they will integrate with and enhance the classical theoretical basis of molecular dynamics, providing even more powerful tools for understanding and predicting molecular behavior.

Molecular dynamics (MD) simulation serves as a computational microscope, enabling researchers to observe the motion and interactions of atoms in biomolecular systems over time. The thermodynamic state of an MD simulation is defined by its statistical ensemble, which specifies the macroscopic variables held constant during the run. For researchers focusing on binding affinity refinement, selecting the appropriate ensemble is not merely a technical detail but a fundamental decision that directly impacts the biological relevance and thermodynamic accuracy of the simulation results.

This application note details the three essential ensembles—NVE (microcanonical), NVT (canonical), and NPT (isothermal-isobaric)—for biomolecular simulations. We provide theoretical foundations, practical implementation protocols, and specific applications to binding free energy calculations, with a special focus on membrane proteins like G-protein coupled receptors (GPCRs), crucial targets in modern drug discovery.

Theoretical Foundations of MD Ensembles

The Statistical Mechanical Basis

Thermodynamic ensembles are artificial constructs that facilitate the derivation of system properties through statistical mechanics. Although averages of basic thermodynamic properties become consistent across different ensembles in the thermodynamic limit, fluctuations—which are critical for calculating properties like specific heat or compressibility—vary significantly between ensembles [8]. This ensemble equivalence breaks down for finite systems, making the choice of ensemble critical for accurate biomolecular simulation.

Ensemble Selection for Biomolecular Systems

The selection of an ensemble should be guided by the experimental conditions being modeled and the free energy potential of interest. In principle, different ensembles sample different free energies: NVE corresponds to internal energy, NVT to Helmholtz free energy, and NPT to Gibbs free energy [9]. For biomolecular simulations in aqueous environments, which naturally occur at constant temperature and pressure, the NPT ensemble most closely mimics experimental conditions [9] [10]. A standard MD protocol typically proceeds through multiple ensembles, beginning with NVT for thermal equilibration, followed by NPT for density equilibration, and concluding with a production run in the NPT ensemble [10].

Essential Ensembles: Characteristics and Applications

Table 1: Key Characteristics of Primary MD Ensembles

Ensemble	Constant Variables	Controlled via	Primary Applications in Biomolecular Simulation
NVE (Microcanonical)	Number of particles (N), Volume (V), Energy (E)	Newton's equations without temperature/pressure control [8]	- Exploring constant-energy conformational space [8]- Spectral calculations from velocity correlation functions [9]- Testing conservation properties of integrators
NVT (Canonical)	Number of particles (N), Volume (V), Temperature (T)	Thermostats (e.g., Nosé-Hoover, Langevin, Berendsen) [11] [12]	- Conformational searches in vacuum [8]- Simulations where volume change is negligible [11]- Ion diffusion in solids, surface adsorption/reactions [11]
NPT (Isothermal-Isobaric)	Number of particles (N), Pressure (P), Temperature (T)	Thermostat + Barostat (e.g., Parrinello-Rahman, Berendsen) [13] [14]	- Mimicking standard laboratory conditions [9] [10]- Simulating thermal expansion, phase transitions [14]- Binding affinity calculations for drug discovery [2]

NVE (Microcanonical) Ensemble

The NVE ensemble models an isolated system that cannot exchange energy or matter with its surroundings. Energy conservation is its defining characteristic, though minor numerical drift occurs due to integration errors [8]. In the Verlet leapfrog algorithm, potential and kinetic energies are calculated half a step out of synchrony, contributing to energy fluctuations [8].

Biomolecular Application Considerations: NVE simulations are not recommended for system equilibration due to the lack of controlled energy flow needed to reach a desired temperature [8]. However, they are crucial for specific properties; for instance, infrared spectrum calculation from velocity autocorrelation functions requires an NVE production run because thermostats decorrelate velocities, which would destroy the signal [9].

NVT (Canonical) Ensemble

The NVT ensemble models a system at constant temperature, allowing heat exchange with an external thermostat. This is typically achieved by scaling atomic velocities or using more sophisticated algorithms like Nosé-Hoover or Langevin dynamics [11] [12].

Biomolecular Application Considerations: NVT is the default ensemble for simulations without periodic boundary conditions, such as conformational searches of molecules in vacuum, where volume, pressure, and density are undefined [8]. It is also appropriate for systems where volume changes are negligible, such as ion diffusion in solids or reactions on stable surface structures [11]. A significant limitation is that the average pressure depends on the initial simulation box volume, which may not match experimental conditions [12].

NPT (Isothermal-Isobaric) Ensemble

The NPT ensemble maintains constant temperature and pressure, requiring both a thermostat and a barostat. The barostat adjusts the simulation box volume to maintain constant pressure [13] [10]. Popular methods include the Parrinello-Rahman and Berendsen barostats [14].

Biomolecular Application Considerations: NPT is the ensemble of choice for most solution-phase biomolecular simulations as it correctly reproduces experimental densities and pressures [8] [14]. It is particularly critical for binding affinity calculations where the correct system density is essential for accurate free energy estimation [2]. Special care is needed for systems with limited long-range order to prevent irreversible cell deformation [13].

Implementation Protocols

Thermostat and Barostat Selection

Table 2: Common Temperature and Pressure Control Methods

Method	Type	Key Parameters	Advantages	Limitations
Nosé-Hoover	Thermostat	SMASS (virtual mass) [12]	- Deterministic, reproduces correct NVT ensemble [11]	- Can lack ergodicity in simple systems [9]
Langevin	Thermostat	LANGEVIN_GAMMA (friction coefficient) [11] [13]	- Good for mixed phases; controls atoms individually [11]	- Stochastic; trajectories not reproducible [11]
Berendsen	Thermostat/Baristat	taut (time constant), taup, compressibility [11] [14]	- Simple, efficient, good convergence [11] [14]	- Can produce unnatural phenomena, does not yield correct ensemble [11]
Parrinello-Rahman	Barostat	pfactor (barostat parameter) [14]	- High flexibility; controls all cell degrees of freedom [14]	- Requires careful parameter tuning (e.g., pfactor) [14]

Standard Protocol for Biomolecular Simulation

A typical MD procedure for binding affinity studies utilizes multiple ensembles in sequence to properly equilibrate the system [10]:

Diagram 1: Standard MD Protocol for Binding Affinity Studies

• Energy Minimization: Removes steric clashes and unfavorable contacts in the initial structure.
• NVT Equilibration: Heats the system to the target temperature while keeping volume fixed. Typical duration: 100-500 ps. Key parameter: temperature coupling constant (ttime ~ 20 fs in ASE [14]).
• NPT Equilibration: Adjusts system density to the target pressure. Typical duration: 1-5 ns. Key parameters: pressure coupling constant (pfactor ~ 10⁶-10⁷ GPa·fs² for solids [14]) and target pressure (1 bar for biomolecules).
• NPT Production Run: Generates trajectory for analysis. Duration depends on the biological process (ns to μs).

VASP Implementation Examples

NVT Ensemble with Nosé-Hoover Thermostat [12]:

NPT Ensemble with Parrinello-Rahman Barostat [13]:

Application in Binding Affinity Refinement

Case Study: GPCR Binding Affinity Prediction

A recent study demonstrated the successful application of NPT-MD simulations for binding free energy calculations on GPCR targets using a re-engineered Bennett Acceptance Ratio (BAR) method [2]. The protocol achieved significant correlation (R² = 0.7893) with experimental binding affinities (pK_D) for various agonists bound to β1 adrenergic receptor (β1AR) in both active and inactive states [2].

Key Methodology [2]:

• System Preparation: Built membrane-protein-ligand systems in explicit lipid bilayers and water.
• Equilibration: Applied multi-step NPT protocol to stabilize the membrane-solvent environment.
• Production: Conducted multiple independent NPT-MD simulations (≥ 50 ns) for each ligand-receptor complex.
• Free Energy Calculation: Utilized the BAR method over multiple intermediate states (λ windows) to compute binding free energies.

The Scientist's Toolkit: Essential Research Reagents

Table 3: Essential Computational Tools for Biomolecular MD

Tool Category	Specific Examples	Function in Biomolecular Simulation
Simulation Software	GROMACS [2], CHARMM [2], AMBER [2], VASP [15] [13] [12]	- Provides force fields, integrators, and analysis tools for running MD simulations
Thermostats	Nosé-Hoover [12], Langevin [11] [13], Berendsen [11], CSVR [12]	- Controls system temperature by scaling velocities or applying stochastic forces
Barostats	Parrinello-Rahman [13] [14], Berendsen [14]	- Maintains constant pressure by adjusting simulation box volume and shape
Force Fields	ASAP3-EMT [11] [16], PFP [14], AMBER/CHARMM lipid force fields	- Defines potential energy functions and parameters for interatomic interactions
Free Energy Methods	BAR (Bennett Acceptance Ratio) [2], FEP (Free Energy Perturbation) [2]	- Calculates binding free energies from alchemical perturbation simulations

The strategic selection and proper implementation of MD ensembles form the foundation of reliable biomolecular simulations, particularly in binding affinity refinement studies. The NVE ensemble serves specialized purposes, NVT is ideal for volume-constrained systems, while NPT most accurately mimics standard laboratory conditions for solution-phase biology. Through careful parameter selection and a structured equilibration protocol, researchers can generate thermodynamically valid ensembles for precise binding free energy calculations. The successful application of these principles to GPCR drug targets highlights their critical role in advancing computational structural biology and rational drug design.

Molecular dynamics (MD) simulations have emerged as a powerful tool in structural biology and drug discovery, capable of bridging the gap between static snapshots from crystallography and the dynamic reality of protein-ligand interactions. Whereas X-ray structures provide crucial atomic-resolution information, they represent merely a single conformational state, often missing the intricate pathways and transitional conformations that define molecular recognition processes. MD simulations address this limitation by modeling the time-dependent evolution of molecular systems, offering unique insights into binding mechanisms, conformational plasticity, and the thermodynamic and kinetic parameters governing molecular interactions.

The importance of capturing these dynamic processes is particularly critical for binding affinity refinement research, where understanding the complete binding pathway—not just the final bound state—can inform more effective drug design strategies. This application note details how MD simulations, especially when enhanced by advanced sampling and analysis techniques, can reveal the critical binding conformations and pathways that underlie biological function and therapeutic intervention.

Theoretical Framework: From Static Structures to Dynamic Ensembles

Proteins are inherently dynamic molecules that continuously fluctuate and transition between various conformational states while performing their biological functions. These dynamics occur across multiple timescales, from femtosecond bond vibrations to millisecond domain motions, and are essential for mechanisms like induced-fit binding and allosteric regulation [17]. MD simulations explicitly model these motions by numerically solving Newton's equations of motion for all atoms in the system, typically using a 2-femtosecond integration time step [18].

The fundamental MD algorithm involves three core steps repeated for thousands to millions of iterations:

• Force computation: Calculating forces on each atom based on the potential energy function (force field)
• Configuration update: Updating atomic positions and velocities by integrating equations of motion
• Data output: Periodically saving coordinates, velocities, and energies for analysis [18]

This process generates trajectories—sequences of snapshots representing the system's conformational evolution over time. Analysis of these trajectories reveals not only the most probable conformations but also the transitions between them, providing a comprehensive view of the energy landscape governing binding events.

Advanced MD Methodologies for Capturing Binding Events

Enhanced Sampling Techniques

Conventional MD simulations often struggle to capture rare events like ligand binding and unbinding within practical computational timeframes. Enhanced sampling methods address this limitation by accelerating the exploration of conformational space.

Hypersound-accelerated MD introduces high-frequency ultrasound perturbation (e.g., 625 GHz) to accelerate slow biomolecular processes. This method generates local high-temperature/pressure regions that promote barrier crossing without altering macroscopic system properties. In CDK2-inhibitor binding studies, hypersound acceleration increased binding event observation by 9.6-17.7 times compared to conventional MD, enabling detailed pathway analysis in 100-200 ns simulations [19].

Other Enhanced Sampling Approaches include replica-exchange MD (REMD) for improved conformational sampling [17] and Markov State Models (MSMs) that aggregate many short simulations to reconstruct long-timescale kinetics [19]. These methods collectively enhance the efficiency of capturing binding-related conformational transitions.

Machine Learning-Enhanced Analysis

The high-dimensional data produced by MD simulations presents significant analysis challenges. Dimensionality reduction techniques are essential for identifying meaningful patterns in these complex datasets.

Principal Component Analysis (PCA) operates on atomic coordinates from MD trajectories to identify dominant modes of collective motion. By projecting high-dimensional data onto orthogonal principal components that capture maximum variance, PCA reveals large-scale conformational changes relevant to binding, such as hinge bending and domain motions [17]. When applied to HIV-1 protease simulations, PCA revealed flap-closing motions critical for inhibitor binding, differentiating the mechanisms of different inhibitors [17].

Self-Organising Maps (SOMs) combined with hierarchical clustering provide a two-level approach for analyzing conformational ensembles. SOMs create a topological mapping of conformational space onto a 2D grid, where similar structures cluster in neighboring neurons. Subsequent clustering of prototype vectors identifies functionally relevant conformations and highlights differences in conformational dynamics across protein variants [20].

Table 1: Performance Comparison of MD Methodologies in Capturing Binding Events

Methodology	Simulation Time	Binding Probability	Key Advantages	Representative Application
Conventional MD	100 ns	0.5-0.7%	Baseline sampling	Simple protein-ligand systems
Hypersound-Accelerated MD	100-200 ns	4.8-12.4% (9.6-17.7x acceleration)	Direct observation of binding pathways	CDK2-inhibitor binding [19]
SOM with Hierarchical Clustering	40 ns per system	N/A	Topological mapping of conformational space; comparison of multiple variants	α-spectrin SH3 domain dynamics [20]
PCA combined with MD	Varies by system	N/A	Identification of collective motions; separation of timescales	HIV-1 protease flap dynamics [17]

Application Notes: Protocol for Binding Pathway Analysis

System Setup and Simulation Parameters

Initial Structure Preparation Begin with high-resolution crystal structures of the target protein, ensuring completeness of missing residues and loops. For ligand binding studies, obtain accurate force field parameters for small molecules using tools like CGenFF or GAFF. Place the protein in an appropriate simulation box (e.g., rectangular or dodecahedral) with a minimum 1.2 nm distance between the protein and box boundaries, solvate with explicit water models (e.g., SPC, TIP3P), and add ions to neutralize system charge [20] [18].

Simulation Parameters

• Force Field: Select appropriate force fields (e.g., GROMOS96 43a2, AMBER, CHARMM) based on system requirements [20] [21]
• Electrostatics: Employ Particle Mesh Ewald (PME) summation with 9 Å direct space cutoff [20]
• Thermostat: Use Berendsen or Nosé-Hoover thermostats with separate coupling for protein and solvent (0.1 ps coupling constant) [20]
• Barostat: Apply Parrinello-Rahman or Berendsen barostat for pressure coupling (1.0 ps coupling constant) [20]
• Constraints: Implement LINCS algorithm for bond constraints and SETTLE for water geometry [20]
• Integration: Utilize leap-frog integrator with 2-fs time step [18]

Enhanced Sampling for Binding Events

Hypersound-Accelerated MD Protocol

• Parameter Selection: Set hypersound frequency to 625 GHz (1.6 ps period) with amplitude parameter (N = 50, vmax = 400 m/s) for optimal acceleration [19]
• Simulation Setup: Run multiple independent simulations (100-200 ns) with ligands initially placed randomly in solution
• Trajectory Analysis: Identify binding events by monitoring ligand-protein distance and interaction stability
• Pathway Clustering: Group observed binding pathways based on structural similarity and energy profiles

Binding Pathway Analysis Workflow The following diagram illustrates the comprehensive workflow for capturing and analyzing binding pathways using enhanced MD simulations:

Conformational Clustering and Pathway Analysis

Self-Organising Maps for Conformational Analysis

• Input Data Preparation: Extract protein conformations from MD trajectories and project onto essential dynamics space using PCA [20]
• SOM Training: Optimize SOM parameters (map size, learning rate) using experimental design approaches like Taguchi methods [20]
• Prototype Vector Clustering: Perform hierarchical clustering (e.g., complete linkage) on SOM prototype vectors to identify major conformational states [20]
• Representative Structure Selection: Extract central structures from each cluster for detailed analysis of binding-relevant conformations

Principal Component Analysis Protocol

• Coordinate Preparation: Align all trajectory frames to a reference structure to remove global translation and rotation [17]
• Covariance Matrix Construction: Calculate the covariance matrix of atomic fluctuations (typically Cα atoms only) [17]
• Diagonalization: Solve eigenvalue problem to obtain principal components (PCs) and corresponding variances [17]
• Projection: Project trajectories onto dominant PCs to visualize conformational space and identify metastable states [17]
• Mode Analysis: Interpret physical significance of dominant PCs through visualization and cross-correlation analysis

Case Studies: MD in Binding Affinity Research

CDK2-Inhibitor Binding Pathways

Hypersound-accelerated MD simulations of CDK2 with slow-binding inhibitors CS3 and CS242 revealed conformationally and energetically diverse binding pathways, challenging the conventional single-pathway model. Analysis of 67 (CS3) and 14 (CS242) independent binding trajectories showed:

• Multiple Transition States: Each binding pathway featured distinct energy barriers at different structural locations
• Pathway-Specific Barriers: The highest-energy transition state varied between pathways, occurring either during pocket entry or internal conformational rearrangement
• Asymmetric Binding/Unbinding: Unbinding trajectories often followed different pathways than binding, contradicting simple kinetic models [19]

These insights enabled estimation of activation energies (3.9 ± 1.8 kcal/mol for CS3, 6.7 ± 2.4 kcal/mol for CS242) and Arrhenius parameters, connecting microscopic simulation data with macroscopic kinetic measurements [19].

HIV-1 Protease Inhibitor Binding Mechanism

PCA analysis of MD simulations revealed distinct flap-closing mechanisms in HIV-1 protease upon binding different inhibitors. The study demonstrated:

• Stronger Inhibitor Binding: 4UY formed more hydrogen bonds and hydrophobic interactions than DRV, resulting in 4.3-6.4 kcal/mol stronger binding affinity [17]
• Collective Motions: PCA identified specific flap-closing motions correlated with inhibitor potency [17]
• Residue-Specific Interactions: Analysis highlighted key residues differentiating binding mechanisms to wild-type and mutant variants [17]

α-Spectrin SH3 Domain Conformational Dynamics

The two-level SOM/hierarchical clustering approach applied to MD simulations of α-spectrin SH3 domain and six mutants successfully:

• Mapped Conformational Space: Created 2D topological visualization of complex conformational distributions [20]
• Identified Functional States: Highlighted differences in conformational dynamics directly related to biological function [20]
• Compared Multiple Systems: Enabled direct comparison of conformational ensembles across wild-type and mutant proteins [20]

Table 2: Key Research Reagents and Computational Tools for MD Binding Studies

Resource Category	Specific Tools/Reagents	Function/Application	Key Features
MD Simulation Software	GROMACS [20] [18], AMBER [21], DESMOND [21]	Biomolecular MD simulations	Optimized force fields; GPU acceleration; Enhanced sampling methods
Enhanced Sampling Methods	Hypersound-accelerated MD [19], Replica Exchange MD [17]	Accelerating rare events; Improved conformational sampling	External field perturbation; Parallel tempering
Analysis Algorithms	Principal Component Analysis [17], Self-Organising Maps [20]	Dimensionality reduction; Conformational clustering	Identification of collective motions; Topological mapping
Force Fields	GROMOS96 43a2 [20], AMBER, CHARMM	Potential energy functions	Atomistic parameters; Bonded and non-bonded terms
Visualization Software	VMD, PyMOL, Chimera	Trajectory visualization; Structure analysis	Molecular graphics; Animation; Plugin support

Data Analysis and Interpretation

Quantitative Analysis of Binding Trajectories

Kinetic Parameter Extraction From observed binding events in accelerated MD, calculate association rate constants (kon) using: kon = (Nbound / Ntotal) × (1 / [L] × t) where Nbound is trajectories with stable binding, Ntotal is total trajectories, [L] is ligand concentration, and t is simulation time [19]. For CDK2 inhibitors, this yielded kon = 3.68×10⁶ M⁻¹s⁻¹ for CS3 and 1.92×10⁶ M⁻¹s⁻¹ for CS242 [19].

Energy Landscape Construction Compute free energy landscapes by projecting trajectories onto collective variables (e.g., PC1 vs. PC2) and calculating potential of mean force: G(x) = -kBT ln P(x) where P(x) is probability distribution along coordinate x [17]. This identifies stable conformations and transition states along binding pathways.

Conformational Analysis Metrics

Essential Dynamics Convergence Evaluate PCA convergence using scree plots (eigenvalue vs. mode index) and ensure sufficient sampling of conformational space [17]. Well-converged simulations typically show rapid eigenvalue decay with first few PCs capturing majority of variance.

Cluster Validation For SOM/hierarchical clustering, assess quality using Davies-Bouldin index or silhouette coefficients to determine optimal cluster number [20]. Higher values indicate better separation between conformational states.

Integration with Binding Affinity Refinement

The true value of MD-derived binding conformations and pathways lies in their integration with binding affinity prediction methods. Studies demonstrate that incorporating MD-based structural dynamic information (4D descriptors) can improve machine learning predictions of binding affinity, particularly for targets with considerable structural flexibility like TAF1-BD2 [22]. However, the initial structure quality critically influences the final predictive performance, emphasizing the need for adequate conformational sampling [22].

The following diagram illustrates how MD-derived structural dynamics integrate with binding affinity prediction frameworks:

This integrated approach enables researchers to move beyond static structural snapshots and incorporate dynamic information that more accurately reflects the physical processes underlying molecular recognition and binding.

A Practical Guide to Binding Affinity Calculation Methods in MD Simulations

The accurate prediction of binding free energies is a central challenge in computational biophysics and structure-based drug design. Alchemical free energy calculations, which include Free Energy Perturbation (FEP), Thermodynamic Integration (TI), and the Bennett Acceptance Ratio (BAR) method, provide a rigorous, physics-based framework for addressing this challenge by leveraging molecular dynamics (MD) simulations and statistical mechanics [23]. These methods enable the calculation of relative free energy differences by simulating the alchemical transformation of one molecule into another along a non-physical pathway, thus avoiding the need to simulate the actual binding and unbinding processes, which occur on timescales often inaccessible to MD [24].

The integration of these methods into drug discovery pipelines has been transformative, allowing for the efficient optimization of lead compounds with accuracy often approaching experimental error (∼1 kcal/mol) [25] [24]. This application note details the protocols, applications, and recent advancements for FEP, TI, and BAR, providing a practical guide for researchers aiming to implement these techniques for binding affinity refinement.

Theoretical Foundations

The Thermodynamic Cycle and Alchemical Transformations

Alchemical free energy calculations do not directly simulate the physical binding process. Instead, they exploit the concept of a thermodynamic cycle to compute relative binding free energies (ΔΔG) between a wild-type and mutant system or between two similar ligands [26] [24].

Thermodynamic Cycle for Relative Binding Free Energy

Figure 1: The thermodynamic cycle enables the calculation of the relative binding free energy ΔΔG = ΔG_bind1 - ΔG_bind2 by computing the alchemical transformation energies ΔG1 and ΔG2 in the bound and unbound states.

The relative binding free energy is calculated as ΔΔG~bind~ = ΔG~bind,B~ - ΔG~bind,A~ = ΔG~bound~ - ΔG~unbound~, where ΔG~bound~ is the free energy change for transforming Ligand A to Ligand B in the protein's binding site, and ΔG~unbound~ is the equivalent transformation in solution [24].

The three primary alchemical methods share a common goal but differ in their computational approach.

• Free Energy Perturbation (FEP): Based on the Zwanzig equation, FEP estimates the free energy difference between two states by exponentially averaging the energy difference over samples collected from simulations of the initial state [27] [23]. Modern implementations like FEP+ use sophisticated sampling algorithms and force fields to achieve high accuracy [25] [24].
• Thermodynamic Integration (TI): In TI, the free energy difference is computed by integrating the average derivative of the Hamiltonian with respect to a coupling parameter λ, which gradually transforms the system from one state to the other [26]. The integral is typically evaluated numerically over multiple discrete λ windows.
• Bennett Acceptance Ratio (BAR): The BAR method is a potentially more efficient estimator than FEP, as it uses data from both the initial and final states to determine the free energy difference [28] [2]. It is based on the acceptance ratio of Monte Carlo moves between the two states and is considered optimal for minimizing the variance of the estimate for a given amount of sampling [28].

Table 1: Comparison of Key Alchemical Free Energy Methods

Feature	FEP	TI	BAR
Theoretical Basis	Zwanzig equation [23]	Numerical integration of <∂H/∂λ> [26]	Acceptance ratio using data from both states [28]
Key Equation	ΔA = -k~B~T ln<exp(-ΔU/k~B~T)>~A~	ΔA = ∫<∂H/∂λ> dλ	ΔA = k~B~T ln [ <f(U~B~ - U~A~ + C)>~A~ / <f(U~A~ - U~B~ - C)>~B~ ] + C
Sampling Requirement	Samples from one state (traditional)	Samples from multiple λ windows	Samples from both end states (A and B)
Computational Efficiency	High with modern REST-enhanced sampling [24]	Moderate to High, depends on number of λ windows [26]	High, considered statistically optimal [28] [2]
Typical Output	Relative free energy change (ΔΔG)	Relative free energy change (ΔΔG)	Relative free energy change (ΔΔG)

Protocols and Workflows

System Setup and Common Prerequisites

A successful alchemical calculation requires careful initial setup. The following steps are critical across FEP, TI, and BAR:

• Structure Preparation: Protein and ligand structures are obtained from crystallography, homology modeling [24], or computational predictions. Structures must be protonated with correct ionization and tautomeric states assigned using tools like the Protein Preparation Wizard [24].
• Solvation and Embedding: The system is solvated in an explicit water box (e.g., TIP3P) with added ions to neutralize charge and achieve physiological concentration. For membrane proteins like GPCRs, the system must be embedded within an explicit lipid bilayer [2].
• Force Field Assignment: All atoms must be assigned parameters from a modern force field such as OPLS2.1 [24], CHARMM, or AMBER [2]. The choice of force field is a critical determinant of accuracy.
• Alchemical Transformation Definition: The atoms that mutate between the two end-states are identified. The protocol for this transformation can follow a single-topology, dual-topology, or hybrid-topology approach [27].

Detailed Method-Specific Protocols

Free Energy Perturbation (FEP) Protocol

The FEP+ protocol as implemented in Schrödinger's suite represents a state-of-the-art approach [24].

FEP+ Calculation Workflow

Figure 2: Typical workflow for an FEP+ calculation, featuring enhanced sampling and cycle closure for robust results.

• Step 1: Ligand Preparation and Mapping. A congeneric series of ligands is aligned, and a common core structure is defined. Non-core atoms are designated as changing and mapped for the alchemical transformation [24].
• Step 2: λ Window Setup. The transformation is divided into multiple λ windows (often 12-24). In FEP+, the Hamiltonian is scaled differently for the ligand's bonded, van der Waals, and electrostatic terms across these windows.
• Step 3: Enhanced Sampling Simulation. MD simulations are performed at each λ window. The Replica Exchange with Solute Tempering (REST2) [24] or Hamiltonian REMD [26] method is often applied to improve conformational sampling and convergence by allowing exchanges between adjacent λ windows.
• Step 4: Free Energy Analysis and Cycle Closure. The Zwanzig equation is used to calculate the free energy change between adjacent windows, which are summed for the total ΔG. To improve precision, multiple transformations are performed in a closed cycle (e.g., A→B, B→C, C→A), and a consistency correction is applied to minimize the total error around the cycle [24].

Thermodynamic Integration (TI) Protocol

An optimized TI protocol for antibody-antigen binding affinity prediction demonstrates key steps for robustness [26].

• Step 1: System Setup and λ Schedule. A system is built with 12-16 λ windows. A soft-core potential is typically used for van der Waals interactions to avoid singularities as atoms appear or disappear.
• Step 2: Simulation and Hamiltonian REMD. Each λ window is simulated for 3-5 ns. Incorporating Hamiltonian REMD, where replicas at different λ values can exchange, significantly improves convergence and accuracy by overcoming energy barriers [26].
• Step 3: Monitoring and Post-Processing. The derivative ⟨∂H/∂λ⟩ is collected for each window. Spikes in energy derivative (dV/dL) can occur, for instance, when mutating to charged residues. These can be mitigated using a smooth step function and by identifying and excluding outlier windows where the dV/dL deviation between bound and unbound simulations is large (e.g., >1 standard deviation) [26].
• Step 4: Integration. The values of ⟨∂H/∂λ⟩ are integrated numerically (e.g., using the trapezoidal rule) over λ from 0 to 1 to yield the total free energy change.

Bennett Acceptance Ratio (BAR) Protocol

A re-engineered BAR protocol for GPCR-ligand binding affinity involves the following steps [2].

• Step 1: End-State Simulations. Unlike TI, which uses intermediate λ windows, the BAR method typically requires simulations only at the two end-states (λ=0 and λ=1). However, for large perturbations, intermediate states may still be necessary to ensure sufficient overlap.
• Step 2: Energy Collection. The potential energy difference ΔU = U~B~ - U~A~ is collected from trajectories of both state A and state B.
• Step 3: BAR Iteration. The BAR equation is solved iteratively for the free energy difference ΔA and a constant C, which is ideally close to ΔA itself [28]. The optimal function for the BAR method is f(x) = 1 / [1 + exp(x)] [28].
• Step 4: Uncertainty Estimation. The statistical error of the BAR estimate can be computed by bootstrapping or other methods to ensure reliability.

Performance and Validation

The predictive performance of alchemical methods has been rigorously validated across diverse protein classes and ligand series.

Table 2: Performance Benchmarks of Alchemical Methods Across Different Systems

Method / Protocol	Test System	Performance (vs. Experiment)	Key Requirements / Notes	Source
FEP+ (QresFEP-2)	Protein Stability (10 proteins, ~600 mutations)	Excellent accuracy, high computational efficiency	Hybrid-topology; open-source	[27]
FEP+	Diverse Targets (Tyk2, BRD4, A~2A~AR, MCL1)	Pearson R ~0.8, RMSE ~1.0 kcal/mol	Works with crystal structures and homology models (>22% seq. identity)	[24]
Optimized TI	Antibody-Antigen Complexes	Pearson R = 0.74, RMSE = 1.05 kcal/mol (with HREMD)	Requires HREMD for convergence on large systems; 3ns MD per window	[26]
Re-engineered BAR	GPCRs (β~1~AR agonists)	R² = 0.79 with experimental pK~D~	Tailored for membrane proteins; efficient sampling	[2]

Table 3: Key Software, Force Fields, and Computational Resources

Item Name	Type	Function / Application	Examples / Notes
FEP+	Commercial Software Suite	Integrated platform for automated FEP setup, simulation, and analysis.	Schrödinger; industry standard for drug discovery [25] [24]
QresFEP-2	Open-Source Protocol	Hybrid-topology FEP for protein stability and protein-ligand binding.	Integrated with MD software Q; high efficiency [27]
GROMACS	MD Simulation Engine	High-performance open-source software for running MD and free energy calculations.	Used with PMX, TI, and BAR protocols [26] [2]
Desmond	MD Simulation Engine	GPU-accelerated MD engine optimized for FEP+ calculations.	Part of the Schrödinger suite [24]
OPLS Force Fields	Force Field	Family of molecular force fields providing parameters for proteins and ligands.	OPLS2.1 and OPLS3 show high accuracy in FEP+ [24]
CHARMM/AMBER	Force Field	Alternative force fields for biomolecular simulations.	Compatible with various simulation engines like GROMACS [2]
Hamiltonian REMD (HREMD)	Sampling Algorithm	Enhanced sampling technique that exchanges replicas at different λ values.	Critical for improving convergence in TI and FEP [26]
REST2	Sampling Algorithm	Replica Exchange with Solute Tempering; scales the Hamiltonian of a selected region.	Used in FEP+ to improve ligand and binding site sampling [24]

Applications in Drug Discovery

Alchemical methods have moved from theoretical exercises to practical tools that directly impact drug discovery projects.

• Lead Optimization: The primary application is optimizing the potency of lead compounds by predicting the relative binding affinities of proposed analogues before synthesis. This has been successfully applied to numerous targets, including kinases, GPCRs, and protein-protein interactions [24]. For example, FEP+ was used to discover a novel, highly potent A~2A~ adenosine receptor inhibitor [24].
• Antibody Engineering: Optimized TI and FEP protocols can identify beneficial mutations in antibodies to enhance their binding affinity and neutralization potency against antigens, as demonstrated for SARS-CoV-2 neutralizing antibodies [26].
• Protein Engineering and Stability: Protocols like QresFEP-2 are benchmarked on large datasets of protein stability mutations, allowing the prediction of the thermodynamic stability changes resulting from point mutations, which is valuable for enzyme design and understanding genetic diseases [27].
• Selectivity Profiling: These methods can predict a ligand's relative affinity against related off-targets (e.g., within a kinase family), enabling the rational design of selective compounds to minimize side effects [25].

Alchemical free energy calculations using FEP, TI, and BAR have matured into powerful, reliable tools for predicting binding affinities. Driven by advances in force fields, sampling algorithms, and computational hardware, these physics-based methods can now deliver accuracy that rivals experimental measurements. The detailed protocols and benchmarks provided in this application note offer a roadmap for researchers to implement these techniques, thereby accelerating drug discovery and broadening our understanding of molecular recognition. As these methods continue to evolve, their integration with machine learning and application to ever more challenging systems like protein-protein interactions and membrane proteins will further expand their impact on scientific research and therapeutic development.

Molecular recognition, the process by which biomolecules interact with specific partners, is fundamental to nearly all biological processes. The binding affinity between a protein and a small molecule ligand, quantified as the binding free energy (ΔGbind), serves as a crucial determinant in drug discovery and design, dictating the potency and efficacy of potential therapeutic compounds [29]. Among the computational techniques developed to predict this key parameter, end-point methods have emerged as a popular balance between computational efficiency and theoretical rigor. Unlike more computationally intensive alchemical methods that require simulation of intermediate states, end-point methods focus solely on the initial (unbound) and final (bound) states of the binding process [30] [31].

The most widely used end-point approaches are the Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) and Molecular Mechanics Generalized Born Surface Area (MM/GBSA) methods [30]. First introduced by Kollman and colleagues in the late 1990s, these methods combine molecular mechanics energy calculations with implicit solvent models to estimate binding free energies [30] [32]. Their appealing balance of efficiency and reasonable accuracy has led to widespread application in studying protein-ligand interactions, protein-protein complexes, and virtual screening campaigns in drug discovery [30] [32].

This application note details the theoretical foundations, recent methodological enhancements, and practical protocols for applying MM/PBSA and MM/GBSA methods within a research framework focused on molecular dynamics simulations for binding affinity refinement. We provide structured guidelines, validation data, and implementation workflows to assist researchers in leveraging these powerful techniques effectively.

Theoretical Framework

Fundamental Equations

MM/PBSA and MM/GBSA methods estimate the binding free energy (ΔGbind) for a receptor (R) and ligand (L) forming a complex (RL) using the following thermodynamic relationship [30] [32]:

ΔGbind = Gcomplex - (Greceptor + Gligand)

Each free energy term (Gi) is calculated as an ensemble average from molecular dynamics (MD) trajectories and decomposed into specific components:

Gi = ⟨EMM⟩ + ⟨Gsolv⟩ - T⟨S⟩

Where:

• ⟨EMM⟩ represents the average molecular mechanics gas-phase energy
• ⟨Gsolv⟩ denotes the average solvation free energy
• T⟨S⟩ is the entropy contribution at temperature T

The molecular mechanics term is further decomposed as [32]:

EMM = Ebonded + Eelectrostatic + EvdW

Where Ebonded includes bond, angle, and torsion energies; Eelectrostatic represents Coulombic interactions; and EvdW accounts for van der Waals forces.

The solvation free energy is partitioned into polar and non-polar contributions [32]:

Gsolv = Gpolar + Gnon-polar

Key Methodological Variations

MM/PBSA and MM/GBSA differ primarily in their treatment of the polar solvation component:

• MM/PBSA employs numerical solutions to the Poisson-Boltzmann equation for calculating electrostatic solvation energies, generally considered more accurate but computationally demanding [32].
• MM/GBSA utilizes the Generalized Born approximation to estimate electrostatic solvation, offering faster computation with reasonable accuracy for many systems [30] [32].

Both methods typically estimate the non-polar solvation component using a linear relationship with the solvent-accessible surface area:

Gnon-polar = γ × SASA + b

Where SASA represents the solvent-accessible surface area, and γ and b are empirical constants [32].

Recent Methodological Advances

Addressing the Entropy Challenge

A significant limitation of traditional MM/PBSA and MM/GBSA implementations has been the neglect or inadequate treatment of entropy contributions due to the computational expense of normal mode analysis (NMA). Recent work by Dong et al. (2025) has introduced a formulaic entropy approach that computes entropy effects from a single structure based on variations in polar and non-polar solvent-accessible surface areas and rotatable bond counts [33].

This advancement enables incorporation of entropy contributions without incurring additional computational costs. Benchmarking studies demonstrate that integrating formulaic entropy systematically enhances the performance of both MM/PBSA and MM/GBSA, with MM/PBSA_S (including formulaic entropy but excluding dispersion) outperforming other variants across diverse datasets [33].

Performance Across Biological Systems

Recent evaluations have quantified the performance of MM/P(G)BSA methods across various biological targets:

Table 1: Performance of MM/PBSA and MM/GBSA Across Different Systems

System Type	Method	Performance	Conditions	Citation
RNA-ligand complexes	MM/GBSA	Rp = -0.513	GBⁿᵉ², εᵢₙ = 12-16, YIL force field	[34]
RNA-ligand complexes	MM/GBSA	Top-1 success rate: 39.3%	Pose prediction	[34]
General protein-ligand	MM/PBSA_S	Superior performance	With formulaic entropy, no dispersion	[33]
SARS-CoV-2 3CLᵖʳᵒ	MM/PBSA	-100.664 ± 0.691 kJ/mol	For Nirmatrelvir complex	[35]

Membrane Protein Applications

Extending MM/PBSA to membrane proteins presents unique challenges due to the additional complexity of the membrane environment. Recent developments in Amber provide automated membrane parameter calculation and multitrajectory approaches that significantly improve accuracy for systems exhibiting large ligand-induced conformational changes [36].

Zhao et al. (2025) demonstrated this enhanced capability using the human purinergic platelet receptor P2Y12R as a model system, assigning distinct protein conformations (pre- and post-ligand binding) as receptors and complexes in a multitrajectory approach combined with ensemble simulations and entropy corrections [36].

Experimental Protocols

Standard Single-Trajectory Protocol

The most widely used approach employs a single trajectory of the protein-ligand complex:

• System Preparation

  • Obtain 3D structures of the protein-ligand complex from crystallography, docking, or modeling
  • Add hydrogen atoms and assign protonation states appropriate for physiological pH
  • Parameterize the ligand using appropriate force fields (GAFF for small molecules)

• Molecular Dynamics Simulation

  • Solvate the system in explicit water (TIP3P, SPC, etc.)
  • Add ions to neutralize system charge and achieve physiological concentration
  • Energy minimization using steepest descent and conjugate gradient algorithms
  • Gradual heating to target temperature (typically 300K) with position restraints on heavy atoms
  • Equilibration at constant pressure (1 atm) and temperature (300K)
  • Production MD simulation (typically 50-500 ns) with 2 fs time step

• Snapshot Extraction

  • Extract snapshots at regular intervals (typically 50-100 ps) from the stable portion of the trajectory
  • Remove water molecules and ions from each snapshot
  • Ensure adequate sampling by verifying convergence of energy calculations

• Free Energy Calculation

  • Calculate molecular mechanics energies using force fields (AMBER, CHARMM, OPLS)
  • Compute polar solvation energies using PB or GB models
  • Calculate non-polar solvation energies from SASA
  • Estimate entropy contributions using normal mode analysis or formulaic approaches

• Statistical Analysis

  • Calculate ensemble averages and standard errors across all snapshots
  • Perform bootstrapping analysis to assess statistical significance

Specialized Protocol for Membrane Proteins

For membrane protein-ligand systems, the protocol requires modifications to account for the lipid bilayer environment [36]:

• Membrane Placement and System Setup

  • Use specialized tools (Membrane Builder in CHARMM-GUI, PPM server) to embed protein in appropriate lipid bilayer
  • Select lipid composition matching biological context (e.g., POPC for mammalian plasma membranes)
  • Use flexible water models (TIP3P, SPC) and ion parameters compatible with the membrane force field

• Multitrajectory Approach for Conformational Changes

  • Simulate multiple trajectories representing distinct conformational states
  • Assign pre-ligand binding conformations as "receptor" and post-binding as "complex"
  • Combine with ensemble simulations to enhance conformational sampling

• Continuum Membrane Model for MMPBSA

  • Implement heterogeneous dielectric implicit membrane model
  • Automatically determine membrane thickness and location from MD trajectories
  • Ensure consistent treatment of continuum dielectric in electrostatic energy calculations

Entropy Calculation Protocols

Two primary approaches for entropy estimation:

Normal Mode Analysis

• Select representative snapshots from MD trajectory (typically 50-100)
• Minimize each snapshot until root mean square gradient < 0.001 kcal/mol/Å
• Calculate vibrational frequencies for minimized structures
• Compute entropy using quasi-harmonic approximation
• Note: Computationally expensive and may require substantial resources

Formulaic Entropy Approach [33]

• Calculate polar and non-polar SASA for each snapshot
• Count rotatable bonds in the ligand structure
• Apply empirical formula: S = f(SASApolar, SASAnonpolar, Nrotatable)
• Implemented in recent versions of AMBER and other MD packages

Workflow Visualization

MM-PBSA/GBSA Calculation Workflow

Membrane Protein MMPBSA Protocol

The Scientist's Toolkit

Table 2: Essential Research Reagents and Computational Tools

Tool/Resource	Type	Primary Function	Application Notes
AMBER	Software Suite	MD simulations & MMPBSA	Includes MMPBSA.py for automated calculations [36]
CHARMM	Software Suite	MD simulations & analysis	Alternative to AMBER with membrane capabilities
GAFF	Force Field	Small molecule parameterization	General Amber Force Field for drug-like molecules
YIL Force Field	Force Field	RNA-ligand simulations	Specialized for nucleic acid systems [34]
GBⁿᵉ² Model	Solvation Model	Generalized Born solvation	Recommended for RNA-ligand systems [34]
PDBbind	Database	Experimental binding data	Benchmarking and validation [29]
CHARMM-GUI	Web Service	Membrane system preparation	Membrane Builder for bilayer setup [36]

Validation and Best Practices

Method Selection Guidelines

Choosing between MM/PBSA and MM/GBSA depends on the specific research context:

• MM/PBSA is preferred for final, high-accuracy calculations when computational resources permit, particularly for systems with complex electrostatic environments [32].
• MM/GBSA offers advantages for virtual screening and large-scale studies where computational efficiency is prioritized, or for RNA-ligand systems where specific GB models (GBⁿᵉ²) have demonstrated good performance [34].
• The single-trajectory approach is recommended for most systems due to better error cancellation and convergence properties [30].
• The multiple-trajectory approach should be reserved for systems with significant conformational changes upon binding [36].

Statistical Reliability and Error Analysis

Robust implementation requires careful statistical analysis:

• Use bootstrapping or block averaging to estimate standard errors
• Ensure convergence by monitoring cumulative averages of energy components
• Include enthalpy-entropy compensation analysis for complete thermodynamic profiling
• Report Pearson correlation coefficients (Rp) and root mean square errors (RMSE) relative to experimental data

Recent studies indicate that with proper parameterization, MM/GBSA can achieve correlation coefficients of -0.513 to -0.8 with experimental data across diverse systems [33] [34].

MM/PBSA and MM/GBSA methods continue to evolve as valuable tools for binding free energy estimation in drug discovery and molecular recognition studies. Recent advances in entropy calculation, membrane protein applications, and RNA-ligand systems have expanded their applicability and improved their accuracy. The integration of formulaic entropy approaches, automated membrane modeling, and specialized force fields has addressed key limitations while maintaining the computational efficiency that makes these methods particularly valuable for virtual screening and lead optimization campaigns.

When implemented with careful attention to system preparation, sampling adequacy, and statistical validation, these end-point methods provide reliable binding affinity estimates that complement experimental measurements and guide molecular design. Their continued development within the broader context of molecular dynamics simulations for binding affinity refinement promises further enhancements in accuracy, efficiency, and applicability to challenging biological systems.

G Protein-Coupled Receptors (GPCRs) represent one of the most important families of membrane proteins and drug targets, with over 30% of FDA-approved pharmaceuticals targeting these receptors [37]. Predicting the binding affinities and kinetics of GPCR-ligand interactions constitutes a critical challenge in modern drug discovery. Molecular dynamics (MD) simulations have emerged as a powerful computational microscope, providing insights into GPCR-ligand recognition at spatial and temporal scales often inaccessible to experimental techniques [38]. This application note details how MD simulations and complementary computational approaches are revolutionizing the prediction of GPCR-ligand binding parameters within the broader context of affinity refinement research.

The structural flexibility of GPCRs, including pronounced 'breathing' motions observed on nanosecond to microsecond timescales, significantly complicates affinity prediction [39]. Recent advances in dynamic deep transfer learning, efficient sampling algorithms, and large-scale MD datasets have substantially improved our ability to predict and refine binding affinities and understand binding kinetics. These methodologies are increasingly integrated with experimental high-throughput screening to accelerate GPCR drug discovery [40] [37].

Key Methodological Approaches

Researchers employ multiple computational strategies to predict GPCR-ligand binding parameters, each with distinct strengths and applications. The following table summarizes the primary approaches documented in recent literature:

Table 1: Computational Methods for Predicting GPCR-Ligand Binding Affinities and Kinetics

Method	Key Principle	Application in GPCR Research	Performance Highlights
Dynamic Deep Transfer Learning (TrGPCR)	Transfers knowledge from large source domains (e.g., BindingDB) to smaller GPCR-specific datasets (e.g., GLASS database) [41]	Predicts binding affinity using protein sequence and secondary structure features without requiring 3D structure	Improved RMSE by 5.2% and MAE by 4.5% compared to DeepDTA [41]
Molecular Dynamics Refinement	Uses atomistic simulations to refine initial GPCR-ligand models by sampling conformational space [42]	Improves ligand binding mode accuracy and identifies allosteric sites	Improved agreement with experimental ligand binding mode for majority of models [42]
Re-engineered BAR Method	Implements Bennett Acceptance Ratio for efficient sampling and binding free energy calculations [43]	Calculates binding affinities across diverse GPCR targets with varied structural landscapes	Demonstrated correlation with experimental binding affinity data [43]
Large-Scale MD Analysis	Analyzes extensive simulation datasets (e.g., GPCRmd) to uncover conformational dynamics [39]	Reveals allosteric sites, lateral gateways, and breathing motions	Identified lipid insertion sites as markers for allosteric pockets [39]
Free Energy Perturbation (FEP)	Computes relative binding free energies through alchemical transformations [44]	Validates and refines GPCR models by reproducing experimental binding affinities	Accurately predicted binding affinity changes for mutated ligands [44]

Experimental Protocols

MD Refinement of GPCR-Ligand Complexes

Purpose: To improve the accuracy of theoretical GPCR-ligand models through molecular dynamics simulations in a realistic membrane environment.

Materials:

• Initial GPCR-ligand model (from homology modeling or docking)
• Membrane bilayer composition (typically POPC lipids)
• Hydrated simulation system with counterions
• MD simulation software (e.g., GROMACS, NAMD, ACEMD)
• Force field parameters (e.g., CHARMM, OPLS-AA)

Procedure:

• System Setup: Embed the GPCR-ligand complex in an explicit lipid bilayer (e.g., POPC) and hydrate the system with water molecules. Add counterions to neutralize system charge [42] [44].
• Energy Minimization: Perform energy minimization to remove steric clashes using steepest descent or conjugate gradient algorithms.
• Equilibration: Conduct gradual equilibration with position restraints on protein and ligand heavy atoms, first in the NVT ensemble then in the NPT ensemble.
• Production Simulation: Run unrestrained MD simulations for tens to hundreds of nanoseconds. Multiple independent replicates are recommended to improve sampling [42].
• Trajectory Analysis: Cluster simulation snapshots based on ligand RMSD. Compare centroid structures of the largest clusters to reference crystal structures for binding mode validation [42].

Key Considerations: Simulation length should sufficient to capture relevant conformational changes. Weak restraints on transmembrane helices may improve predictions of ligand binding mode and extracellular loop conformations [42].

Binding Affinity Prediction Using Deep Transfer Learning

Purpose: To predict GPCR-ligand binding affinities without requiring protein 3D structures.

Materials:

• Source domain database (e.g., BindingDB)
• Target domain database (e.g., GLASS GPCR database)
• Protein sequence and secondary structure features
• Deep learning framework (e.g., TensorFlow, PyTorch)

Procedure:

• Data Preparation: Extract protein sequences and ligand information from source and target domains. Generate protein pocket features including secondary structure information [41].
• Feature Engineering: Represent proteins using sequence embeddings and secondary structure attributes. Encode ligands using molecular descriptors or fingerprints.
• Model Architecture: Implement a deep neural network with domain adaptation components to transfer knowledge from source to target domain.
• Training: Pre-train the model on the source domain (general protein-ligand complexes), then fine-tune on the target domain (GPCR-specific complexes) using dynamic transfer learning [41].
• Validation: Evaluate model performance using root mean square error (RMSE) and mean absolute error (MAE) on held-out test sets from the GPCR domain.

Key Considerations: Incorporating protein secondary structure (pocket features) significantly improves prediction accuracy compared to sequence-only approaches [41].

Research Reagent Solutions

Table 2: Essential Research Reagents and Computational Tools for GPCR-Ligand Binding Studies

Item	Function/Application	Examples/Specifications
Stabilized GPCR Constructs	Enables cell-free screening with purified receptors; improved stability in micelles [40]	Mutant neurotensin receptor type 1 (NTS1) with enhanced intramolecular stability
Fluorescence Polarization Assay Kits	High-throughput binding screening in 384-well format [40]	DMSO-tolerant systems for library screening
Surface Plasmon Resonance (SPR)	Orthogonal validation of binding hits; measures binding kinetics [40]	Biacore systems for confirming compound binding
GPCRmd Database	Provides curated MD trajectories for comparative analysis and method validation [39]	Contains over 190 GPCR structures with cumulative simulation time >500 μs
Specialized Force Fields	Accurate parameterization of GPCRs, lipids, and ligands for MD simulations [42] [44]	CHARMM27, OPLS-AA with lipid parameters
Genetically Encoded Biosensors	Live-cell reporting of GPCR activation via secondary messengers [37]	GCaMP for calcium, cAMP EPAC sensors

Signaling Pathways and Workflows

Diagram 1: Binding Affinity Prediction Workflow

Diagram 2: GPCR Conformational Dynamics and Pathways

Data Analysis and Interpretation

Performance Benchmarks

Table 3: Quantitative Performance Metrics of GPCR-Ligand Prediction Methods

Method	Dataset/System	Performance Metrics	Comparative Improvement
TrGPCR	GLASS GPCR database	RMSE: 5.2% improvement\nMAE: 4.5% improvement [41]	Superior to DeepDTA baseline
MD Refinement	D3 dopamine receptor models from GPCR Dock 2010	Improved ligand binding mode accuracy for majority of models [42]	Weak TM restraints enhanced EL2 and binding mode prediction
Re-engineered BAR	Multiple GPCR targets	Correlation with experimental binding affinities [43]	Efficient sampling enabled accurate free energy calculations
Large-Scale MD	GPCRmd (190 structures)	Identification of breathing motions: 9.07% intermediate states, 0.5% open states in apo receptors [39]	Revealed lipid insertion sites as conserved allosteric pockets

Technical Considerations

The performance of GPCR-ligand binding affinity prediction is significantly influenced by several technical factors. Data quality and diversity in training sets critically impact model generalizability, with recent studies highlighting how train-test leakage in common benchmarks like PDBbind and CASF has inflated reported performance metrics [45]. Simulation timescale directly affects the ability to capture relevant conformational motions, with GPCR breathing motions occurring on nanosecond to microsecond timescales [39]. Membrane environment representation is essential, as lipid interactions influence ligand binding pathways and allosteric site accessibility [39] [38].

Advanced sampling techniques substantially improve computational efficiency. The re-engineered Bennett Acceptance Ratio method achieves accurate binding free energy calculations through more efficient phase space sampling [43]. Transfer learning approaches address data scarcity issues in GPCR-specific affinity prediction by leveraging larger general protein-ligand datasets [41]. Integration of experimental data, such as mutagenesis studies and binding assays, provides crucial constraints for validating and refining computational models [44].

Aptamers, short single-stranded oligonucleotides, have emerged as powerful alternatives to antibodies in therapeutic and diagnostic applications due to their high specificity, stability, and cost-effectiveness [46]. The detailed understanding of aptamer-protein interaction mechanisms remains a significant challenge in molecular recognition research. This application note presents a case study on employing Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA) calculations to unravel the binding interactions between a truncated DNA aptamer (TrAptm-1) and Macrobrachium rosenbergii nodavirus (MrNV) capsid protein, demonstrating how this computational approach provides quantitative insights into binding affinity and structural basis of molecular recognition [47].

The MM/PBSA method has established itself as a valuable end-point free energy calculation technique that strikes a balance between computational efficiency and theoretical rigor [48]. When applied to aptamer-protein systems, it enables researchers to quantify binding energies, identify key interacting residues, and rationalize experimental observations at the atomic level [49] [47]. This case study details the integration of MM/PBSA within a comprehensive molecular dynamics workflow to investigate how oligomeric state influences aptamer binding affinity, providing a protocol applicable to diverse nucleic acid-protein systems [47].

Background

Aptamer-Protein Interactions

Aptamers are typically selected through Systematic Evolution of Ligands by EXponential enrichment (SELEX), but the conventional process is labor-intensive and often generates numerous candidates requiring costly experimental validation [46]. Molecular dynamics simulations with binding free energy calculations offer a powerful complementary approach to study the structural and dynamic behavior of aptamer-protein interactions at the atomic level [47]. Unlike rigid drug targets, single-stranded nucleic acids like DNA aptamers are inherently flexible molecules with numerous configurations, making molecular dynamics simulation essential for accurately sampling this vast conformational space and understanding their true structure and function [46].

MM/PBSA Methodology

MM/PBSA is an end-point free energy method that calculates binding free energies by combining molecular mechanics energies with implicit solvation models [48]. The binding free energy (ΔGbind) is calculated as follows [48]:

ΔGbind = ΔEMM + ΔGsolv - TΔS

Where ΔEMM represents the gas-phase molecular mechanical energy change, ΔGsolv denotes the solvation free energy change, and -TΔS accounts for the conformational entropy change upon binding [48]. The method is particularly valuable in drug discovery as it offers intermediate accuracy and computational effort between empirical scoring and strict alchemical perturbation methods [30]. For aptamer-protein systems, MM/PBSA provides a computationally efficient approach to evaluate binding affinities and identify key interaction residues, supporting rational aptamer design and optimization [47].

Case Study: TrAptm-1 Aptamer Binding to MrNV Capsid Protein

This case study investigated the binding properties of a truncated DNA aptamer, TrAptm-1, against both monomeric and trimeric forms of the MrNV capsid protein [47]. The primary objectives were to quantitatively compare binding affinities to different oligomeric states, characterize structural rearrangements induced by aptamer binding, and evaluate the antiviral potential of TrAptm-1 to interfere with the capsid protein self-assembly process [47]. The integration of molecular dynamics simulations with MM/PBSA calculations provided a robust framework to address these questions at the molecular level.

Key Findings

Molecular dynamics simulations coupled with MM/PBSA binding free energy calculations revealed that TrAptm-1 exhibited significantly higher binding affinity to the trimeric capsid protein (-153.95 ± 6.74 kcal/mol) compared to the monomeric form (-120.77 ± 2.46 kcal/mol) [47]. This substantial energy difference of approximately 33 kcal/mol indicates strong oligomeric state-dependent binding. TrAptm-1 binding induced significant conformational changes and structural rearrangements in the capsid protein, particularly affecting the protruding (P) domain and shell (S) domain, which are critical for viral capsid assembly [47].

The binding free energy decomposition provided insights into the energetic components driving complex formation. The molecular mechanics component (ΔEMM), comprising van der Waals (ΔEvdW) and electrostatic (ΔEelec) interactions, contributed favorably to binding, while solvation effects (ΔGsolv) presented a significant barrier that was overcome by strong protein-aptamer direct interactions [47]. The stability of the complexes was further validated through root mean square deviation (RMSD) and root mean square fluctuation (RMSF) analyses over 1000 ns simulation trajectories [47].

Table 1: MM/PBSA Binding Free Energy Components for TrAptm-1 Aptamer Complexes (kcal/mol)

System	ΔEvdW	ΔEelec	ΔGpolar	ΔGnonpolar	ΔGbinding
Monomeric MrNV Capsid Protein	-158.32 ± 8.45	-895.64 ± 35.72	941.19 ± 36.15	-7.99 ± 0.55	-120.77 ± 2.46
Trimeric MrNV Capsid Protein	-245.83 ± 12.61	-1124.51 ± 48.93	1235.42 ± 49.87	-19.03 ± 1.15	-153.95 ± 6.74

The study demonstrated that extended simulations up to microseconds are required to capture the slow conformational rearrangements characteristic of large oligomeric protein complexes [47]. These findings provided the molecular basis for developing aptamer-based antiviral strategies and designing biosensors for early detection of MrNV in aquaculture settings [47].

Experimental Protocol

The following workflow outlines the comprehensive protocol for studying aptamer-protein interactions using molecular dynamics simulations and MM/PBSA calculations:

Detailed Methodology

Initial Structure Preparation and Molecular Docking

• Aptamer Structure Preparation: Predict the secondary structure of the truncated aptamer using UNAFold webserver. Based on this secondary structure, construct the three-dimensional structure using 3dRNA web tool. Prior to molecular docking, refine the structure by adding hydrogen atoms [47].

• Protein Structure Preparation: Retrieve the crystal structure of the target protein from the Protein Data Bank (for MrNV capsid protein, PDB ID: 6C6B). Add missing residues using Swiss-PdbViewer software. Perform initial visual inspection of the enzyme structure using UCSF Chimera to remove nonstandard residues [47] [50].

• Molecular Docking: Perform docking using HDOCK, with the refined 3D aptamer structure as the ligand and the target protein as the receptor. Select the aptamer-protein complex with the most stable predicted binding interactions based on the docking energy score for subsequent molecular dynamics simulations [47].

Molecular Dynamics Simulations

• System Preparation: Perform molecular dynamics simulations using GROMACS with the CHARMM36 force field. Refine initial atomic coordinates obtained from molecular docking using CHARMM-GUI. Neutralize the system with counterions as needed. For the MrNV capsid protein-aptamer complex, the final solvated and neutralized systems contained approximately 221,984 atoms for the monomeric complex and 606,104 atoms for the trimeric complex, including solute, water molecules, and counterions [47].

• Energy Minimization: Perform energy minimization using the steepest descent algorithm for 50,000 iterations to ensure steric conflict resolution and structural stability. During this phase, apply no temperature or pressure coupling [47].

• System Equilibration: Equilibrate the system in two successive phases:

  • NVT Equilibration: Equilibrate at 300 K using the velocity-rescale (v-scale) thermostat to stabilize temperature.
  • NPT Equilibration: Equilibrate at 1 bar utilizing the Berendsen barostat to ensure proper system density and pressure stabilization.

  Perform each equilibration phase with a time step of 1 fs (0.001 ps) for 125,000 steps to allow solvent molecules to reach equilibrium [47].

• Production Simulation: Subject the system to production molecular dynamics simulation under NPT ensemble conditions for a duration of 1 μs (1,000 ns), with a time step of 2 femtoseconds (0.002 picoseconds). Apply the v-rescale thermostat to maintain a stable temperature, while the C-rescale barostat regulates pressure at 1 bar. Use the Particle Mesh Ewald (PME) method to compute long-range electrostatic interactions, while calculating van der Waals interactions with a 1.2 nm cutoff. Employ the LINCS algorithm to constrain bond lengths [47].

MM/PBSA Binding Free Energy Calculations

• Snapshot Extraction: Extract a total of 1000 snapshots at 1 ns intervals from the 1000 ns production MD trajectory to ensure statistically robust and representative sampling [47].

• Energy Calculations: Perform energy calculations using the FF99SB force field. Model solvation effects using the Poisson-Boltzmann (PB) approach. Set the solute and solvent dielectric constants to 1.0 and 80.0, respectively. Compute nonpolar solvation contributions using a surface tension model with surften = 0.0072 and surfresc = 0.00 [47].

• Single-Trajectory Approach: Use a single-trajectory approach for consistency across the complex, protein, and aptamer components. Calculate the binding free energy using the following equation [47]:

  ΔGbinding = Gcomplex - (Greceptor + Gligand)

  Where:

  • Gcomplex = free energy of the bound aptamer-protein complex
  • Greceptor = free energy of the unbound viral capsid protein
  • Gligand = free energy of the unbound aptamer

• Energy Decomposition: Further decompose the binding free energy to evaluate contributions from molecular mechanics and solvation effects [47]:

  ΔGbinding = (ΔEvdW + ΔEelec) + (ΔGpolar + ΔGnonpolar)

  Where:

  • ΔEvdW = van der Waals interaction energy
  • ΔEelec = electrostatic interaction energy
  • ΔGpolar = polar solvation free energy
  • ΔGnonpolar = non-polar solvation free energy

The Scientist's Toolkit

Table 2: Essential Research Reagent Solutions for Aptamer-Protein Interaction Studies

Reagent/Software	Function	Application Notes
GROMACS	Molecular dynamics simulation package	Primary software for running MD simulations; compatible with CHARMM36 force field [47]
CHARMM36 Force Field	Empirical energy function for biomolecules	Provides parameters for molecular interactions; suitable for protein-nucleic acid systems [47]
CHARMM-GUI	Web-based graphical user interface	System preparation, solvation, and neutralization; generates proper input files [47]
gmx_MMPBSA	MM/PBSA calculation tool	Implements MM/PBSA method; integrates with GROMACS trajectories [47]
HDOCK	Molecular docking server	Predicts aptamer-protein binding modes; uses refined 3D structures [47]
UNAFold	DNA/RNA folding server	Predicts secondary structure of aptamers; informs 3D structure construction [47]
3dRNA	RNA 3D structure modeling	Constructs three-dimensional aptamer structures from secondary structure [47]
Swiss-PdbViewer	Protein structure analysis	Adds missing residues to protein structures; prepares proteins for simulation [47]
UCSF Chimera	Molecular visualization	Visual inspection of enzyme structures; removes nonstandard residues [47]

Critical Parameters and Troubleshooting

Simulation Quality Control

• Convergence Monitoring: Monitor root mean square deviation (RMSD) of protein and aptamer backbone atoms throughout the simulation to ensure system stability. The MrNV capsid protein-aptamer complex showed RMSD values of approximately 0.5-1.5 nm over 1000 ns simulations [47].

• Sampling Adequacy: Use principal component analysis (PCA) to identify essential dynamics and confirm sufficient sampling of conformational space. For large oligomeric protein complexes, extended simulations up to microseconds are often required to capture slow conformational rearrangements [47].

• Energy Conservation: Monitor potential and kinetic energy throughout the simulation to ensure proper energy conservation and thermostat/barostat function.

MM/PBSA Calculation Parameters

• Dielectric Constants: Set solute and solvent dielectric constants to 1.0 and 80.0, respectively, for polar solvation energy calculations [47].

• Snapshot Selection: Extract snapshots at regular intervals (e.g., every 1 ns) from stable trajectory regions after system equilibration. Avoid using snapshots from the initial equilibration phase [47].

• Entropy Considerations: Note that entropy calculations are computationally demanding and often omitted due to high uncertainty [51]. When included, they are typically estimated via normal-mode or quasi-harmonic analysis [48].

This application note demonstrates the successful application of MM/PBSA methodology to investigate aptamer-protein interactions, using the TrAptm-1-MrNV capsid protein complex as a case study. The protocol outlined provides researchers with a comprehensive framework for quantifying binding affinities, identifying key interaction residues, and rationalizing experimental observations at the atomic level. The integration of molecular dynamics simulations with MM/PBSA calculations offers a powerful approach to accelerate aptamer refinement and candidate validation, effectively reducing the number of sequences requiring expensive experimental testing [46] [47]. As computational resources continue to improve and force fields become more refined, MM/PBSA is poised to play an increasingly important role in the rational design and optimization of aptamer-based therapeutics and diagnostics.

Molecular dynamics (MD) simulations are a cornerstone of computational chemistry and structural biology, providing atomic-level insights into biomolecular interactions, dynamics, and energetics. In binding affinity refinement research, MD simulations enable the precise calculation of free energy differences associated with ligand-receptor binding, which directly correlates with experimental inhibition constants (Ki) and half-maximal inhibitory concentrations (IC50) [31]. The accuracy of these predictions hinges on rigorous system preparation, thorough equilibration, and appropriate production simulations followed by sophisticated analysis. Recent advances have integrated machine learning potentials and automated workflow frameworks to enhance the scalability, efficiency, and accuracy of these simulations, making them more accessible to researchers in drug discovery [52] [53] [54]. This protocol details a comprehensive workflow for MD simulations, with particular emphasis on binding affinity calculation methods like the Bennett Acceptance Ratio (BAR), which has demonstrated superior performance in correlating with experimental binding data for challenging membrane protein targets such as G-protein coupled receptors (GPCRs) [31].

Experimental Workflow and Signaling Pathways

The following diagram illustrates the complete MD simulation workflow for binding affinity studies, from initial system preparation through to final analysis.

Diagram 1: Complete MD simulation workflow from system preparation to binding affinity calculation. The workflow proceeds sequentially through preparation, equilibration, production simulation, and analysis phases, with critical sub-steps at each stage necessary for obtaining reliable binding affinity estimates.

Detailed Methodologies and Protocols

System Preparation Protocol

Objective: Generate physiologically realistic initial structures of the ligand-receptor complex for simulation.

• Structure Preparation: Obtain protein structures from the Protein Data Bank and ligand structures from databases such as PubChem [52]. Process structures using molecular modeling software to add missing hydrogen atoms, assign protonation states, and optimize hydrogen bonding networks.
• Forcefield Selection: Employ biomolecular forcefields such as CHARMM36 or AMBER for protein and ligands. The CHARMM36 force field is particularly well-tested for membrane proteins including GPCRs [52] [54].
• System Building: Place the protein-ligand complex in an appropriate simulation box with adequate padding (typically ≥10 Å from the protein surface). For membrane proteins like GPCRs, embed the system in a lipid bilayer using membrane building tools.
• Solvation and Ionization: Solvate the system with explicit water models (e.g., TIP3P) and add ions to neutralize the system and achieve physiological salt concentration (e.g., 150 mM NaCl) [31] [52].

System Equilibration Protocol

Objective: Gradually relax the system to remove steric clashes and achieve stable temperature and pressure conditions before production simulation.

• Energy Minimization: Perform 50,000 steps of energy minimization using the steepest descent algorithm until the maximum force is < 1000 kJ/mol/nm to remove bad contacts and high-energy configurations [52].
• NVT Equilibration: Equilibrate the system for 100 ps in the NVT ensemble (constant Number of particles, Volume, and Temperature) at 310 K using the V-rescale thermostat (τ = 0.1 ps) with position restraints applied to protein and ligand heavy atoms [52].
• NPT Equilibration: Further equilibrate the system for 100 ps in the NPT ensemble (constant Number of particles, Pressure, and Temperature) at 310 K and 1 bar using the Parrinello-Rahman barostat (τ = 2.0 ps) with maintained position restraints [52].
• Restraint Release: Gradually release position restraints in stages to allow full relaxation of the system while maintaining overall structure.

Production Simulation Protocol

Objective: Generate extensive conformational sampling for statistical analysis of structural dynamics and binding interactions.

• Simulation Duration: Conduct production simulations for 50-100 ns, with longer timescales often required for complex conformational changes or slow binding events [52].
• Integration Parameters: Use a 2 fs time step with LINCS constraints applied to all bonds involving hydrogen atoms to enable longer time steps while maintaining numerical stability [52].
• Electrostatics Treatment: Calculate electrostatic interactions using the Particle Mesh Ewald (PME) method with a Fourier spacing of 1.0 Å for accurate long-range electrostatics [52].
• Data Collection: Save trajectory frames every 10 ps for subsequent analysis, balancing storage requirements with temporal resolution needs [52].

Binding Affinity Calculation Protocol

Objective: Quantitatively estimate ligand binding free energy using alchemical methods.

• Bennett Acceptance Ratio Method: Implement the BAR method, which has shown superior predictive performance in correlating with experimental binding affinities for membrane protein systems [31].
• Lambda Scheduling: Define multiple intermediate states (typically 10-20) represented by scaling factors (lambda values) to smoothly transition between the bound and unstates [31].
• Explicit Solvent Modeling: Use explicit water models for soluble proteins or explicit membrane models for membrane proteins to achieve higher accuracy compared to implicit solvent models, despite the increased computational cost [31].
• Convergence Assessment: Ensure sufficient sampling by monitoring the convergence of free energy estimates over simulation time, potentially requiring re-sampling of new trajectories [31].

Quantitative Data and Parameters

Table 1: Key Simulation Parameters for MD Workflow Stages

Workflow Stage	Key Parameters	Typical Values	Purpose
System Preparation	Forcefield; Water model; Box size; Ion concentration	CHARMM36/AMBER; TIP3P; ≥10 Å padding; 150 mM NaCl	Create physiologically relevant simulation environment
Energy Minimization	Algorithm; Max steps; Force tolerance	Steepest descent; 50,000 steps; <1000 kJ/mol/nm	Remove steric clashes and high-energy configurations
NVT Equilibration	Duration; Temperature; Position restraints	100 ps; 310 K; Applied to protein/ligand	Stabilize temperature while maintaining structure
NPT Equilibration	Duration; Pressure; Barostat	100 ps; 1 bar; Parrinello-Rahman	Stabilize pressure and density to experimental conditions
Production Simulation	Duration; Time step; Constraint algorithm	50-100 ns; 2 fs; LINCS	Generate conformational sampling for analysis
Trajectory Output	Save frequency; Frame interval	Every 10 ps	Balance temporal resolution with storage requirements
Binding Affinity (BAR)	Lambda states; Sampling per window	10-20 states; Adequate for convergence	Calculate free energy differences between states

Table 2: Performance Metrics for Binding Affinity Prediction Methods

Method	Computational Cost	Accuracy	Best Use Cases	Limitations
BAR (Bennett Acceptance Ratio)	High	High (R²=0.7893 with exp. data) [31]	Membrane proteins, GPCRs [31]	Requires extensive sampling; Multiple intermediate states
FEP (Free Energy Perturbation)	High	Moderate to High	Soluble proteins; Lead optimization	Sensitive to initial conditions; Convergence issues
TI (Thermodynamic Integration)	High	Moderate to High	Well-defined binding sites	Numerical integration errors
MM-PBSA/GBSA	Low to Moderate	Moderate	High-throughput screening; Ranking compounds	Implicit solvent limitations; Entropy estimation
LIE (Linear Interaction Energy)	Low	Moderate	Initial screening; Large compound libraries	Empirical parameterization; Transferability issues

Research Reagent Solutions

Table 3: Essential Research Reagents and Computational Tools

Reagent/Tool	Function/Purpose	Example Sources/Implementations
GROMACS	Molecular dynamics simulation package	Open-source package for MD simulation [31] [52]
LAMMPS	Large-scale Atomic/Molecular Massively Parallel Simulator	MD simulation package with MLIP support [54]
CHARMM36 Force Field	Empirical energy function for biomolecules	Parameter set for proteins, lipids, nucleic acids [52] [54]
PyTorch ML Models	Machine learning interatomic potentials	ML-IAP-Kokkos interface for LAMMPS [54]
AutoDock Vina	Molecular docking and virtual screening	CB-Dock2 platform for protein-ligand docking [52]
PubChem Database	Chemical structure and bioactivity database	Source for ligand structures and properties [52]
SwissTarget Prediction	Target prediction for bioactive compounds	Database for potential protein targets [52]
DynaMate Framework	AI-agent for automating simulation workflows	Multi-agent framework for MD setup and analysis [53]

Advanced Integration and Automation

The field of molecular dynamics simulation is rapidly evolving with the integration of machine learning approaches and automated workflow frameworks. The ML-IAP-Kokkos interface enables fast and scalable MD simulations by integrating PyTorch-based machine learning interatomic potentials with the LAMMPS MD package, facilitating end-to-end GPU acceleration [54]. This interface uses Cython to bridge Python and C++/Kokkos LAMMPS, ensuring optimal performance on NVIDIA GPU architectures [54].

Furthermore, frameworks like DynaMate leverage AI agents to automate repetitive tasks in molecular simulations, providing a modular template for building multi-agent systems that can generate, run, and analyze MD simulations through natural language prompts [53]. This approach significantly reduces the barrier for researchers to implement complex simulation workflows and enhances reproducibility in binding affinity studies.

For binding affinity refinement research, these advanced computational approaches provide unprecedented ability to predict and optimize drug-target interactions, with particular value for challenging target classes such as GPCRs, which represent approximately 34% of all approved drugs [31]. The continuous development of more accurate forcefields, enhanced sampling algorithms, and automated workflow integration promises to further accelerate drug discovery pipelines in the coming years.

Overcoming Computational Hurdles: Strategies for Robust and Accurate Affinity Predictions

In molecular dynamics (MD) simulations, the accurate prediction of biomolecular processes, such as protein-ligand binding, is fundamentally constrained by the insufficient sampling of rare events. Conventional MD simulations often fail to adequately explore the energy landscape because high energy barriers can trap the simulation in local minima, preventing the observation of biologically relevant transitions. This is particularly problematic for binding affinity predictions, where the calculated free energies often fail to validate experimental results due to this limited sampling [31]. Enhanced sampling techniques, such as metadynamics and umbrella sampling, directly address this limitation by accelerating the exploration of the free energy landscape along carefully chosen degrees of freedom, thereby providing more reliable and computationally efficient estimates of thermodynamic and kinetic properties.

Enhanced Sampling Methodologies: Core Principles and Protocols

Metadynamics

Metadynamics is an enhanced sampling method that facilitates the escape from energy minima by employing a history-dependent bias potential. This bias is constructed as a sum of Gaussian functions deposited along predefined Collective Variables (CVs), which are descriptors of the system's slowest degrees of freedom (e.g., a key protein-ligand distance) [55]. Over time, this bias "fills up" the visited energy basins, forcing the system to explore new regions and allowing for the reconstruction of the underlying free energy surface.

Protocol for Determining GPCR-Ligand Binding Modes using Metadynamics [55]

The following workflow outlines a protocol for predicting the binding mode of a ligand to a G Protein-Coupled Receptor (GPCR), a major pharmaceutical target.

• Step 1: System Preparation.

  • Obtain the GPCR structure, often from a crystal structure. The starting structure can be crystallized with a different ligand or in a different activation state.
  • Place the unbound ligand of interest at a random position within the solvent above the receptor.
  • Equilibrate the system using conventional MD (e.g., microsecond-scale simulations) in the absence of the ligand to relax the protein structure [55].

• Step 2: Define the Collective Variable (CV).

  • A common and effective CV for GPCR-ligand binding is the z-component of the distance between a key atom in the ligand (e.g., the ammonium nitrogen) and the Cα atom of the conserved residue Trp6.48 at the bottom of the binding pocket [55].
  • To improve sampling efficiency, a funnel-shaped restraint can be applied to keep the ligand within a defined region near the binding path, preventing it from drifting into the bulk solvent [55].

• Step 3: Perform Initial Metadynamics.

  • Conduct one or several short metadynamics simulations starting from the unbound state. The goal is to observe multiple entries of the ligand into the orthosteric binding site. For a system like histamine binding to the H1 receptor, 25 ns may be sufficient to observe reversible binding events [55].

• Step 4: Set Up Multiple Walker Metadynamics.

  • Extract sample conformations from the initial trajectories along the entire binding path.
  • Use these conformations as starting points for a multiple walker metadynamics simulation. In this approach, several simulations ("walkers") run in parallel, sharing a common bias potential, which leads to more efficient exploration and convergence of the free energy profile [55].

• Step 5: Analyze Results and Identify Binding Mode.

  • Calculate Free Energy Profile: From the multiple walker trajectories, compute the free energy as a function of the CV. The global minimum of this profile corresponds to the most stable bound state.
  • Cluster Analysis: Extract all simulation frames located within a small interval (e.g., ±3 Å) around the CV value of the global minimum. After aligning the protein backbone, cluster these frames based on the ligand's conformation. The most populated cluster often represents the preferred binding mode [55].
  • Refinement with MD (Optional): To allow for final conformational adjustments, run short, conventional MD simulations starting from the representative structure of the top cluster(s). The most stable pose from this unbiased simulation is the final predicted binding mode.

Umbrella Sampling

Umbrella Sampling is a powerful method for calculating the free energy profile along a specific reaction coordinate. The core idea is to divide the reaction path into multiple overlapping "windows." In each window, a harmonic restraint potential is applied to maintain the system at a specific value of the reaction coordinate. By simulating all windows independently and then combining the data, a complete picture of the free energy landscape is obtained.

Protocol for Calculating Free Energy of Membrane Insertion using Umbrella Sampling [56]

This protocol is adapted from studies calculating the free energy of transfer of small solutes into a lipid membrane.

• Step 1: Define the Reaction Coordinate.

  • For a process like membrane insertion, the reaction coordinate is typically the distance between the center of mass of the solute (e.g., a polyethylene oligomer) and the center of mass of the lipid bilayer.

• Step 2: Generate Initial Configurations.

  • Perform a preliminary simulation, such as a "pull" simulation, to generate a series of configurations along the reaction coordinate, from the water phase to the membrane center.

• Step 3: Set Up and Run Umbrella Windows.

  • Window Selection: Define a set of simulations (windows) that span the entire reaction coordinate, ensuring sufficient overlap between the probability distributions of adjacent windows.
  • Apply Restraint: In each window, run an MD simulation with a harmonic potential (the "umbrella") applied to the reaction coordinate. The force constant of the restraint must be strong enough to ensure good sampling within the window but allow for overlap with neighboring windows.

• Step 4: Free Energy Reconstruction.

  • Use the Weighted Histogram Analysis Method (WHAM) to combine the data from all umbrella sampling windows. WHAM removes the bias introduced by the restraint potentials and yields the unbiased free energy profile along the reaction coordinate.

Comparative Analysis and Practical Application

Quantitative Comparison of Techniques

The choice between enhanced sampling methods depends on the scientific question and available computational resources. The table below summarizes a direct comparison between Metadynamics and Umbrella Sampling for a specific benchmark system.

Table 1: Comparison of Metadynamics and Umbrella Sampling for calculating the free energy of transfer of small solutes into a model lipid membrane [56].

Feature	Metadynamics	Umbrella Sampling
Core Principle	History-dependent bias (Gaussians) fills visited states.	Harmonic restraints applied in sequential windows.
Free Energy Estimate	Directly from the applied bias potential.	Post-processing via WHAM.
Statistical Efficiency	More efficient in the tested system; yielded lower statistical uncertainties for the same simulation time [56].	Less efficient for the same simulation time.
Result Quality	Yielded the same free energy profile estimate as Umbrella Sampling [56].	Yielded the same free energy profile estimate as Metadynamics [56].
Key Advantage	Efficient exploration without prior knowledge of the full path.	Well-established, robust method for obtaining profiles along a pre-defined coordinate.

Application in Binding Affinity Calculation

Beyond identifying binding poses, enhanced sampling is crucial for calculating binding free energies. Alchemical methods, such as the Bennett Acceptance Ratio (BAR), are widely used for this purpose. A recent study addressed sampling limitations in BAR calculations for membrane proteins like GPCRs by developing a protocol that achieves efficient sampling through re-engineering the classic BAR method [31]. This BAR-based approach demonstrated a significant correlation (R² = 0.7893) with experimental binding affinity data (pK_D) for ligands like isoprenaline and salbutamol bound to the β1 adrenergic receptor (β1AR), validating its performance [31]. This highlights how improving sampling directly enhances the accuracy of binding affinity predictions.

Table 2: Key research reagents and computational tools for implementing enhanced sampling protocols.

Item / Resource	Function / Description	Relevance to Protocol
GPCR Crystal Structure	Provides the initial atomic coordinates of the target protein.	Starting point for system setup; can be from a different ligand/state [55].
Collective Variable (CV)	A low-dimensional descriptor of the process of interest (e.g., binding, insertion).	Critical for guiding metadynamics; a poor CV leads to inefficient sampling [55].
Funnel Restraint	A soft boundary potential that keeps the ligand near the binding site.	Increases sampling efficiency in metadynamics by preventing ligand escape [55].
Multiple Walker Algorithm	Parallel simulations that share and contribute to a common bias potential.	Accelerates convergence of metadynamics simulations [55].
Weighted Histogram Analysis Method (WHAM)	A statistical method for combining data from biased simulations.	Essential for unbinding the data from umbrella sampling windows to get the free energy profile [56].
Molecular Dynamics Engine	Software to perform the simulations (e.g., GROMACS, CHARMM, AMBER).	Execution platform for both conventional and enhanced sampling simulations [31].

Enhanced sampling techniques like metadynamics and umbrella sampling are indispensable for overcoming the timescale limitations of conventional MD, enabling accurate determination of binding modes and free energies. The developed protocols provide a clear roadmap for researchers to apply these methods to biologically critical systems such as GPCRs. The field continues to evolve rapidly, with a major future direction being the integration of machine learning. ML techniques are now being used to data-drive the construction of optimal collective variables, improve biasing schemes, and even enable automated strategies for rare-event sampling through reinforcement learning and generative approaches [57]. This synergy between physical simulation and machine learning promises to further automate and enhance the predictive power of molecular simulations in drug discovery and biomolecular research.

Accurate prediction of protein-ligand binding affinity is a critical objective in computational drug discovery. The selection of an appropriate force field—the mathematical function and parameter set that describes the potential energy of a molecular system—serves as the foundation for reliable molecular dynamics (MD) simulations and binding free energy calculations. Force fields inherently represent a balance between computational efficiency and physical accuracy, with different classes of force fields occupying distinct positions on this spectrum. Traditional molecular mechanics (MM) force fields, which use fixed atomic charges and pre-defined bonding terms, enable the simulation of large biological systems over microsecond timescales but may lack the quantum mechanical precision necessary for accurately modeling bond formation and breaking or complex electronic phenomena. Recent advances have introduced more sophisticated potentials, including reactive force fields (ReaxFF) and neural network potentials (NNPs), which bridge the accuracy gap while introducing new computational considerations. This application note provides a structured framework for force field selection and validation, with a specific focus on optimizing this process for binding affinity refinement within drug discovery pipelines.

Force Field Landscape and Performance Comparison

The current landscape of force fields and binding affinity prediction methods spans a wide range of the accuracy-computational cost spectrum. Understanding the performance characteristics of each approach is essential for making informed selections based on project requirements and resource constraints.

Table 1: Comparison of Binding Affinity Prediction Methods and Their Performance Characteristics

Method	Accuracy (RMSE)	Correlation with Experiment	Computational Cost	Typical Use Case
Molecular Docking [51]	2–4 kcal/mol	~0.3	Minutes on CPU	High-throughput virtual screening
MM/PBSA & MM/GBSA [51]	Variable, often high	Variable, can be low	Medium (hours-days)	Post-docking rescoring, trend analysis
Free Energy Perturbation (FEP) [51] [58]	~1 kcal/mol	0.65+	High (12+ hours GPU per calculation)	Lead optimization, relative binding affinity
Bennett Acceptance Ratio (BAR) [31]	High (correlation R²=0.79 shown)	High	High (similar to FEP)	Membrane protein targets (e.g., GPCRs)
Neural Network Potentials (NNPs) [59]	DFT-level accuracy (e.g., MAE ~0.1 eV/atom)	High for trained systems	Lower than DFT, higher than classical MD	Reactive processes, material properties
Physics-Informed ML [58]	Comparable to FEP	High	~1000x less than FEP	Early-stage screening, diverse chemical space

The choice of method often involves a direct trade-off between speed and accuracy. Docking provides rapid results but with significant error, making it suitable for initial screening but unreliable for quantitative predictions [51]. At the other extreme, rigorous alchemical methods like FEP and BAR provide high accuracy and strong correlation with experimental binding affinities but require substantial computational resources, limiting their throughput [51] [31]. The middle ground has traditionally been occupied by methods like MM/GBSA, but their predictive performance can be variable. Emerging approaches, such as physics-informed machine learning and neural network potentials, aim to disrupt this trade-off by offering near-FEP accuracy at a fraction of the cost [59] [58].

Table 2: Comparison of Force Field Types and Their Characteristics

Force Field Type	Physical Rigor	Transferability	Computational Cost	Key Applications
Classical MM (e.g., AMBER, CHARMM)	Medium	High	Low	Standard protein-ligand MD, MM/PBSA
Reactive (ReaxFF) [59] [60]	Medium-High	Low-Medium	Medium	Chemical reactions, bond breaking/formation
Neural Network (NNP) [59]	High (DFT-level)	Medium (with training)	Medium-High	High-accuracy energy/force prediction
Physics-Informed ML [58]	High (Physics-grounded)	High	Very Low (after training)	High-throughput binding affinity prediction

Force Field Validation Methodologies

Benchmarking against Experimental and High-Level Theoretical Data

Validation is a critical step to ensure that a selected force field is appropriate for the biological system and property of interest. A robust validation protocol involves multiple layers of comparison.

For binding affinity prediction, the primary validation metric is the correlation between calculated and experimental binding free energies (ΔG). As demonstrated in studies on GPCRs, a well-validated protocol using the BAR method can achieve a strong correlation (R² = 0.79) with experimental pKD values [31]. Similarly, for NNPs like EMFF-2025, validation involves comparing predicted energies and forces against Density Functional Theory (DFT) calculations, with benchmarks showing mean absolute errors (MAE) for energy within ±0.1 eV/atom and for forces within ±2 eV/Å [59].

Protein-ligand complex structures should also be validated by assessing the stability of key interactions observed in crystal structures during MD simulations. For instance, in β1 adrenergic receptor simulations, the persistence of hydrogen bonds with residues like S211 and S215 can be monitored as an indicator of simulation quality [31].

Addressing Data Dependency and Generalization

A critical, often overlooked aspect of validation is assessing a method's performance on unseen data. Machine learning-based approaches, in particular, are susceptible to data leakage and overfitting if not carefully validated. Studies show that models achieving high correlation (Pearson R up to 0.70) under random data partitioning can see significantly degraded performance when evaluated using more rigorous, protein-based partitioning schemes that better simulate real-world prediction scenarios [61]. Therefore, validation should always use strict data separation rules, such as UniProt-based splits or the anchor-query framework, which leverages limited reference data to improve predictions for novel proteins [61].

Detailed Experimental Protocols

Protocol 1: MM/GBSA Calculation for Binding Affinity Ranking

This protocol outlines the steps for performing an MM/GBSA calculation to estimate relative binding affinities, based on insights from failed attempts that highlight critical pitfalls [51].

• System Preparation:

  • Obtain the protein-ligand complex structure from docking or experimental sources.
  • Prune the protein to a fixed radius (e.g., 10-12 Å) around the ligand to reduce computational cost. Retain crystallographic waters if they are part of the binding site.
  • Parameterize the ligand using a tool like antechamber (from AMBER tools) with the GAFF force field and assign partial charges (e.g., using the AM1-BCC method).
  • Add hydrogen atoms to the protein and assign protonation states of titratable residues (e.g., using H++ or PROPKA) appropriate for the physiological pH.

• Molecular Dynamics Simulation:

  • Solvate the system in a TIP3P water box with a minimum 10 Å buffer. Add ions to neutralize the system and achieve a physiological salt concentration (e.g., 150 mM NaCl).
  • Perform energy minimization in two stages: first, restraining the heavy atoms of the protein and ligand; second, a full minimization without restraints.
  • Gradually heat the system from 0 K to 300 K over 50-100 ps under constant volume (NVT ensemble) with restraints on solute heavy atoms. Failure to heat gradually can cause large initial forces and convergence issues [51].
  • Equilibrate the system at 300 K and 1 atm for at least 100 ps (NPT ensemble) without restraints.
  • Run a production MD simulation for a sufficient duration (e.g., 4-20 ns, NPT ensemble). The length depends on system size and convergence.

• Trajectory Processing and Snapshot Extraction:

  • After discarding the equilibration phase, extract snapshots from the production trajectory at regular intervals (e.g., every 10-100 ps). For the system described in [51], 300 snapshots were extracted.
  • Process the trajectory to correct for periodicity and align frames to a reference (e.g., the protein backbone) to remove global rotation and translation. Tools like PyTraj (for AMBER) or MDTraj can be used.

• Free Energy Calculation:

  • For each snapshot, calculate the gas-phase interaction energy (ΔH_gas) using the molecular mechanics force field.
  • Calculate the solvation free energy (ΔG_solvent), which is decomposed into:
    • Polar contribution: Solved using an implicit solvent model like Generalized Born (GB). This can be done using MD packages like OpenMM or MMPBSA.py (AMBER).
    • Non-polar contribution: Approximated as a linear function of the Solvent-Accessible Surface Area (SASA), calculated using an algorithm like Shrake-Rupley (available in MDTraj).
  • The binding free energy for each snapshot is estimated as: ΔG ≈ ΔHgas + ΔGsolvent. The entropic term (-TΔS) is computationally demanding and noisy to calculate and is often omitted for relative ranking [51].
  • Average the ΔG values over all snapshots to obtain the final estimate.

Workflow for MM/GBSA Binding Affinity Calculation

Protocol 2: Validation of Force Field Parameters using Bayesian Inference

This protocol describes the use of Bayesian Inference of Conformational Populations (BICePs) for force field validation and refinement against ensemble-averaged experimental data, such as NMR observables [62].

• Define the Prior and Observables:

  • Prior Generation: Run molecular simulations (e.g., MD) using the force field to be validated. This generates a conformational ensemble, which serves as the prior distribution, p(X).
  • Experimental Data: Collect ensemble-averaged experimental measurements (D), such as J-couplings or chemical shifts from NMR. Define the forward model, f(X), which is a function that computes the theoretical value of the observable from a molecular structure.

• BICePs Sampling:

  • Configure the BICePs algorithm with a likelihood function (e.g., Gaussian or Student's t-model) and a non-informative Jeffreys prior, p(σ), for the uncertainty parameters.
  • Use a replica-averaged forward model. The number of replicas (N_r) should be large enough to provide a good estimate of the ensemble average.
  • Perform Markov Chain Monte Carlo (MCMC) sampling to sample the posterior distribution p(X, σ | D), which includes both conformational states and the uncertainties.

• Calculate the BICePs Score and Analyze:

  • Compute the BICePs score, which is a free energy-like quantity that reflects the total evidence for the model given the data. A more negative BICePs score indicates a better agreement between the force field's prior ensemble and the experimental data.
  • Analyze the posterior distributions of the uncertainties (σ). Parameters with high inferred uncertainty may be outliers or subject to systematic error.
  • The BICePs score can be used for model selection (choosing between different force fields) or as an objective function for automated parameter refinement via variational optimization [62].

Workflow for Bayesian Force Field Validation (BICePs)

Protocol 3: High-Throughput Screening with Physics-Informed Machine Learning

This protocol leverages physics-informed machine learning models for high-throughput binding affinity prediction, as a cost-effective pre-screening step before FEP [58].

• Model Selection and Setup:

  • Select a physics-informed ML method that uses multiple-instance learning and explicitly models physical factors like shape, electrostatics, and hydrogen bonding [58].
  • Prepare a database of protein-ligand complexes with experimentally determined binding affinities for model training and validation. Ensure strict data partitioning (e.g., by protein sequence) to avoid overfitting [61].

• Ligand Pose Generation and Affinity Prediction:

  • For a new ligand, the model dynamically identifies and refines optimal ligand poses within the protein binding pocket, analogous to molecular docking.
  • The model outputs a predicted binding affinity. Due to its low computational cost (approximately 0.1% of FEP), this can be applied to thousands or millions of compounds [58].

• Triaging and Refinement:

  • Rank all screened compounds based on their predicted affinity.
  • Select the top-ranking candidates (e.g., a few hundred) for more accurate and computationally expensive validation using FEP or BAR methods. This sequential approach allows for the exploration of a much wider chemical space with the same computational budget [58].

The Scientist's Toolkit: Essential Research Reagents and Materials

Table 3: Key Software Tools and Resources for Force Field Applications

Tool Name	Type/Category	Primary Function	Application Context
GROMACS [31]	MD Simulation Engine	High-performance molecular dynamics	Running production MD simulations for trajectory generation.
AMBER [31]	MD Simulation & Analysis Suite	Molecular dynamics, parameterization, MM/PBSA	System setup, simulation, and endpoint free energy calculations.
OpenMM	MD Simulation Library	GPU-accelerated molecular dynamics	Rapid energy minimization, equilibration, and production MD.
Deep Potential (DP) [59]	Neural Network Potential Framework	Training and running NNPs	Performing MD with DFT-level accuracy for complex reactions.
BICePs [62]	Bayesian Inference Software	Force field validation and refinement	Assessing and optimizing force fields against experimental data.
Optibrium (StarDrop) [58]	Drug Discovery Platform	Physics-informed ML for affinity prediction	High-throughput binding affinity screening and prioritization.
Rosetta (flex ddG) [63]	Protein Design Suite	Predicting binding affinity changes upon mutation	Assessing impacts of protein mutations on protein-peptide interactions.
DP-GEN [59]	Automated Training Workflow	Generating general-purpose NNPs	Automating the exploration and training of neural network potentials.
PLINDER	Dataset	Curated protein-ligand complex data	Training and benchmarking machine learning models for binding affinity [51].
BindingDB	Database	Experimental binding data	Source of experimental IC50/KD values for validation [51].

The effective selection and validation of force fields are paramount for achieving reliable results in molecular dynamics simulations aimed at binding affinity refinement. No single force field or method is universally superior; the optimal choice is a strategic decision based on the specific question, system, and available computational resources. The protocols and analyses provided here offer a practical roadmap for researchers. By rigorously benchmarking methods against relevant experimental data, employing robust validation techniques like BICePs to avoid overfitting, and strategically combining high-throughput physics-informed ML with rigorous FEP calculations, scientists can significantly enhance the efficiency and predictive power of their drug discovery pipelines. The future of force field development lies in the continued integration of physical principles with data-driven machine learning approaches, promising further improvements in both accuracy and computational accessibility.

Molecular dynamics (MD) simulation is a pivotal tool in modern computational biology and drug discovery, providing atomic-level insight into the behavior of biomolecular systems. However, the study of biological processes such as protein folding, ligand binding, and cellular-scale phenomena requires the examination of systems and timescales that are often prohibitively expensive for all-atom molecular dynamics (AAMD). Drug action is inherently multiscale, connecting molecular interactions to emergent properties at cellular and larger scales. Coarse-grained molecular dynamics (CG-MD) addresses this challenge by reducing the number of degrees of freedom in the system, thereby extending the accessible spatial and temporal scales by orders of magnitude [64] [65]. In CG-MD, groups of atoms are consolidated into single "beads" or interaction sites, creating a simplified representation that captures essential physics while enabling the simulation of biologically relevant scales [65].

Multiscale modeling represents a complementary approach that explicitly connects different levels of biological organization. These methods aim to model and analyze how changes on one scale lead to changes at another, creating a fluid knowledge landscape that synthesizes disparate modeling and experimental data across spatiotemporal scales [64]. Multiscale techniques can be constructed in different ways, including combining different levels of theory or resolution (e.g., QM/MM, MD/BD) or through the hierarchical integration of sets of approaches carried out at different scales [64]. The convergence of improved biophysical data, advanced algorithms, and increasingly powerful computing architectures is driving advances in multiscale modeling methods capable of bridging chemical and biological complexity from the atom to the cell [64].

Theoretical Foundations of Coarse-Graining

Physical Basis and Equations of Motion

The motion of coarse-grained sites is governed by the potential of mean force and the friction and stochastic forces that result from integrating out the secondary degrees of freedom [65]. The physical basis for coarse-grained force fields stems from statistical mechanics, where the effective interactions between CG sites are derived from the underlying all-atom system. According to the Mori-Zwanzig formalism, the equation of motion for CG degrees of freedom can be expressed as:

[ \frac{d}{dt}Pi(t) = -\nabla{Ri}W(R) - \frac{1}{kB T} \sumj \int0^t \langle \delta Fi^Q(\tau) \delta Fj^Q(\tau) \rangle \cdot vi(\tau) d\tau + \delta Fi^Q(t) ]

where ( W(R) ) is the potential of mean force, the second term represents friction forces, and the final term corresponds to stochastic forces [65]. In practical implementations, the friction and stochastic force terms are often simplified using Langevin dynamics, which assumes that fine-grained degrees of freedom move much faster than coarse-grained ones [65].

Machine Learning Revolution in Coarse-Graining

Recent advances have integrated deep learning methods with large and diverse training sets of all-atom protein simulations to develop bottom-up CG force fields with chemical transferability. These machine-learned models can be used for extrapolative molecular dynamics on new sequences not used during model parameterization, demonstrating the feasibility of universal and computationally efficient CG models for proteins [66]. By leveraging the variational force-matching approach, such force fields have been shown to accurately reproduce all-atom distributions of CG observables, enabling the simulation of complex biomolecular processes with accuracy approaching all-atom methods but at substantially lower computational cost [66].

Table 1: Comparison of Molecular Dynamics Approaches

Feature	All-Atom MD (AAMD)	Traditional CG-MD	Machine-Learned CG-MD
Resolution	Atomic detail	Beads representing multiple atoms	Beads representing multiple atoms
Computational Speed	Reference (1x)	10-1000x faster than AAMD	Several orders of magnitude faster than AAMD
Accuracy	High for small systems	System-specific, limited transferability	High, with demonstrated transferability
Training Requirement	Not applicable	Parameterized for specific systems	Requires training on diverse all-atom datasets
Typical Applications	Ligand binding, enzyme catalysis	Membrane formation, large complexes	Protein folding, conformational dynamics, binding

Key Coarse-Grained Models and Force Fields

Established Coarse-Grained Force Fields

Several CG force fields have been developed with different strengths and applications. The MARTINI force field, developed by Marrink's group, is particularly effective for modeling intermolecular interactions including membrane structure formation and protein interactions, though it inaccurately models intramolecular protein dynamics [66] [65]. UNRES and AWSEM are CG force fields developed specifically to model protein folding and conformational dynamics, though they often struggle to capture alternative metastable states [66]. The UNICORN model of biological macromolecules represents another approach developed for specific laboratory applications [65].

These traditional CG models face several challenges. The most significant is the difficulty in efficiently modeling multi-body interaction terms, which are essential to realistically represent correct protein thermodynamics and implicit solvation effects [66]. In contrast, classical all-atom force fields model most non-bonded interactions as a sum of two-body terms, simplifying computation but potentially limiting accuracy for complex biomolecular processes.

Emerging Machine-Learned Coarse-Grained Models

Recent advances in machine learning have enabled the development of neural network-based CG models that are truly transferable in sequence space. These models are learned using bottom-up approaches from atomistic simulations of a set of proteins and can successfully simulate conformational dynamics of proteins never seen during any learning stage, even with low (16-40%) sequence similarities to training proteins [66]. For example, the CGSchNet model demonstrates the ability to predict metastable states of folded, unfolded and intermediate structures, fluctuations of intrinsically disordered proteins, and relative folding free energies of protein mutants, while being several orders of magnitude faster than all-atom models [66].

Table 2: Performance Metrics of Machine-Learned CG Model on Test Proteins

Protein System	Fraction of Native Contacts (Q) in Folded State	Cα RMSD to Native (nm)	Ability to Predict Metastable States
Chignolin (2RVD)	Close to 1	Low	Yes, including misfolded state
TRPcage (2JOF)	Close to 1	Low	Yes
BBA (1FME)	~0.75	~0.5 nm	Partial (difficulty with β-sheet motifs)
Villin Headpiece (1YRF)	Close to 1	Low	Yes
Engrailed Homeodomain (1ENH)	~0.75	~0.5 nm	Yes, with similar terminal flexibility to AAMD

Application Notes: Binding Affinity Predictions

Alchemical Methods for Binding Free Energy Calculations

In binding free energy calculations, the use of an explicit solvent model for soluble proteins or a membrane model for membrane proteins generally provides higher accuracy than implicit solvent models, though at significantly increased computational cost [31]. Various computational methods have been utilized to estimate binding free energies, including the Linear Interaction Energy (LIE) method that focuses on simple endpoint energies, implicit solvent approaches such as MM-PB/GBSA, and alchemical perturbation methods such as Free Energy Perturbation (FEP), Thermodynamic Integration (TI), and Bennett Acceptance Ratio (BAR) [31].

The BAR method has shown particular promise for predicting binding affinities of ligands to membrane protein targets such as G-protein coupled receptors (GPCRs). Recent implementations have demonstrated significant correlation with experimental pK_D values (R² = 0.7893) for diverse β1 adrenergic receptor ligands, including full agonists, partial agonists, and weak partial agonists [31]. This approach successfully captured the differential binding affinity of selective agonists to active versus inactive states of β1AR, demonstrating its sensitivity to receptor conformational states.

Protocol: BAR-Based Binding Free Energy Calculations for Membrane Proteins

Objective: Predict binding affinity of small molecules to membrane protein targets using the Bennett Acceptance Ratio method.

Step 1 - System Preparation

• Obtain protein structure from experimental sources or AlphaFold2 prediction [67]
• For membrane proteins: embed in appropriate lipid bilayer using insane or CHARMM-GUI protocols
• Add explicit solvent (e.g., TIP3P water) and ions to neutralize system
• Generate ligand parameters using appropriate tools (e.g., CGenFF, GAFF)

Step 2 - Equilibration

• Energy minimization using steepest descent algorithm (5,000-10,000 steps)
• Gradual heating from 0 to 300 K over 100-200 ps with position restraints on protein and ligand
• Equilibrium MD without restraints for 10-20 ns until system properties stabilize

Step 3 - Lambda Sampling Setup

• Divide the perturbation range into 10-20 intermediate states (λ values)
• For each λ state, generate initial coordinates and topology
• Set up simulations for both forward and backward transitions

Step 4 - Production Simulations

• Run parallel simulations for each λ window (1-10 ns each)
• Employ replica exchange or Hamiltonian replica exchange if possible
• Save trajectories at regular intervals for analysis

Step 5 - Free Energy Analysis

• Calculate free energy differences using BAR estimator between adjacent λ windows
• Sum contributions across all windows to obtain total binding free energy
• Perform error analysis using bootstrapping or block averaging

Step 6 - Validation

• Compare with experimental binding data if available
• Analyze convergence by examining time-dependent free energy estimates
• Calculate thermodynamic cycles to check for consistency

Application Notes: Protein-Protein Interaction Modulation

Docking at Protein-Protein Interfaces

Protein-protein interactions (PPIs) represent promising therapeutic targets but present unique challenges due to their large, flat contact surfaces and limited structural data [67]. Advances in docking protocols have significantly enhanced the field of PPI modulation, with AlphaFold2 (AF2) and MD refinements playing pivotal roles. Recent studies demonstrate that AF2 models perform comparably to experimentally solved structures in docking protocols targeting PPIs, validating their use when experimental data are unavailable [67].

Benchmarking studies have revealed that local docking strategies outperform blind docking, with TankBind_local and Glide providing the best results across different structural types [67]. MD simulations and other ensemble generation algorithms such as AlphaFlow can refine both native and AF2 models, improving docking outcomes, though with significant variability across conformations [67]. This suggests that while structural refinement can enhance docking in some cases, overall performance appears to be constrained by limitations in scoring functions and docking methodologies rather than model quality.

Protocol: MD-Refined Docking for PPI Modulators

Objective: Identify and optimize small molecules that modulate protein-protein interactions.

Step 1 - Target Selection and Characterization

• Identify key residue contacts at the PPI interface from experimental structures or AF2 models
• Characterize interface properties (hydrophobicity, electrostatic potential, cavity size)
• Assess druggability using pocket detection algorithms

Step 2 - Ensemble Generation

• Generate conformational ensemble using 500 ns all-atom MD simulations [67]
• Alternatively, use AlphaFlow sequence-conditioned generative model for ensemble generation
• Cluster trajectories to identify representative conformations

Step 3 - Molecular Docking

• Prepare ligand library focused on PPI-appropriate chemotypes
• Perform docking against multiple ensemble conformations
• Use local docking strategies focused on the interface region
• Employ multiple docking algorithms (TankBind_local, Glide) for consensus

Step 4 - MD Refinement and Binding Analysis

• Run MD simulations (50-100 ns) of top-ranked complexes
• Analyze binding stability through RMSD, contact maps, and interaction fingerprints
• Calculate binding free energies using MM-PB/GBSA or alchemical methods
• Identify key interaction residues and structural water molecules

Step 5 - Experimental Validation

• Select top candidates for biochemical PPI inhibition assays
• Validate binding mode through mutagenesis of key interface residues
• Optimize hits using structure-activity relationship (SAR) guided by simulation insights

Table 3: Key Research Reagent Solutions for CG-MD and Multiscale Modeling

Resource	Type	Function	Example Applications
GROMACS	Software Package	High-performance MD simulation	All-atom and coarse-grained simulations [31] [68]
MARTINI	Coarse-Grained Force Field	Modeling biomolecular interactions	Membrane proteins, lipid bilayers [66] [65]
AWSEM	Coarse-Grained Force Field	Protein folding and dynamics	Protein structure prediction, folding mechanisms [66]
CGSchNet	Machine-Learned CG Model	Transferable coarse-grained simulations	Protein folding, conformational dynamics [66]
AlphaFold2	Structure Prediction	Protein 3D structure prediction	Generating starting structures for docking [67]
CHARMM-GUI	Web-Based Tool	Biomolecular system preparation	Membrane protein insertion, solution building [31]
Bennett Acceptance Ratio	Algorithm	Binding free energy calculation	Ligand binding affinity prediction [31]

Advanced Applications and Future Directions

Cellular and Subcellular Scale Modeling

The ongoing surge of high-resolution structural data is enabling the development of realistic molecular models of cells and subcellular organelles. New tools such as CellPack, which interoperates with CellView and LipidWrapper, can model complex biomolecular systems at the mesoscale, reaching atomic resolution and spatial ranges up to the micron dimension [64]. Combining different datasets across scales of resolution enables direct multiscaling from the standpoint of data integration, facilitating the creation of detailed models of cellular environments that incorporate structural information from cryo-electron tomography, soft X-ray tomography, and other emerging experimental techniques [64].

Future Outlook

The field of coarse-grained and multiscale modeling is rapidly evolving, driven by advances in several key areas. The convergence of simulation and data science is creating new opportunities through improved data sharing, integration, and analytics [64]. The development of increasingly accurate integrative models as initial frameworks for multiscale simulation presents a powerful emerging paradigm for drug discovery [64]. Additionally, the expanding computing landscape, including GPU acceleration, cloud computing, and the emerging horizon of exascale computing, is extending the scope and range of possible simulations [64] [68].

Machine-learned coarse-grained models represent a particularly promising direction, potentially enabling universal, quantitative, and predictive simulations at computational costs several orders of magnitude lower than all-atom methods [66]. As these models continue to improve in accuracy and transferability, they may eventually achieve the long-sought goal of combining the computational efficiency of coarse-grained approaches with the predictive power of all-atom simulations, ultimately providing researchers with a comprehensive toolkit for addressing biological complexity across all relevant scales.

Molecular dynamics (MD) simulations are a cornerstone of modern computational chemistry and drug design, providing atomic-level insights into biomolecular interactions. However, a significant limitation of conventional MD is its treatment of protonation states—the protonation states of titratable amino acid residues and small molecule ligands remain fixed throughout the simulation [69]. This simplification overlooks a crucial aspect of biological reality: proteins and ligands exist in environments with varying pH, and their electrostatic properties dynamically respond to these conditions. Constant pH Molecular Dynamics (CpHMD) has emerged as an advanced technique that overcomes this limitation by incorporating pH as an external thermodynamic parameter, allowing protonation states to fluctuate in response to both conformational dynamics and environmental pH [70] [71].

The importance of accurately modeling protonation states cannot be overstated, particularly in the context of binding affinity refinement research. Experimental studies reveal that in approximately 60% of protein-small molecule complexes, at least one titratable residue in the protein undergoes a change in protonation state between the free and bound forms [71]. Furthermore, an estimated 60-80% of orally administered drugs are weak acids or bases whose protonation states are influenced by cellular pH and the electrostatic environment of their protein targets [71]. This widespread prevalence of protonation-coupled binding necessitates computational approaches that can accurately capture these complex electrostatic dependencies.

The theoretical foundation for modeling pH-dependent binding lies in the binding polynomial formalism developed by Wyman and Tanford [71]. This framework describes the relationship between the apparent binding constant ((K{app})) and pH through the following equation: [ \Delta G^\circ(pH) = \Delta G^\circ{ref} - RT\ln\left(1 + 10^{pKa^C - pH}\right) + RT\ln\left(1 + 10^{pKa^F - pH}\right) ] where (\Delta G^\circ{ref}) represents the reference binding free energy at a specific protonation state, while (pKa^F) and (pK_a^C) denote the pKa values of the ligand in its free and bound states, respectively [71]. This formalism enables researchers to compute binding free energies across a range of pH values, provided the pKa shifts upon binding are known. The CpHMD methodology directly provides these critical pKa values, making it an indispensable tool for predicting pH-dependent binding affinities.

Table 1: Key Terminology in CpHMD and pH-Dependent Binding

Term	Description	Theoretical Importance
Titratable Sites	Functional groups that can gain or lose protons (e.g., Asp, Glu, His, ligand amine groups)	Their protonation states dictate electrostatic interactions and binding energetics [70] [71].
pKa Shift ((\Delta pK_a))	Change in pKa value of a titratable group between free and bound states	A direct measure of the change in electrostatic environment; crucial for accurate (\Delta G) calculation [71].
Binding Polynomial	A statistical mechanical framework that sums over all possible protonation states of a system	Enables calculation of the pH-dependent apparent binding constant, (K_{app}) [71].
Proton-Linked Binding	Binding events that involve a net uptake or release of protons	Results in binding affinities that are highly dependent on the pH of the environment [71].

Key Applications and Case Studies in Drug Discovery

Application to BACE1 Inhibitors for Alzheimer's Disease

BACE1 (β-secretase 1) is a prime therapeutic target for Alzheimer's disease, but its inhibition is complicated by a strong pH dependence. The enzyme functions within acidic cellular compartments (endosomes, pH ~3.5–5.5), and its active site contains a catalytic dyad (Asp32 and Asp228) that must be in a specific protonation state for activity [70]. CpHMD simulations have been instrumental in elucidating why structurally similar BACE1 inhibitors exhibit dramatically different binding affinities within the enzyme's active pH range.

A compelling case study involves two heterocyclic inhibitors with nearly identical binding modes but significantly different inhibitory activities (IC₅₀ of inhibitor 2 is 30-fold higher than inhibitor 1) [70]. Conventional analysis could not explain this discrepancy, as the inhibitor with the higher basic amine pKa (8.7 for inhibitor 2 vs. 7.5 for inhibitor 1) was expected to form a stronger salt bridge with the catalytic aspartates. CpHMD simulations revealed the critical mechanism: binding of inhibitor 1 significantly lowered the pKa of Asp32 from 4.1 (apo state) to 2.5, ensuring it was deprotonated and negatively charged in the active pH range. In contrast, inhibitor 2 produced a negligible pKa shift (Asp32 pKa = 3.7), leaving it partially protonated and weakening key hydrogen bonds [70]. This insight demonstrates how CpHMD can reveal subtle electrostatic mechanisms that dictate inhibitor efficacy.

Host-Guest Systems as Models for pH-Dependent Binding

The cucurbit[7]uril (CB[7]) host-guest system provides a well-characterized model for studying fundamental principles of proton-linked binding. Studies on benzimidazole derivatives have shown that encapsulation within the CB[7] host can induce massive pKa shifts of up to 4.0 pK units (e.g., for Thiabendazole, pKa shifts from 4.6 to 8.6) [71]. This means a guest molecule that is largely deprotonated at neutral pH when free in solution can become predominantly protonated upon binding.

CpHMD simulations have successfully reproduced these large pKa shifts [71]. The subsequent application of the binding polynomial formalism, using these CpHMD-derived pKa values, allows for the accurate calculation of binding free energies across the pH spectrum. This approach has highlighted the potential errors in conventional simulations with fixed protonation states, which can miscalculate binding free energies by >2 kcal/mol at neutral pH for these systems [71]. This error margin is substantial enough to severely mislead lead optimization efforts in drug discovery.

Table 2: Experimentally Measured pKa Shifts for Benzimidazole Guests upon Binding to CB[7] Host [71]

Guest Molecule	Free Guest pKa (pK(_a^F))	Bound Guest pKa (pK(_a^C))	Observed pKa Shift ((\Delta pK))
Benzimidazole (BZ)	5.5	9.0	3.5
Thiabendazole (TBZ)	4.6	8.6	4.0
Fenbendazole (FBZ)	4.8	8.6	3.8
Albendazole (ABZ)	3.5	6.1	2.6
Carbendazim (CBZ)	4.5	7.0	2.5

Detailed Computational Protocols

General Workflow for CpHMD Simulations

The following diagram outlines the standard workflow for setting up and running a CpHMD simulation to study a protein-ligand system, integrating steps from multiple sources [70] [71] [72].

Protocol for CpHMD of Membrane-Protein Drug Complexes

This protocol is adapted from studies investigating the membrane partitioning of ionizable drugs, which utilizes pH replica exchange (pHRE) to enhance sampling [72].

• Step 1: System Setup. For each drug molecule, build a simulation system containing a lipid bilayer (e.g., 128 DMPC lipids). Place the drug in the bulk water phase approximately 15 Å away from the membrane surface to allow for unbiased sampling of the insertion process. Solvate the system with explicit water (e.g., SPC model) and add ions to achieve the desired ionic strength (e.g., 0.1 M) [72].

• Step 2: Parameter and Topology Generation. Use tools like the Automated Topology Builder (ATB) to generate initial topologies for the drug molecules. Manually curate the topologies, particularly for conjugated aromatic rings, to ensure proper 1-4 interaction exclusion according to the force field rules (e.g., GROMOS 54A7). Derive atomic charges for the drugs using quantum mechanical calculations (e.g., Merz-Singh-Kollman analysis based on electrostatic potentials from Gaussian 09) [72]. A critical step is to prepare the drug molecules in both neutral and protonated (cationic) states for subsequent pKa analysis.

• Step 3: Configure pH Replica Exchange CpHMD. Set up multiple replicas (e.g., 4 replicas per system) spanning a molecule-specific pH range (e.g., pH 5.5 to 8.5 in 1.0 unit steps). Configure the simulation to attempt replica exchanges at a fixed frequency (e.g., every 20 ps) to improve sampling efficiency. Run multiple independent replicates (e.g., 5 replicates of 200 ns each) to ensure statistical robustness. For more polar molecules that exhibit slower membrane insertion, extend the simulation time accordingly (e.g., 300-500 ns) [72].

• Step 4: Trajectory Analysis and pKa Calculation. For each conformation sampled, calculate the relative insertion depth of the drug's titratable group relative to the membrane's phosphate atoms. Bin the conformations and their corresponding protonation states by this insertion depth. In each bin, fit the protonation fraction vs. pH data to the Henderson-Hasselbalch equation to compute the pKa as a function of membrane insertion depth. Apply quality control criteria, such as requiring data from multiple replicates and observing monotonic decreases in protonation with increasing pH [72].

Protocol for Calculating pH-Dependent Binding Free Energies

This protocol describes how to use CpHMD results to compute binding free energies across a range of pH values [71].

• Perform Reference Binding Free Energy Calculation. Use a rigorous alchemical method such as Thermodynamic Integration (TI) or Free Energy Perturbation (FEP) to calculate the binding free energy, (\Delta G^\circ_{ref}), for a reference protonation state. This is typically the state where no net proton transfer occurs upon binding (e.g., the deprotonated ligand binding to the receptor) [71].

• Obtain pKa Values via CpHMD. Run separate CpHMD simulations for the free ligand in solution and the protein-ligand complex to determine the pKa values of all relevant titratable groups in both states. These are denoted as (pKa^F) (free) and (pKa^C) (complex) [71].

• Apply the Binding Polynomial. Input the (\Delta G^\circ{ref}), (pKa^F), and (pK_a^C) values into the Wyman binding polynomial equation (provided in Section 1) to calculate the pH-dependent binding free energy profile, (\Delta G^\circ(pH)) [71].

Table 3: Key Software and Computational Tools for CpHMD Simulations

Tool Name	Type/Category	Primary Function in CpHMD Research
GROMACS [2] [72]	MD Simulation Engine	A highly optimized package for performing MD and CpHMD simulations; widely used for its performance and extensive toolset.
CHARMM / AMBER [2]	MD Simulation Engine	Alternative MD engines to which CpHMD methods can be applied; offer different force fields and capabilities.
Automated Topology Builder (ATB) [72]	Parameter Generation	Web server for generating molecular topologies and parameters for a wide range of molecules, compatible with GROMOS force fields.
Gaussian 09 [72]	Quantum Chemistry Software	Used for quantum mechanical geometry optimization and electrostatic potential calculation to derive partial charges for novel ligands.
PETIT Program [72]	Monte Carlo Sampler	Used to sample protonation states based on free-energy terms from Poisson-Boltzmann calculations.
Delphi [72]	Electrostatics Solver	Solves the Poisson-Boltzmann equation to calculate electrostatic free energies, which can be used in pKa predictions and analysis.
Binding Polynomial Formalism [71]	Theoretical Framework	The mathematical framework for integrating CpHMD-derived pKa data into pH-dependent binding free energies.

Critical Insights and Best Practices

Successful implementation of CpHMD for binding affinity studies requires attention to several critical factors. First, protonation state prediction prior to conventional MD is crucial for maintaining structural integrity, especially for membrane proteins involved in proton transport [69]. Simulations initiated with incorrect standard protonation states can rapidly diverge from the experimental structure. Second, researchers must be aware of the systematic deviations in absolute pKa values predicted by some CpHMD methods. The focus should be on the calculated pKa shifts ((\Delta pK_a)), which are often more accurate and are the direct input for the binding polynomial analysis [70]. Finally, for drug-design projects targeting enzymes like BACE1, it is essential to simulate within the enzyme's active pH range, as the protonation states and key interactions (e.g., with the catalytic dyad) are highly pH-dependent [70].

The integration of CpHMD with the binding polynomial formalism provides a powerful and theoretically sound framework for predicting binding affinities across physiological and pathological pH ranges. This approach moves computational chemistry beyond the fixed-protonation-state paradigm, enabling more accurate and predictive simulations that account for the dynamic electrostatic nature of biological systems. As force fields and sampling algorithms continue to improve, CpHMD is poised to become a standard tool in computational drug discovery, particularly for targets where pH-dependent binding is a critical factor.

Leveraging Machine Learning for Force Field Development and Trajectory Analysis

Molecular dynamics (MD) simulations serve as a cornerstone of modern computational materials science and drug discovery, providing atomistic insights into dynamical processes [73]. The accuracy of these simulations is fundamentally governed by the force field (FF), which calculates the potential energy of a system based on atomic coordinates [74]. Traditional force fields, while computationally efficient, often lack quantum-level accuracy, whereas ab initio methods, though accurate, are computationally prohibitive for large systems and long timescales [75]. Machine learning force fields (MLFFs) have emerged as a transformative solution to this accuracy-efficiency trade-off, leveraging data-driven approaches to achieve quantum-chemical accuracy while maintaining computational efficiency for molecular dynamics simulations [75] [73]. This paradigm shift enables the reliable simulation of complex processes such as protein-ligand binding, which is essential for rational drug design [29].

The development of robust MLFFs is particularly impactful for binding affinity refinement in drug discovery. Accurate prediction of protein-ligand binding affinity, which characterizes the strength of biomolecular interactions, is essential for tackling challenges in therapeutic design, protein engineering, and elucidating biological mechanisms [29]. With the U.S. Food and Drug Administration's (FDA) move to phase out animal testing, AI-driven in silico models, including those powered by MLFFs, are poised to become indispensable tools for predictive toxicology and efficacy testing [29]. This application note details advanced protocols and methodologies for developing, validating, and applying MLFFs, with a specific focus on workflows relevant to drug development professionals.

Machine Learning Force Fields: Core Concepts and Workflows

Machine learning force fields aim to address the system-size limitations of accurate ab initio methods by learning the potential energy surface directly from quantum mechanical data [75]. Unlike conventional force fields that use a fixed analytical approximation, MLFFs employ flexible, data-driven models to map atomic configurations to energies and forces [75]. During training and simulation, the local atomic environments are converted into a set of generic descriptors or features. These descriptors are fed into a machine learning algorithm—such as a neural network or kernel regression model—to predict the system's energy. The model parameters are optimized by minimizing the difference between predicted and reference ab initio energies, forces, and stresses [75].

Key Methodological Frameworks

Several advanced frameworks have been developed to enhance the accuracy and data-efficiency of MLFF construction:

• DPmoire for Moiré Systems: This specialized methodology and open-source software package addresses the significant computational challenge of simulating twisted moiré structures, such as transition metal dichalcogenides (TMDs), where lattice relaxation profoundly impacts electronic properties [76]. The standard protocol involves:

  • Constructing 2x2 supercells of non-twisted bilayers with in-plane shifts to generate diverse stacking configurations for training.
  • Performing structural relaxations and Molecular Dynamics (MD) simulations under constraints to collect a robust training dataset.
  • Utilizing a test set composed of large-angle moiré patterns relaxed with ab initio methods to prevent overfitting and ensure transferability [76]. This approach has demonstrated high accuracy, with root mean square errors in forces as low as 0.007 eV/Å for systems like WSe₂ [76].

• Fused Data Learning Strategy: A significant challenge in MLFF development is the inherent inaccuracy of the source data. Density Functional Theory (DFT), a common data source, does not always quantitatively agree with experimental observations [77]. A fused learning strategy concurrently trains an ML potential on both DFT calculations (bottom-up) and experimentally measured properties (top-down), such as mechanical properties and lattice parameters [77]. This hybrid approach results in a molecular model of higher fidelity compared to models trained on a single data source, successfully correcting inaccuracies of DFT functionals for target experimental properties [77].

• End-to-End Differentiable Simulation: This novel framework employs automatic differentiation for both property prediction and force-field optimization. It computes gradients of the simulation output analytically with respect to force field parameters, enabling direct optimization of the force field to reproduce complex target properties like elastic constants, vibrational density of states, and radial distribution functions [74]. This method allows for efficient multi-objective optimization and has been shown to identify optimal parameters for both classical and machine-learned potentials in just a few iterations [74].

The following diagram illustrates a generalized, high-level workflow for developing and applying MLFFs, integrating concepts from these advanced frameworks.

Quantitative Performance Benchmarking

Rigorous benchmarking, such as the TEA Challenge 2023, provides critical insights into the performance of various MLFF architectures across diverse chemical systems [73]. The table below summarizes key quantitative metrics for common MLFF models, highlighting their accuracy in energy and force predictions as well as computational efficiency.

Table 1: Benchmarking MLFF Models from the TEA Challenge 2023 [73]

Model Architecture	Energy MAE (meV/atom)	Force MAE (meV/Å)	Relative Speed (M steps/hr)	Number of Parameters	Key Application Notes
MACE	~1.5	~25	~1.2	2,983,184	High-accuracy, heavier neural network
SO3krates	~2.5	~40	~1.0 (Baseline)	~487,613	Lighter neural network
sGDML	~3.0	~45	~0.7	Not Specified	Suited for molecular systems
SOAP/GAP	~4.0	~55	~2.5	123,000	Kernel-based, strong on periodic materials
FCHL19*	~5.0	~65	~0.5	Not Specified	Kernel-based model

These benchmarks reveal a trade-off between model size, accuracy, and computational speed. While universal MLFFs can achieve broad applicability, their precision may be insufficient for specific applications like moiré systems, where energy scales are on the order of meV. This necessitates the development of MLFFs specifically tailored to individual material systems, where algorithms like NequIP and Allegro can achieve errors of a fraction of a meV per atom [76].

Application Notes & Experimental Protocols

Protocol 1: Building a Specialized MLFF with DPmoire

This protocol outlines the steps for constructing a robust MLFF for complex systems like twisted bilayer materials using the DPmoire package [76].

• Objective: To generate an accurate and transferable MLFF for simulating moiré superlattices.
• Principle: Leverage data from non-twisted structures and targeted ab initio calculations of large-angle moiré patterns to create a model that generalizes to small-angle twisted systems [76].

Step-by-Step Procedure:

• System Preparation (DPmoire.preprocess):

  • Provide the unit cell structures of each layer. DPmoire will automatically combine two layers and generate shifted structures of a 2x2 supercell to create various stacking configurations.
  • The module will also prepare twisted structures for the test set and generate appropriate input files for DFT software (e.g., VASP) based on user-provided templates.

• Data Generation (DPmoire.dft):

  • Perform structural relaxations for each stacked configuration. Constrain the x and y coordinates of a reference atom from each layer to prevent drift.
  • Run Molecular Dynamics (MD) simulations, using an on-the-fly MLFF algorithm (e.g., in VASP) to explore a wider configuration space. Start from a baseline MLFF trained on single-layer structures for stability.
  • Collect data solely from the DFT calculation steps within the MD simulation to ensure data quality.

• Dataset Curation (DPmoire.data):

  • Collect the DFT-calculated data (energy, forces, stress) from all relaxation and MD runs.
  • Merge this data into a comprehensive training set.
  • Compile the data from the relaxed large-angle moiré patterns into a separate test set.

• Model Training (DPmoire.train):

  • Use the DPmoire.train module to configure the training process for a specific MLFF algorithm (e.g., Allegro or NequIP) based on a template.
  • Submit the training job. The resulting trained MLFF can be used in ASE or LAMMPS to perform structural relaxation of target moiré systems [76].

Protocol 2: Fused Data Learning for Experimental Accuracy

This protocol describes how to integrate experimental data with DFT calculations to create an MLFF that corrects for known DFT inaccuracies [77].

• Objective: To train an MLFF that concurrently satisfies DFT-derived and experimentally measured target properties.
• Principle: Use a dual-training approach, alternating between a DFT-data trainer and an experimental-data trainer to optimize the model parameters [77].

Step-by-Step Procedure:

• Initial DFT Pre-training:

  • Train an initial MLFF model using a standard database of DFT calculations (energies, forces, and virial stresses). This serves as a physically reasonable starting point for subsequent optimization.

• Define Experimental Loss Function:

  • Select target experimental properties (e.g., temperature-dependent elastic constants, lattice parameters).
  • Use a method like Differentiable Trajectory Reweighting (DiffTRe) to compute the gradients of these macroscopic properties with respect to the MLFF parameters without backpropagating through the entire MD trajectory [77].

• Alternating Optimization:

  • DFT Trainer Step: For one epoch, optimize the MLFF parameters (θ) to minimize the loss between the model's predictions and the DFT database values (energy, forces, stress).
  • EXP Trainer Step: For one epoch, optimize θ to minimize the loss between simulation-computed observables and the target experimental values, using gradients from DiffTRe.
  • Iterate between the DFT and EXP trainers until convergence. Batch-wise switching can also be employed [77].

• Validation:

  • Validate the final "fused" model on a held-out DFT test set and on experimental properties not used in training (e.g., phonon spectra, properties of different phases) to assess its generalizability [77].

Protocol 3: Trajectory Analysis for Binding Affinity Prediction

This protocol applies MLFF-driven MD simulations to the critical task of predicting protein-ligand binding affinity, a key step in drug discovery [29] [78].

• Objective: To estimate the binding strength and stability of a protein-ligand complex from an MD simulation.
• Principle: Analyze the dynamic trajectory of the complex to compute interaction energies, identify key residues, and estimate binding free energy [78].

Step-by-Step Procedure:

• System Setup and Equilibration:

  • Obtain the 3D structure of the protein-ligand complex, ideally from a reliable database like the PDB.
  • Solvate the complex in an explicit water model (e.g., TIP3P) within a periodic boundary box. Add counterions to neutralize the system.
  • Energy minimize the system and perform equilibration MD runs in the NVT and NPT ensembles at the target temperature (e.g., 310 K) and pressure (1 bar).

• Production MD Simulation:

  • Run a production MD simulation using a validated MLFF (or a classical FF if an accurate MLFF is unavailable) for a sufficient duration (e.g., 50-100 ns) to achieve adequate sampling of the bound state [78]. Save trajectory frames at regular intervals (e.g., every 10-100 ps).

• Trajectory Analysis:

  • Structural Stability: Calculate the Root Mean Square Deviation (RMSD) of the protein backbone and ligand heavy atoms to assess the stability of the simulation.
  • Key Interactions: Analyze hydrogen bonding, hydrophobic contacts, and salt bridges between the ligand and protein residues over time. Identify persistent interactions that stabilize the complex.
  • Binding Free Energy: Employ methods like Molecular Mechanics with Generalized Born and Surface Area solvation (MM/GBSA) on snapshots from the trajectory to estimate the binding free energy (ΔG_bind) [78]. This provides a quantitative measure of binding affinity.

• Validation:

  • Compare the computational predictions with experimentally measured binding constants (e.g., IC₅₀, K_d) where available to validate the protocol [29] [78].

Table 2: Key Software and Databases for MLFF Development and Application

Tool Name	Type	Primary Function	Relevance to MLFF/Binding Affinity
VASP	Software	Ab initio DFT calculator	Generates high-fidelity reference data (energy, forces) for training MLFFs [76] [77].
LAMMPS	Software	Molecular dynamics simulator	A widely used engine for running production MD simulations with trained MLFFs [76].
QuantumATK	Platform	Atomistic modeling platform	Implements ML-FF training workflows (e.g., Moment Tensor Potential) and various simulation engines in a unified platform [75].
JAX-MD	Software	Differentiable MD simulator	Enables end-to-end differentiable simulations for direct force field optimization [74].
Allegro/NequIP	Software	MLFF model architectures	State-of-the-art equivariant neural network models for building accurate MLFFs [76].
PDBbind	Database	Curated protein-ligand complexes and binding affinities	Provides essential data for training and validating binding affinity prediction models [29].
ChEMBL	Database	Bioactive molecules with drug-like properties	Source of bioactivity data for machine learning-guided virtual screening [78].
DPmoire	Software	Automated MLFF construction for moiré systems	Specialized tool for building MLFFs for twisted van der Waals heterostructures [76].

The integration of machine learning into force field development and trajectory analysis represents a paradigm shift in atomic-scale simulation. Methodologies such as specialized training pipelines (e.g., DPmoire), fused data learning, and end-to-end differentiable simulations are pushing the boundaries of accuracy and efficiency [76] [77] [74]. For researchers in drug development, these advances translate directly into more reliable predictions of critical parameters like binding affinity, thereby de-risking and accelerating the early stages of drug discovery [29] [78]. As benchmarked by community challenges, the continuous improvement of MLFF architectures ensures that these tools are not only increasingly accurate but also computationally accessible for simulating biologically and industrially relevant systems [73]. The protocols and resources outlined in this application note provide a foundation for researchers to leverage these powerful techniques in their molecular dynamics workflows for binding affinity refinement.

Benchmarking and Validation: Ensuring Predictive Power in Real-World Scenarios

Accurately predicting the binding affinity between a protein and a ligand is a fundamental challenge in structure-based drug design. The development and validation of computational methods for this task—ranging from classical physics-based approaches to modern machine learning (ML) algorithms—rely heavily on standardized benchmark datasets of protein-ligand complexes with experimentally determined structures and binding affinities. These datasets provide the essential ground truth for training models and fairly evaluating their predictive performance, or "scoring power." Among the most prominent resources are PDBbind, BindingDB, and the Comparative Assessment of Scoring Functions (CASF) benchmark. Understanding their composition, proper use, and inherent limitations is critical for researchers aiming to develop generalizable models, particularly when integrating molecular dynamics (MD) for binding affinity refinement. Recent research highlights that the way these datasets are curated and split significantly impacts a model's ability to predict affinities for novel protein-ligand complexes, an essential requirement for real-world drug discovery applications [79] [80].

Dataset Specifications and Quantitative Comparison

The following table summarizes the core characteristics of the major datasets discussed, highlighting their distinct purposes, scales, and common applications.

Table 1: Key Benchmark Datasets for Protein-Ligand Binding Studies

Dataset Name	Primary Content	Key Metrics & Organization	Primary Use Case	Notable Strengths	Documented Limitations
PDBbind [79] [80]	Curated protein-ligand complexes from the PDB with binding affinity data.	General set (~19,500 complexes), Refined set (higher quality), Core set (used for CASF benchmarking).	Training and benchmarking scoring functions (both classical and ML).	High-quality curation; basis for the standardized CASF benchmark.	Data leakage due to high similarity between training/test splits; structural artifacts in some entries [79] [81].
CASF [79] [81]	A benchmark built upon the PDBbind Core set.	Designed to assess scoring, ranking, docking, and screening powers of scoring functions.	Rigorous, standardized comparison of different scoring functions.	Provides a consistent framework for comparing new methods against established benchmarks.	Inherits potential data quality and leakage issues from the PDBbind Core set [79].
BindingDB [82] [80]	Database of experimental binding affinities (e.g., Ki, Kd, IC50).	Contains over 2.9 million binding data points for 9,300+ targets [82].	Source of binding affinity data for model training and validation.	Vast amount of binding data; covers a wide range of targets and compounds.	Requires pairing with structural data (e.g., from PDB) for structure-based modeling.
BindingNet v2 [82]	Computationally modeled protein-ligand complex structures.	689,796 complexes across 1,794 targets, with high/moderate/low confidence scores.	Augmenting training data for binding pose prediction, especially for novel ligands.	Massive scale and diversity; improves generalization of deep learning models [82].	Contains modeled, not experimentally determined, structures.
LP-PDBBind [79]	A reorganized, "leak-proof" version of PDBbind.	Training, validation, and test splits controlled for protein and ligand similarity.	Training scoring functions for better generalization to novel complexes.	Mitigates data leakage, leading to more realistic performance estimates on new data [79].	A derived dataset, not a primary source.
HiQBind [80]	A high-quality dataset created via the HiQBind-WF workflow.	>18,000 unique PDB entries and >30,000 complexes from BioLiP, Binding MOAD, and BindingDB.	Training and validation requiring high-fidelity, non-covalent complex structures.	Open-source, automated curation corrects common structural errors in PDBbind [80].	A derived dataset, not a primary source.

Experimental Protocols for Dataset Application

Protocol 1: Creating a Leak-Proof Benchmark from PDBbind

Application Note: This protocol is essential for researchers training machine learning scoring functions to ensure evaluated performance reflects true generalization to novel protein-ligand complexes, rather than memorization of similar samples in the training data [79].

• Dataset Acquisition: Download the PDBbind general set (e.g., v2020).
• Data Cleaning (HiQBind-WF Principles):
  • Apply filters to remove covalent protein-ligand complexes, identified via "CONECT" records in PDB files [80] [81].
  • Filter out ligands containing rare elements (beyond H, C, N, O, F, P, S, Cl, Br, I) to avoid learning challenges from infrequent atomic environments [81].
  • Exclude complexes with severe steric clashes (e.g., heavy atom pairs < 2 Å) as they represent non-physical interactions [80] [81].
  • Employ structure-fixing modules (LigandFixer and ProteinFixer) to correct bond orders, protonation states, and add missing protein atoms [80].
• Data Splitting (LP-PDBBind Strategy):
  • Partition the cleaned dataset into training, validation, and test sets.
  • Critical Step: Control for data leakage by ensuring that proteins in different splits have low sequence similarity and ligands in different splits have low chemical structural similarity [79].
  • Avoid random or time-based splits, as they cannot guarantee the absence of highly similar proteins or ligands across splits [79].
• Independent Validation: For a truly independent test, create a benchmark set like BDB2020+ by matching high-quality binding data from BindingDB with PDB structures deposited after a specific cutoff date (e.g., 2020), applying the same similarity control criteria [79].

Protocol 2: Augmenting Training Data with BindingNet v2 for Pose Prediction

Application Note: This protocol is designed for improving the generalization ability of deep learning models, such as Uni-Mol, for predicting the binding poses of novel ligands, where similar crystal structures are unavailable [82].

• Data Integration: Combine the experimentally determined complexes from a dataset like PDBbind with the computationally modeled complexes from BindingNet v2.
• Confidence-Based Sampling: When data volume is a constraint, prioritize complexes from BindingNet v2 based on their confidence scores. Start with high-confidence complexes (hybrid score ≥ 1.2), which show a success rate (ligand RMSD < 2 Å) of 73.79% for the top-ranked pose [82].
• Model Training and Validation:
  • Train the model (e.g., Uni-Mol) on the augmented dataset.
  • Evaluate on a hold-out test set containing novel ligands (e.g., with Tanimoto coefficient < 0.3 to any training ligand) using the PoseBusters benchmark [82].
  • The expected outcome is a significant improvement in success rate. For example, Uni-Mol's success rate for novel ligands increased from 38.55% (trained on PDBbind alone) to 64.25% (augmented with BindingNet v2) [82].
• Physics-Based Refinement (Optional): Further refine the top-ranked poses generated by the DL model using a physics-based method coupled with rescoring (e.g., MM-GB/SA). This combined approach has been shown to increase success rates further, to 74.07%, while also passing PoseBusters validity checks [82].

Protocol 3: Integrating MD Simulations with ML for Affinity Prediction

Application Note: This hybrid protocol uses molecular dynamics simulations to generate structural ensembles, providing 4D dynamic descriptors that can improve binding affinity prediction for flexible targets, especially when starting from docking poses rather than crystal structures [22] [83].

• System Preparation:
  • Start with a protein-ligand complex, ideally from a high-quality curated source.
  • Solvate the complex in an appropriate water model and add ions to neutralize the system.
• MD Simulation and Trajectory Generation:
  • Perform energy minimization and equilibration of the system.
  • Run production MD simulations for a sufficient duration to capture relevant conformational dynamics of the protein and ligand.
• Descriptor Extraction:
  • 2D Descriptors: Calculate traditional molecular fingerprints (e.g., ECFP4) from the ligand's structure [22].
  • 3D Descriptors: Extract energy terms from scoring functions (e.g., Smina, NNscore) based on a single static structure [22].
  • 4D Descriptors: Derive descriptors from the MD trajectory that encode structural dynamic information, such as residue-ligand distance fluctuations, hydrogen bond occupancies, or dynamic contact maps [22].
• Model Training and Evaluation:
  • Train ML models (e.g., Random Forest, XGBoost, SVR) using different combinations of 2D, 3D, and 4D descriptors.
  • Evaluate the model's performance on a held-out test set. Studies show that MD refinement is particularly beneficial for targets with considerable structural flexibility when the initial structure is a docking pose. Its value for rigid targets or when starting from high-resolution crystal structures may be less pronounced [22].

Workflow Visualization

The following diagram illustrates a robust, integrated workflow for binding affinity prediction that leverages curated datasets, deep learning, and molecular dynamics refinement.

Diagram Title: Integrated Workflow for Binding Affinity Prediction

Table 2: Key Software and Data Resources for Binding Affinity Studies

Resource Name	Type	Primary Function in Research	Relevant Context of Use
PDBbind [79] [80]	Dataset	Provides a curated set of protein-ligand complexes with binding affinities for training and benchmarking scoring functions.	The foundational dataset for CASF; often used as a starting point for model development.
CASF [79] [81]	Benchmark	Standardized benchmark for the comparative assessment of scoring functions' scoring, ranking, docking, and screening power.	Used to objectively compare the performance of a newly developed scoring function against existing methods.
BindingDB [82] [80]	Database	Public repository of experimental binding affinities, useful for expanding training data or creating independent test sets.	Sourced for creating the BDB2020+ independent benchmark [79].
HiQBind-WF / PDBBind-Opt [80] [81]	Workflow	An open-source, automated workflow for curating high-quality, non-covalent protein-ligand complex structures from PDB files.	Applied to PDBbind and BioLiP to create improved datasets (HiQBind, BioLiP2-Opt) by fixing structural errors.
LP-PDBBind [79]	Dataset	A reorganized version of PDBbind with minimized data leakage between training and test splits.	Used for training scoring functions to achieve better generalization on novel protein-ligand complexes.
BindingNet v2 [82]	Dataset	A large-scale set of computationally modeled protein-ligand complexes used to augment training data.	Effectively used to improve the generalization ability of the Uni-Mol model for binding pose prediction of novel ligands.
AutoDock Vina [79] [84]	Software	A widely used program for molecular docking, pose generation, and scoring.	Commonly used as a baseline docking tool and classical scoring function in benchmark studies.
MM-GB/SA [82]	Method	A physics-based method for calculating binding free energies, often used for refinement and rescoring.	Coupled with DL model outputs to refine binding poses and improve success rates, as in the BindingNet v2 study [82].

In the context of molecular dynamics (MD) simulations for binding affinity refinement, demonstrating predictive power is not merely a statistical exercise but a critical validation step that bridges computational models with biochemical reality. The correlation between calculated binding affinities and experimentally determined data serves as the primary metric for assessing the practical utility of simulation protocols. Quantitative measures, most commonly the Coefficient of Determination (R²) and Mean Absolute Error (MAE), provide standardized benchmarks for comparing different computational approaches across diverse protein-ligand systems [31] [85]. This protocol outlines standardized procedures for evaluating and validating MD-based binding affinity predictions, ensuring robust assessment of methodological performance in drug discovery pipelines.

Key Metrics for Predictive Power Evaluation

Interpretation Guidelines for R² and MAE

The table below outlines the core metrics used for evaluating predictive power in binding affinity studies:

Table 1: Key Metrics for Evaluating Predictive Power in Binding Affinity Studies

Metric	Definition	Interpretation in Binding Affinity Context	Optimal Range
R² (Coefficient of Determination)	Proportion of variance in experimental data explained by the model [86]	Measures how well computational predictions track with experimental trends; higher values indicate better predictive accuracy	>0.7 (Strong) [85]
MAE (Mean Absolute Error)	Average absolute difference between predicted and experimental values [85]	Quantifies typical error magnitude in affinity units (e.g., pKD, pIC50); lower values indicate higher precision	<0.4 (pIC50/pKD units) [85]
RMSE (Root Mean Square Error)	Square root of the average squared differences [85]	Similar to MAE but penalizes larger errors more heavily; useful for identifying outlier performance	<0.5 (pIC50/pKD units) [85]

Performance Benchmarks from Recent Studies

Recent studies employing advanced MD and machine learning approaches have established performance benchmarks for binding affinity prediction:

Table 2: Performance Benchmarks from Recent Binding Affinity Prediction Studies

Study System	Methodology	R²	MAE	RMSE	Reference
GPCR-ligand complexes	BAR-based binding free energy calculations	0.789	N/R	N/R	[31]
FAK inhibitors (protein-level)	Machine Learning (LightGBM) with combined fingerprints	0.892	0.331	0.467	[85]
FAK inhibitors (cell-level)	Machine Learning (LightGBM) on U87-MG cells	0.789	0.395	0.536	[85]
PLAS-20k dataset	MD/MMPBSA on diverse protein-ligand complexes	Good correlation reported (better than docking)	N/R	N/R	[87]

Figure 1: Workflow for Evaluating Predictive Power of Binding Affinity Predictions

Experimental Protocols for Correlation Validation

MD Simulation Protocol for Binding Affinity Calculation

Application Context: This protocol is adapted from recent studies demonstrating successful correlation with experimental data for membrane protein targets, particularly GPCRs [31].

Materials:

• Protein structure (experimental or predicted)
• Ligand topology files
• MD simulation software (GROMACS, AMBER, or CHARMM)
• Explicit solvent/membrane environment

Procedure:

• System Preparation
  • Generate ligand parameters using GAFF2 force field [87]
  • Prepare protein using AMBER ff14SB or similar force field [87]
  • Solvate system in TIP3P water box with 10Å padding [87]
  • Add ions to neutralize system charge

• Simulation Equilibration

  • Perform energy minimization with backbone restraints (1000 steps) [87]
  • Gradually heat system from 50K to 300K over 200ps [87]
  • Equilibrate in NVT ensemble for 1-2ns [87]
  • Equilibrate in NPT ensemble for 1-2ns [87]

• Production Simulation & Analysis

  • Run production simulation for sufficient duration to capture relevant dynamics (4-100ns) [31] [87]
  • Extract multiple conformational snapshots for binding affinity calculation [87]
  • Calculate binding free energies using BAR, MM/PBSA, or related methods [31]
  • Repeat with multiple independent trajectories for statistical robustness [87]

Machine Learning Enhancement Protocol

Application Context: Integrating machine learning with MD simulations significantly enhances predictive performance for challenging systems with large, flexible ligands [85] [83].

Materials:

• MD simulation trajectories
• Molecular fingerprinting software (PaDEL, RDKit)
• Machine learning frameworks (Scikit-learn, LightGBM, XGBoost)
• Curated experimental binding affinity data

Procedure:

• Feature Generation
  • Compute molecular descriptors (CDK, CDK extended fingerprints) [85]
  • Calculate substructure fingerprint counts [85]
  • Extract dynamic features from MD trajectories [83]

• Model Training & Validation

  • Split data into training/test sets (80:20 ratio) [85]
  • Implement 10-fold cross-validation [85]
  • Optimize hyperparameters using grid search [85]
  • Train multiple algorithms (LightGBM, Random Forest, XGBoost) [85]

• Performance Evaluation

  • Calculate R², MAE, and RMSE on independent test set [85]
  • Compare performance against physical models (MM/PBSA, docking) [87]
  • Validate on external datasets to assess generalizability

Figure 2: Integrated MD/ML Workflow for Enhanced Binding Affinity Prediction

Table 3: Essential Research Reagents and Computational Resources for Binding Affinity Validation

Category	Specific Tool/Resource	Application Context	Key Function
MD Software	GROMACS [31], AMBER [87], OpenMM [87]	Binding free energy calculations	Molecular dynamics simulation engine
Free Energy Methods	BAR (Bennett Acceptance Ratio) [31], MM/PBSA [87], FEP	Binding affinity calculation	Thermodynamic estimation of binding free energies
Machine Learning Libraries	Scikit-learn [85], LightGBM [85], XGBoost [85]	Predictive model development	Implementation of ML algorithms for affinity prediction
Molecular Descriptors	PaDEL [85], RDKit, CDK fingerprints [85]	Feature generation	Transformation of molecular structures to numerical vectors
Validation Datasets	PLAS-20k [87], ChEMBL [85]	Method benchmarking	Curated protein-ligand complexes with experimental affinities
Visualization Tools	Matplotlib, ggplot2 [88]	Correlation analysis	Creation of scatter plots, correlation matrices [86]

Critical Considerations for Robust Correlation Analysis

Statistical Significance and Practical Relevance

When reporting correlation metrics, it is essential to distinguish between statistical significance and practical relevance:

• Assess statistical significance using p-values to determine whether observed correlations are unlikely to occur by chance [89]
• Evaluate practical significance by considering whether the prediction accuracy is sufficient for the specific research context (e.g., lead optimization vs. virtual screening) [86]
• Report confidence intervals for correlation coefficients to communicate estimation precision [86]

Data Quality Considerations

The reliability of correlation analysis depends fundamentally on data quality:

• Experimental data variability: Acknowledge and account for experimental uncertainty in reference measurements [31]
• Structural diversity: Ensure test sets encompass sufficient chemical and structural diversity to assess method generalizability [85]
• Domain applicability: Clearly define the method's applicability domain and avoid extrapolation beyond validated chemical space [83]

Visualization Best Practices

Effective visualization enhances interpretation of correlation results:

• Utilize scatter plots with trend lines to visualize relationships between predicted and experimental values [86]
• Implement correlation matrices with color coding to identify patterns across multiple variables [86]
• Include uncertainty visualization through confidence bands or error bars on correlation plots [88]

• Introduction: Overview of binding affinity prediction and methodological evolution.
• Methodological overview: Classification and characteristics of different scoring approaches.
• Quantitative performance comparison: Tables comparing accuracy, speed, and applicability.
• Experimental protocols: Step-by-step workflows for MD simulations, ML scoring, and conventional methods.
• Research toolkit: Key resources, databases, and software for binding affinity prediction.
• Integration strategies: Hybrid approaches and future directions for method combination.

---

Comparative Analysis of MD Methods vs. Conventional and Machine Learning Scoring Functions

The accurate prediction of protein-ligand binding affinity represents a cornerstone of computational drug discovery, enabling researchers to identify and optimize potential therapeutic compounds with greater efficiency and reduced costs. As the pharmaceutical industry faces increasing pressure to accelerate development timelines while managing expenses, in silico methods have emerged as indispensable tools for prioritizing candidate molecules for experimental validation [29]. The recent phasing out of animal testing by the FDA has further amplified the importance of reliable computational approaches, positioning AI-driven models as transformative elements in modern drug development pipelines [29]. This review examines three fundamental methodological paradigms for binding affinity prediction: molecular dynamics (MD) simulations, conventional scoring functions, and machine learning (ML)-based approaches, with particular emphasis on their comparative strengths, limitations, and optimal application domains.

The evolution of binding affinity prediction methods reflects broader trends in computational chemistry and biology, beginning with physics-based models, expanding to include empirical and knowledge-based statistical approaches, and culminating in today's sophisticated machine learning algorithms [90]. Despite these advances, each method demonstrates characteristic performance profiles, with trade-offs between computational expense, predictive accuracy, and domain applicability that must be carefully considered in research design. Molecular dynamics simulations offer a physics-based foundation for modeling biomolecular interactions but require substantial computational resources [91]. In contrast, conventional scoring functions provide computational efficiency suited for high-throughput screening but may sacrifice accuracy for speed [90]. Machine learning approaches, particularly deep learning models, have demonstrated remarkable predictive performance on benchmark datasets yet face challenges in generalization and interpretability [92].

Within the context of a broader thesis on molecular dynamics simulations for binding affinity refinement, this analysis seeks to establish a comprehensive framework for method selection and integration. By providing detailed experimental protocols, performance comparisons, and practical implementation guidelines, we aim to equip researchers with the knowledge necessary to leverage these complementary approaches effectively in drug discovery campaigns.

Classification of Scoring Approaches

Binding affinity prediction methods can be broadly categorized into three distinct classes based on their underlying theoretical foundations and implementation strategies. Conventional scoring functions represent the historical foundation of computational binding affinity assessment and are further subdivided into physics-based, empirical, and knowledge-based approaches [90]. Physics-based functions utilize molecular mechanics force fields to calculate interaction energies between protein and ligand atoms, typically incorporating van der Waals, electrostatic, hydrogen bonding, and solvation terms [90]. Empirical scoring functions adopt a linear regression approach to fit a combination of energy terms to experimental binding affinity data, while knowledge-based methods employ statistical potentials derived from frequency analyses of atom-atom contacts in known protein-ligand complexes [93].

Machine learning-based scoring functions represent a more recent development in the field, leveraging statistical learning algorithms to infer complex relationships between structural features and binding affinities without relying on predetermined physical equations [93]. These approaches range from traditional machine learning methods like Random Forests and Support Vector Machines to sophisticated deep learning architectures including Graph Neural Networks (GNNs) and Convolutional Neural Networks (CNNs) [92] [93]. The hallmark of ML-based approaches is their data-driven nature, with performance heavily dependent on the quality and quantity of training data [93].

Molecular dynamics (MD) simulations offer a distinct approach based on explicit sampling of the conformational space and energy landscape of protein-ligand complexes [91]. Unlike scoring functions that provide static assessments, MD simulations model the temporal evolution of molecular systems, capturing dynamic processes such as conformational changes, water-mediated interactions, and entropic effects that significantly influence binding affinities [87]. Advanced MD-based methods include alchemical free energy perturbation (FEP) calculations, which provide rigorous binding free energy estimates by simulating non-physical pathways between states, and MM/PBSA (Molecular Mechanics/Poisson-Boltzmann Surface Area) approaches that combine molecular mechanics energy calculations with implicit solvation models [87].

Key Characteristics and Physical Foundations

Each methodological category embodies different theoretical assumptions and physical approximations that fundamentally shape their application profiles. Conventional scoring functions typically rely on simplified representations of molecular interactions, often neglecting or approximating critical effects such as solvation, entropy, and protein flexibility to maintain computational efficiency [93]. While this enables rapid screening of large compound libraries, it can limit accuracy in systems where neglected effects contribute significantly to binding thermodynamics.

Machine learning scoring functions excel at pattern recognition within their training datasets but may struggle with generalization to novel target classes or chemical scaffolds not represented in training data [92]. The black-box nature of many complex ML models can also complicate mechanistic interpretation of predictions, though methods for explainable AI are increasingly addressing this limitation [90]. The performance of ML approaches is intimately tied to the quality and diversity of training data, with models trained on limited or biased datasets potentially learning spurious correlations rather than physically meaningful relationships [29].

Molecular dynamics simulations implement the most physically realistic models among the three approaches, explicitly accounting for solvent effects, full atomic flexibility, and dynamic interactions [91]. By sampling ensembles of configurations, MD methods can capture entropic contributions and conformational rearrangements that occur upon binding, providing insights not only into binding affinities but also into the fundamental mechanisms of molecular recognition [87]. However, this physical fidelity comes with substantial computational costs, limiting throughput and requiring specialized hardware for practical application in drug discovery timelines [91].

Table 1: Fundamental Characteristics of Binding Affinity Prediction Methods

Method Category	Theoretical Basis	Key Advantages	Inherent Limitations
Physics-Based Scoring	Molecular mechanics force fields	Physical interpretability; Force field transferability	Limited accuracy; Implicit solvent approximations
Empirical Scoring	Multivariate linear regression	Computational efficiency; Optimized for affinity prediction	Parameter correlation; Limited physical basis
Knowledge-Based Scoring	Statistical potentials from structural databases	No required parameter fitting; Capture complex interactions	Database dependence; Limited to observed interactions
Machine Learning Scoring	Statistical learning from structural and affinity data	No predetermined functional form; High benchmark performance	Black-box nature; Generalization concerns
Molecular Dynamics	Newtonian mechanics with empirical force fields	Physical realism; Explicit solvation and dynamics	Computational expense; Sampling limitations

Quantitative Performance Comparison

Accuracy Metrics and Benchmarking

Rigorous comparison of binding affinity prediction methods requires standardized benchmarks and evaluation metrics that enable direct performance assessment across diverse protein-ligand systems. The CASF-2016 benchmark has emerged as a widely adopted standard for evaluating scoring functions, measuring performance through several complementary metrics: scoring power (correlation between predicted and experimental affinities), ranking power (ability to rank congeneric ligands by affinity), and docking power (identification of native binding poses) [92]. On this benchmark, modern machine learning scoring functions typically achieve Pearson correlation coefficients (PCC) in the range of 0.85-0.90 with root mean square errors (RMSE) of 1.5-2.0 kcal/mol between predicted and experimental binding affinities [92] [90].

Molecular dynamics-based approaches, particularly free energy perturbation (FEP) methods, have demonstrated performance approaching chemical accuracy (∼1 kcal/mol) for certain protein families and congeneric series [92]. A recent comparative study reported that FEP+ achieved a weighted mean PCC of 0.68 and Kendall's τ of 0.49 on a challenging test set featuring pharmaceutically relevant targets with congeneric ligands [92]. Under similar conditions, the best ML scoring functions achieved PCC values of 0.59 and Kendall's τ of 0.42, representing a meaningful but narrowing performance gap [92]. Importantly, ML methods accomplish these predictions approximately 400,000 times faster than FEP calculations, highlighting the profound trade-off between computational efficiency and accuracy [92].

The integration of MD simulations with machine learning has yielded promising approaches that leverage the strengths of both paradigms. The Dynaformer model, which learns geometric features of protein-ligand interactions from MD trajectories rather than static structures, has demonstrated state-of-the-art performance on the CASF-2016 benchmark, outperforming methods that rely solely on crystallographic structures [91]. This enhancement stems from the model's ability to capture dynamic interaction patterns and conformational ensembles that more accurately represent the thermodynamic determinants of binding affinity [91].

Computational Requirements and Applicability

The practical utility of binding affinity prediction methods depends not only on their accuracy but also on their computational demands, scalability, and applicability to different stages of the drug discovery pipeline. Molecular dynamics simulations are notably resource-intensive, with production-level FEP calculations requiring days to weeks of specialized GPU computing time per compound pair [92]. In contrast, machine learning scoring functions provide affinity predictions in seconds to minutes on standard GPU hardware, enabling virtual screening of extremely large compound libraries [92] [90]. Conventional scoring functions offer the fastest execution, with predictions achievable in seconds or less on standard CPUs, making them suitable for the initial stages of virtual screening where rapid filtering of millions of compounds is necessary [90].

The applicability domains of these methods also vary significantly. Molecular dynamics approaches are particularly valuable during lead optimization phases, where accurate relative binding affinity predictions for congeneric series can guide medicinal chemistry efforts with high confidence [92]. Machine learning scoring functions show strong performance when applied to targets and chemotypes well-represented in training data but may exhibit degraded performance on novel target classes or unusual scaffolds [92]. Conventional scoring functions, while generally less accurate overall, provide reasonable performance across diverse systems without specific training requirements, making them suitable for preliminary screening where compound diversity is high [90].

Table 2: Performance Comparison Across Binding Affinity Prediction Methods

Method	Accuracy (PCC)	Speed	Best Application Context	Key Limitations
Molecular Dynamics (FEP)	0.68-0.80 [92]	Days-weeks per calculation	Lead optimization; Selectivity assessment	Extreme computational cost; Expert parameterization
Machine Learning (AEV-PLIG)	0.59-0.85 [92]	Seconds-minutes per prediction	Virtual screening; Benchmark performance	Limited generalization; Training data dependence
Conventional Scoring	0.40-0.65 [90]	Seconds or less per prediction	Initial screening; Pose prediction	Limited physical model; Systematic errors
MD-ML Hybrid (Dynaformer)	State-of-the-art on CASF-2016 [91]	Minutes per prediction	Affinity prediction when dynamic features matter	Requires MD trajectories; Complex workflow

Experimental Protocols

Molecular Dynamics Simulation Workflow

Molecular dynamics simulations provide a physics-based approach for studying protein-ligand interactions at atomic resolution with explicit solvation. The following protocol outlines the key steps for implementing MD simulations to assess binding affinity, based on established methodologies from the PLAS-20k dataset generation [87]:

System Preparation: Begin with the experimental protein-ligand complex structure from the Protein Data Bank. For missing protein residues, employ loop modeling tools such as those implemented in UCSF Chimera. Add hydrogen atoms appropriate for physiological pH (7.4) using the H++ server or similar tools. Generate ligand topology parameters using the General AMBER Force Field (GAFF2) via the antechamber program, while employing Amber ff14SB for the protein. Solvate the system in an orthorhombic TIP3P water box with a minimum 10 Å extension beyond the protein surface. Add counter ions to neutralize system charge [87].

Energy Minimization and Equilibration: Perform initial minimization using the L-BFGS algorithm with harmonic restraints (10 kcal/mol/Ų) on protein backbone atoms, gradually reducing restraint force over 1000 steps. Conduct additional minimization without restraints. Gradually heat the system from 50K to 300K over 200ps while maintaining backbone restraints. Execute 1ns of simulation in the NVT ensemble followed by 2ns in the NPT ensemble at 300K and 1 atm pressure using a Langevin thermostat and Monte Carlo barostat [87].

Production Simulation and Analysis: Run production simulations for a minimum of 4ns in the NPT ensemble, saving coordinates every 100ps. For binding affinity calculation, employ the MMPBSA method using the following equation:

where ΔEMM represents the molecular mechanics interaction energy (sum of electrostatic and van der Waals components) and ΔGsolv represents the solvation free energy (sum of polar and non-polar contributions) [87]. For higher accuracy, consider alchemical free energy perturbation methods, which require more extensive sampling but provide superior predictive accuracy for congeneric series [92].

Machine Learning Scoring Function Protocol

Machine learning scoring functions offer a data-driven alternative to physics-based methods, leveraging patterns in structural data to predict binding affinities. The following protocol details the implementation of AEV-PLIG, an attention-based graph neural network that represents the state-of-the-art in ML-based affinity prediction [92]:

Data Curation and Preprocessing: Collect protein-ligand complexes with experimental binding affinities from curated databases such as PDBbind (approximately 20,000 structures) or BindingDB. For each complex, extract atomic coordinates and generate graph representations where nodes correspond to atoms and edges represent chemical bonds or intermolecular interactions. Employ atomic environment vectors (AEVs) to describe the local chemical environment of each ligand atom using Gaussian functions of interatomic distances. Implement extended connectivity interaction features (ECIF) to define 22 distinct protein atom types for richer featurization [92].

Model Architecture and Training: Construct a graph neural network using GATv2 layers, an enhanced graph attention mechanism that improves expressiveness over standard GATs. Build the network with five sequential GATv2 layers followed by global pooling and fully connected layers for regression output. Partition data into training, validation, and test sets using time-split or cluster-based splitting to avoid data leakage and ensure realistic performance estimation. Implement appropriate regularization strategies including dropout, weight decay, and early stopping to prevent overfitting. Train the model using the Adam optimizer with mean squared error loss between predicted and experimental binding affinities [92].

Data Augmentation Strategies: To address limited training data, employ augmentation techniques such as template-based ligand alignment and molecular docking to generate additional training examples. For systems lacking experimental structures, utilize homology modeling or docking poses as input representations. Fine-tune the model on target-specific data when sufficient affinity measurements are available to enhance performance on pharmaceutically relevant systems [92].

Conventional Scoring Function Implementation

Conventional scoring functions provide rapid binding affinity estimates and remain widely used for molecular docking and virtual screening. The following protocol outlines the development and application of empirical and knowledge-based scoring functions:

Empirical Scoring Function Development: Identify key interaction terms hypothesized to contribute to binding affinity, including hydrogen bonding, hydrophobic contacts, van der Waals interactions, and loss of ligand conformational entropy. Derive functional forms for each term based on physical principles or statistical observations. Collect a training set of protein-ligand complexes with experimental binding affinities. Use multivariate linear regression to fit weighting coefficients for each term against experimental data. Validate the fitted function on an independent test set not used during parameter optimization [90] [93].

Knowledge-Based Scoring Function Implementation: Compile a database of high-quality protein-ligand complex structures from the PDB. Calculate atom-atom pair distribution functions for protein-ligand contacts within a defined cutoff distance (typically 4.0-6.0 Å). Apply the inverse Boltzmann relationship to convert observed frequencies into effective interaction potentials. Sum pairwise potentials across all protein-ligand atom contacts to generate a comprehensive scoring function. Optionally, refine the scoring function by adjusting parameters against experimental affinity data [93].

Application in Virtual Screening: Prepare protein structures by removing crystallographic waters (except conserved structural waters) and adding hydrogen atoms. Generate multiple conformations for each ligand using tools such as Open Babel or LigPrep. Execute molecular docking with sampling algorithms to generate plausible binding poses. Score each pose using the conventional scoring function and rank compounds by predicted affinity. Select top-ranking compounds for experimental validation [94].

The Scientist's Toolkit

Successful implementation of binding affinity prediction methods requires access to specialized datasets, software tools, and computational resources. The following table catalogues essential resources for researchers in this field:

Table 3: Essential Resources for Binding Affinity Prediction Research

Resource Category	Specific Tools/Databases	Key Features	Application Context
Benchmark Datasets	PDBbind [29], CASF-2016 [92], BindingDB [29], PLAS-20k [87]	Curated complexes with experimental affinities	Method development; Comparative benchmarking
MD Simulation Software	OpenMM [87], GROMACS, AMBER, TorchMD [91]	Molecular dynamics engines with enhanced sampling	Physics-based affinity calculation; Trajectory generation
ML Frameworks	TensorFlow, PyTorch, DeepChem [94]	Deep learning libraries with molecular extensions	Developing custom ML scoring functions
Conventional Docking Tools	AutoDock Vina [94], SMINA, GNINA [94]	Rapid pose generation and scoring	High-throughput virtual screening
Specialized ML Models	AEV-PLIG [92], Dynaformer [91], OnionNet [87]	Pre-trained affinity prediction models	Transfer learning; Baseline comparisons
Free Energy Methods	FEP+ [92], PMX, alchemical-analysis	Alchemical transformation workflows	High-accuracy relative binding affinities

Experimental Design Considerations

When planning binding affinity prediction studies, researchers must consider several critical factors that influence method selection and implementation. For target-specific optimization, prioritize methods that can incorporate relevant physical properties of the specific protein-ligand system under investigation. Membrane proteins, for instance, require explicit lipid bilayer representations in MD simulations, while metalloenzymes demand accurate quantum mechanical treatment of metal coordination chemistry [91].

The choice of evaluation metrics should align with the specific application context. For virtual screening, focus on enrichment factors and early recognition metrics rather than absolute correlation coefficients. For lead optimization, prioritize ranking accuracy within congeneric series and mean absolute error in affinity predictions [92]. When comparing methods, ensure consistent training/test splits and implement proper statistical testing to establish significant performance differences.

Consider hybrid approaches that leverage the complementary strengths of multiple methods. For example, use conventional scoring functions for initial library screening, machine learning scoring for intermediate prioritization, and MD-based methods for final candidate validation [94]. This tiered strategy optimizes the trade-off between computational efficiency and predictive accuracy across different stages of the drug discovery pipeline.

Integration and Future Directions

The convergence of molecular dynamics simulations, machine learning, and conventional scoring functions represents a promising frontier in binding affinity prediction. Hybrid methodologies that integrate physical models with data-driven approaches are increasingly demonstrating superior performance compared to any single method alone [91]. For instance, incorporating MD-derived features such as ligand stability metrics and interaction persistence into machine learning models has been shown to significantly enhance prediction accuracy [91]. Similarly, using machine learning potentials to accelerate MD simulations addresses fundamental limitations in both approaches, enabling more extensive conformational sampling without sacrificing physical accuracy [94].

Future methodological advances will likely focus on several key areas. Explainable AI approaches for interpreting machine learning scoring functions will enhance their utility in medicinal chemistry optimization by providing mechanistic insights into predicted affinities [90]. Multi-scale modeling frameworks that combine quantum mechanical, classical MD, and coarse-grained representations will extend accurate affinity prediction to challenging systems such as metalloenzymes and membrane proteins [91]. The development of large-scale MD datasets like PLAS-20k will enable more robust training of machine learning models while providing benchmarks for method evaluation [87].

As these computational methods continue to mature, their integration into automated drug discovery platforms will play an increasingly important role in reducing attrition rates and accelerating the development of novel therapeutics. By understanding the complementary strengths and limitations of MD simulations, conventional scoring functions, and machine learning approaches, researchers can develop optimized workflows that maximize predictive accuracy while managing computational resources effectively.

Molecular dynamics (MD) simulations have become an indispensable tool in the modern drug discovery pipeline, providing atomic-level insights into the behavior of proteins and other biomolecules with fine temporal resolution. Their impact is particularly pronounced in the study of challenging target classes such as membrane proteins and viral entry proteins, where experimental structure determination remains difficult. By capturing the dynamic nature of biomolecular interactions, MD simulations facilitate the evaluation of binding energetics and kinetics between ligands and their targets, thereby guiding the selection and optimization of candidate molecules [68].

This application note highlights specific success stories where MD-driven approaches have delivered crucial insights for drug discovery, with a particular focus on the refinement of binding affinity predictions. We present quantitative data on the performance of these methods, detailed experimental protocols for their implementation, and essential resource information to empower researchers in adopting these powerful computational techniques.

Success Stories and Key Data

Targeting the SARS-CoV-2 Spike Protein Transmembrane Domain

The SARS-CoV-2 spike protein is a class I fusion protein that mediates viral entry into host cells. While much attention has focused on its receptor-binding domain, the transmembrane domain (TMD) is increasingly recognized as a critical functional region. A 2025 study investigated the functional role of the S protein TMD during viral entry, introducing a series of mutations and insertions within its hydrophobic core [95].

Key findings demonstrated that the sequence and structural integrity of the TMD are crucial for viral entry. Determinants critical for viral entry are distributed throughout the TMD, with a more pronounced contribution from its N-terminal region. The TMD mediates homo-oligomerization through a motif enriched in small residues such as glycine and alanine, underscoring its functional importance beyond mere membrane anchoring [95]. This research paved the way for targeting this domain with small molecules.

In a complementary 2025 study, researchers identified compound 261, a thiazolidinedione that inhibits SARS-CoV-2 entry by binding to the aromatic-rich region of the spike protein directly adjacent to the TMD. Using NMR and MD simulations, the binding pocket was mapped to the TMD, revealing interactions with key conserved residues [96]. This inhibitor achieved an IC₅₀ of 0.3 µM in a tissue model and exhibited potential for pan-coronavirus activity due to the high conservation of the target site among human-infecting coronaviruses [96].

Table 1: Efficacy of SARS-CoV-2 TMD-Targeted Compound 261

Assay Type	Result	Significance
In vitro antiviral assay	IC₅₀ = 0.3 µM	High potency in tissue model
Binding Site Mapping	Aromatic-rich juxtamembrane region	Determined via NMR and MD simulations
Conservation Analysis	High among human coronaviruses	Suggests potential for broad-spectrum activity

MD-Based Binding Affinity Prediction with the PLAS-20k Dataset

Large-scale MD simulations have enabled the creation of extensive datasets for training and validating machine learning (ML) models for binding affinity prediction. The PLAS-20k dataset represents a significant advancement, containing 97,500 independent MD simulations on 19,500 different protein-ligand complexes [87].

The binding affinities (ΔG) calculated from these MD simulations using the Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) method showed a good correlation with experimental values and performed better than traditional docking scores. This held true even for a subset of ligands following Lipinski's Rule of Five and across diverse clusters of complex structures, highlighting the dataset's utility for developing predictive ML models [87].

When the OnionNet model was retrained on the PLAS-20k dataset, it established a strong baseline for predicting binding affinities, demonstrating how such resources can form a new synergy between MD simulations and ML to accelerate drug discovery [87].

Table 2: Performance Overview of MD vs. Docking for Affinity Prediction

Method	Correlation with Experiment	Advantages	Limitations
MD/MMPBSA (PLAS-20k)	Good correlation, better than docking	Accounts for dynamics and solvation; Better at classifying strong/weak binders	Computationally expensive
Molecular Docking	Poor correlation with experimental affinity	Very fast screening of large libraries	Limited conformational sampling; Approximate scoring functions

Experimental Protocols

Protocol: MD-Based Refinement of Membrane Protein Structures

This protocol, adapted from a study on membrane protein structure refinement, uses physics-based techniques with molecular dynamics simulations [97].

• Initial Structure Generation: Generate homology models of the target membrane protein using tools like MODELLER, based on templates identified via PSI-BLAST.
• System Setup:
  • Place the initial model in a solvation environment. For membrane proteins, this entails embedding the protein in an explicit lipid bilayer (e.g., DMPC, DPPC, DLPC) solvated with water molecules. A control simulation in aqueous solvent alone can be performed for comparison.
  • Add counterions to neutralize the system's charge.
• Molecular Dynamics Sampling:
  • Perform energy minimization to remove steric clashes.
  • Carry out equilibration runs with positional restraints on the protein backbone, gradually releasing them.
  • Run extensive, unrestrained production MD simulations (typically on the nanosecond to microsecond scale) using a package like GROMACS, AMBER, or NAMD.
• Trajectory Analysis and Scoring:
  • Extract snapshots from the production trajectory at regular intervals.
  • Score these snapshots using either knowledge-based statistical potentials (e.g., DFIRE, RWplus) or implicit membrane-based scoring functions (e.g., HDGB variants).
• Model Selection and Averaging:
  • Select the top-ranked snapshots based on the chosen scoring function.
  • Average the selected structures to generate a refined model.
• Local Refinement: Further refine the local stereochemistry of the averaged structure using a method like locPREFMD (local Protein structure REFinement via Molecular Dynamics) [97].

Protocol: MD-Driven Identification of Viral Entry Inhibitors

This protocol outlines the process of identifying and characterizing a small-molecule inhibitor targeting the TMD of the SARS-CoV-2 spike protein, based on a published study [96].

• Compound Screening:
  • Screen a compound library using an in vitro antiviral assay. For SARS-CoV-2, this involves incubating the virus with compounds at various concentrations, infecting VeroE6 cells, and quantifying infection (e.g., via plaque assay) after 24 hours.
• Hit Validation: Identify hits, such as compound 261, that show significant, dose-dependent inhibition of viral infection.
• Binding Site Mapping via NMR:
  • Use Nuclear Magnetic Resonance (NMR) spectroscopy to map the inhibitor's binding site on the target protein. This involves observing chemical shift perturbations in the presence of the compound to identify interacting residues.
• Molecular Dynamics Simulations:
  • System Setup: Build a simulation system containing the viral protein (e.g., the spike protein TMD trimer) embedded in a lipid bilayer, solvated in water, with the bound inhibitor.
  • Production Simulation: Run an all-atom MD simulation (e.g., using OpenMM or GROMACS) for tens to hundreds of nanoseconds to observe the stability of the binding interaction and the compound's effect on protein dynamics.
  • Analysis: Analyze the trajectory to characterize the binding mode, identify key protein-inhibitor interactions (hydrogen bonds, hydrophobic contacts), and calculate binding energies [96].
• Mechanism of Action Elucidation: Correlate simulation findings with experimental results to propose a mechanism of action, such as the inhibitor disrupting TMD oligomerization or interfering with conformational changes required for membrane fusion.

Diagram 1: Workflow for MD-driven identification of viral entry inhibitors. The process integrates experimental screening and validation with computational modeling and analysis.

The Scientist's Toolkit

Table 3: Essential Research Reagents and Computational Tools

Item/Tool Name	Function/Application	Relevant Context
GROMACS/AMBER/NAMD	Software for performing molecular dynamics simulations.	Core engine for sampling conformational dynamics and calculating interactions [68].
PLAS-20k Dataset	A curated dataset of protein-ligand affinities from MD simulations.	Training and benchmarking machine learning models for binding affinity prediction [87].
MMPBSA/MMGBSA	End-state method to calculate binding free energies from MD trajectories.	Post-processing simulation data to estimate binding affinity [87].
GAFF2 Force Field	General Amber Force Field for small molecules.	Provides parameters for simulating drug-like ligands in MD simulations [87].
AMBER ff14SB	Force field for simulating proteins.	Used for modeling protein dynamics and interactions [87].
VeroE6 Cells	A mammalian cell line permissive to SARS-CoV-2 infection.	In vitro validation of antiviral activity for viral entry inhibitors [96].
Explicit Lipid Bilayers	Pre-equilibrated membranes (e.g., DMPC, DPPC) for simulation setup.	Essential for creating a realistic environment for MD simulations of membrane proteins [97].
Homology Modeling (MODELLER)	Software for building protein models from related structures.	Generating initial structures for targets without experimental structures [97].

Diagram 2: Logical relationship between key computational tools and data objects in an MD workflow for binding affinity prediction. The force fields define the physical model, the MD software executes the simulation, and analysis tools extract meaningful information from the resulting trajectory.

Molecular dynamics (MD) simulations have become an indispensable tool in computational chemistry and biophysics, particularly for binding affinity refinement in drug discovery. This review critically evaluates the performance of various MD approaches, including alchemical free energy methods, endpoint calculations, and emerging machine learning-integrated techniques. By comparing quantitative performance data across different targets and simulation protocols, we provide a structured analysis of the computational efficiency, accuracy, and applicability of each method. The findings demonstrate that while alchemical methods like Free Energy Perturbation (FEP) and Bennett Acceptance Ratio (BAR) offer superior correlation with experimental binding affinities, their performance is highly dependent on sufficient sampling and specialized hardware acceleration. This analysis provides researchers with clear guidelines for selecting appropriate MD strategies based on their specific project requirements, system characteristics, and available computational resources.

Molecular dynamics simulations have transformed drug discovery by enabling researchers to study protein-ligand interactions and predict binding affinities with unprecedented atomic-level detail. In the context of binding affinity refinement, MD approaches help bridge the gap between static structural information and dynamic binding processes, offering insights that are often challenging to obtain through experimental methods alone [31] [21]. The growing importance of MD is reflected in the expanding molecular modeling market, projected to grow from $8.25 billion in 2024 to $9.44 billion in 2025, driven largely by pharmaceutical development activities [98].

The fundamental challenge in binding affinity prediction lies in accurately capturing the complex thermodynamic and kinetic properties of biomolecular interactions while maintaining computational feasibility. Traditional MD approaches often struggle with sufficient sampling of conformational space and overcoming high energy barriers between states [31]. Recent emphasis has been placed on developing more efficient sampling techniques and re-engineering existing algorithms to better validate experimental results [31]. Additionally, the integration of machine learning with physical force fields has opened new avenues for accelerating simulations while maintaining accuracy [99] [100].

This review systematically examines the strengths and limitations of predominant MD methodologies used for binding affinity refinement, focusing on their theoretical foundations, performance metrics, and practical implementation requirements. By providing structured comparisons and detailed protocols, we aim to equip researchers with the necessary knowledge to select and optimize MD approaches for their specific drug discovery applications.

Performance Comparison of MD Approaches

Quantitative Comparison of Method Performance

The performance of different MD approaches varies significantly based on the target system, sampling efficiency, and implementation details. The table below summarizes key performance metrics for major MD methods used in binding affinity prediction:

Table 1: Performance Comparison of Different MD Approaches for Binding Affinity Prediction

Method	Theoretical Basis	Reported R² with Experimental Data	Computational Cost	Key Strengths	Primary Limitations
BAR	Alchemical perturbation with optimal estimator	0.789 (GPCR targets) [31]	High	High accuracy for membrane proteins; Better correlation with experimental data [31]	Requires multiple intermediate states (lambda values); High computational cost [31]
FEP	Alchemical transformation between states	0.65-0.85 (various targets) [31]	High	Direct calculation of free energy differences; Well-established methodology	Slow convergence with large perturbations; Sensitive to overlap between states
MM-PBSA/GBSA	Molecular Mechanics with implicit solvation	0.4-0.7 (highly system-dependent) [31]	Moderate	Faster than alchemical methods; Applicable to standard MD trajectories	Simplified solvation model; Lower accuracy than alchemical methods [31]
LIE	Linear Response Theory with endpoint sampling	0.5-0.7 (empirically parameterized) [31]	Low	Rapid screening possible; Minimal sampling requirements	Requires system-specific parameterization; Limited transferability
Machine Learning Potentials	Neural network-trained force fields	Emerging data (increasing usage) [99] [101]	Variable (high training, lower inference)	Quantum mechanical accuracy; Transferability across systems	Extensive training data requirements; Black box nature

Critical Analysis of Strengths and Weaknesses

Alchemical Methods (FEP, BAR, TI) Alchemical free energy calculations represent the gold standard for binding affinity prediction due to their rigorous theoretical foundation and high accuracy. The Bennett Acceptance Ratio (BAR) method, in particular, has demonstrated exceptional performance for membrane protein targets like GPCRs, with correlation coefficients of R² = 0.7893 with experimental pKD values [31]. A key strength of alchemical methods is their ability to provide direct thermodynamic interpretation of binding processes. However, these approaches require extensive sampling of multiple intermediate states (lambda values) to overcome energy barriers between initial and final states, resulting in substantial computational costs [31]. Additionally, they are sensitive to the overlap between conformational ensembles of adjacent states, potentially leading to slow convergence or inaccurate results with insufficient sampling.

Endpoint Methods (MM-PBSA/GBSA, LIE) Endpoint methods offer a computationally efficient alternative by focusing sampling only on the bound and unbound states. MM-PBSA/GBSA (Molecular Mechanics Poisson-Boltzmann/Generalized Born Surface Area) utilizes implicit solvent models to estimate binding free energies from MD trajectories, providing a balance between speed and accuracy [31]. The Linear Interaction Energy (LIE) method further reduces computational demands by employing empirical relationships based on endpoint energies [31]. The primary weakness of endpoint methods lies in their simplified treatment of solvation and entropy, which can limit accuracy, particularly for systems with complex solvation patterns or significant conformational changes upon binding. These methods are best suited for rapid screening or ranking of compounds with similar chemical scaffolds rather than absolute binding affinity prediction.

Emerging Approaches (Machine Learning Potentials) Machine Learning Interatomic Potentials (MLIPs) represent a paradigm shift in MD simulations, combining the accuracy of quantum mechanical calculations with the efficiency of classical force fields [99]. Trained on large datasets derived from high-accuracy quantum chemistry calculations, MLIPs can predict atomic energies and forces with remarkable precision, enabling MD simulations of complex systems that were previously computationally prohibitive [99]. The main challenges include the need for extensive training datasets and the "black box" nature of some models. As these methods continue to evolve, they are expected to play an increasingly important role in binding affinity refinement, particularly for systems where traditional force fields lack accuracy.

Experimental Protocols for Key MD Approaches

BAR Method Implementation for Membrane Proteins

The BAR method has shown particular promise for binding affinity calculations of membrane proteins such as GPCRs. The following protocol outlines the key steps for implementing BAR calculations:

System Preparation

• Initial Structure Retrieval: Obtain protein structures from the Protein Data Bank (PDB) and ligand structures from databases such as PubChem or ChEMBL. For GPCR targets, ensure the inclusion of appropriate membrane bilayers [31].
• Parameterization: Assign force field parameters (e.g., AMBER, CHARMM) to both protein and ligand atoms. For small molecules, use tools such as antechamber or CGenFF.
• Membrane Embedding: For membrane proteins, embed the protein in an appropriate lipid bilayer using tools such as CHARMM-GUI or MemGen [31].
• Solvation: Solvate the system with explicit water molecules, adding ions to achieve physiological concentration (e.g., 150 mM NaCl) and neutralize system charge.

Simulation Workflow

• System Equilibration: Perform stepwise equilibration starting with positional restraints on protein and ligand heavy atoms, gradually releasing restraints over 1-2 ns.
• Lambda Window Setup: Define multiple intermediate states (typically 12-24 lambda values) for the alchemical transformation. Use soft-core potentials for Lennard-Jones interactions to avoid singularities.
• Sampling Production Runs: Conduct independent MD simulations at each lambda window with appropriate dynamics parameters (e.g., 2 fs timestep with hydrogen mass repartitioning, constant temperature/pressure) [102].
• Energy Extraction: Calculate potential energy differences between adjacent lambda windows throughout the production trajectories.

Free Energy Analysis

• BAR Implementation: Apply the BAR equation to calculate free energy differences between windows using the collected energy difference data [31]: ΔG = kT * ln(⟨f(U_i - U_j + C)⟩_j / ⟨f(U_j - U_i - C)⟩_i) + C where f(x) = 1 / (1 + exp(x/kT)) and C is a constant iteratively determined to equalize the distributions.
• Error Analysis: Estimate uncertainties using bootstrapping or block averaging methods.
• Convergence Assessment: Verify convergence by examining free energy changes over time and overlap statistics between adjacent lambda windows.

Diagram 1: BAR Method Workflow for Binding Affinity Calculation

MM-PBSA/GBSA Binding Free Energy Protocol

MM-PBSA/GBSA provides a more computationally efficient approach for binding affinity estimation, suitable for high-throughput screening applications:

System Setup and Trajectory Generation

• Complex Preparation: Generate the protein-ligand complex structure, ensuring proper protonation states of ionizable residues.
• MD Simulation: Perform a single conventional MD simulation of the protein-ligand complex (typically 50-200 ns). Multiple shorter replicates may be used to improve sampling.
• Trajectory Processing: Remove rotational and translational motions by aligning trajectories to a reference structure.

Free Energy Calculation

• Trajectory Snapshot Selection: Extract evenly spaced snapshots from the stable portion of the MD trajectory (typically 500-1000 frames).
• Energy Components Calculation: For each snapshot, calculate:
  • Gas-phase molecular mechanics energy (bonded and non-bonded terms)
  • Polar solvation energy using Poisson-Boltzmann (PB) or Generalized Born (GB) models
  • Non-polar solvation energy estimated from solvent-accessible surface area
• Entropy Estimation: Estimate conformational entropy using normal mode analysis or quasi-harmonic approximation (note: this step is computationally demanding and often omitted for screening).
• Binding Free Energy Calculation: Compute the binding free energy for each snapshot and average across all snapshots: ΔG_bind = ⟨G_complex⟩ - ⟨G_protein⟩ - ⟨G_ligand⟩ where each G represents the combined molecular mechanics and solvation energies.

Considerations for Method Selection

• Use MM-PBSA for higher accuracy with increased computational cost
• Use MM-GBSA for faster screening with slightly reduced accuracy
• Consider the implicit membrane GB model (IMGB) for membrane-associated systems
• Account for the limitations of implicit solvation models in capturing specific solvent effects

Hardware and Software Considerations

Optimal Hardware Configurations for Different MD Approaches

The computational demands of MD simulations necessitate careful hardware selection to balance performance and resource utilization. The table below summarizes recommended hardware configurations for different MD methods:

Table 2: Hardware Recommendations for Different MD Approaches

MD Method	Recommended GPU	CPU Requirements	Memory/Storage	Key Considerations
BAR/FEP	NVIDIA RTX 6000 Ada (48GB VRAM) [103]	AMD Threadripper PRO (high clock speed) [103]	64-128GB RAM; Fast SSD storage	Large VRAM essential for membrane protein systems with explicit solvent [31] [103]
MM-PBSA/GBSA	NVIDIA RTX 4090 (24GB VRAM) [103]	Mid-range CPU with 8-16 cores	32-64GB RAM; Moderate storage	Post-processing of trajectories; Multiple replicates possible
Machine Learning Potentials	NVIDIA RTX 6000 Ada or A100 [103]	High core count for training	128GB+ RAM; Large fast storage	VRAM critical for model training; Inference less demanding
Standard MD (Equilibration)	NVIDIA RTX 5000 Ada [103]	Balanced CPU/GPU configuration	64GB RAM; Standard storage	Sufficient for system preparation and equilibration phases

Software Implementation and Performance Optimization

MD Software Packages

• GROMACS: Highly optimized for both CPU and GPU execution; excellent for biomolecular systems with strong scaling on multiple GPUs [31] [102]
• AMBER: Specialized for biomolecular simulations with excellent free energy implementation; optimized PMEMD.CUDA for GPU acceleration [31] [102]
• NAMD: Strong scaling on parallel systems; good support for complex membrane systems [103]

Performance Optimization Strategies

• GPU Offloading: Utilize GPUs for non-bonded force calculations, PME, and bond updating as supported by the MD engine [102].
• Multi-GPU Setups: For large systems (>500,000 atoms), multi-GPU configurations can significantly reduce simulation time, particularly with GROMACS or NAMD [103].
• Integration with AI Platforms: Leverage emerging platforms like AlphaFold2 for initial structure prediction and AI-driven molecular dynamics for enhanced sampling [99] [100].

Diagram 2: Hardware Configuration Decision Framework

Successful implementation of MD approaches for binding affinity refinement requires access to specialized software, hardware, and data resources. The following table details essential components of the MD researcher's toolkit:

Table 3: Essential Research Reagents and Resources for MD Simulations

Resource Category	Specific Tools/Platforms	Primary Function	Key Features
MD Simulation Software	GROMACS [31] [102], AMBER [31] [102], NAMD [103]	Core simulation engines	GPU acceleration; Free energy methods; Biomolecular force fields
Force Fields	CHARMM [31], AMBER [31], OPLS [21]	Define interatomic potentials	Parameterized for proteins, lipids, small molecules; Water models
Structure Preparation	PyMOL [101], CHARMM-GUI [31], RDKit [101]	System setup and visualization	Membrane builder; Parameter generation; Structure analysis
Quantum Chemical Software	Gaussian, ORCA	Reference calculations	Training data for ML potentials; Charge parameterization
Specialized Hardware	NVIDIA RTX 6000 Ada [103], NVIDIA RTX 4090 [103]	GPU acceleration	High VRAM for large systems; CUDA cores for parallel processing
Data Resources	Protein Data Bank [99], PubChem [101], Materials Project [99]	Initial structures and compounds	Experimentally determined structures; Chemical compound libraries
Analysis Tools	MDTraj, VMD, PyTraj	Trajectory analysis	RMSD, RMSF, interaction analysis; Visualization
AI/ML Platforms	AlphaFold2 [99], Schrödigner [100], Exscientia [100]	Structure prediction and design	AI-generated structures; Generative chemistry

This critical examination of different MD approaches for binding affinity refinement reveals a diverse landscape of methods with complementary strengths and limitations. Alchemical methods like BAR and FEP provide the highest accuracy for demanding applications such as GPCR drug discovery but require substantial computational resources and careful implementation. Endpoint methods like MM-PBSA/GBSA offer practical alternatives for high-throughput screening despite their simplified treatment of solvation and entropy. Emerging approaches incorporating machine learning interatomic potentials promise to bridge the gap between accuracy and efficiency, though they remain limited by training data requirements.

The performance of any MD approach is highly dependent on proper system preparation, sufficient sampling, and appropriate hardware selection. Researchers should carefully consider their specific accuracy requirements, system characteristics, and computational resources when selecting an MD strategy. As the field continues to evolve, integration of physical models with machine learning approaches, development of enhanced sampling techniques, and improvements in hardware acceleration will likely address current limitations and expand the applicability of MD simulations in binding affinity refinement.

For future work, we recommend focused development in several key areas: (1) multi-scale methods that combine quantum, classical, and coarse-grained approaches; (2) automated workflows that streamline the process from system setup to analysis; (3) improved force fields that more accurately capture polarization effects; and (4) enhanced sampling algorithms that more efficiently explore conformational space. By addressing these challenges, MD simulations will continue to strengthen their role as an indispensable tool in rational drug design and binding affinity refinement.

Conclusion

Molecular dynamics simulations have matured into an indispensable tool for refining binding affinity predictions, offering unparalleled atomic-level insights into the dynamics of biomolecular interactions. The integration of rigorous alchemical methods, robust end-point approximations, and advanced sampling techniques provides a versatile toolkit for researchers. Looking ahead, the synergy of MD with machine learning and AI, the development of more accurate and scalable force fields, and the rise of AI virtual cells (AIVCs) are poised to further transform the field. These advancements will enhance our ability to model complex biological systems dynamically, accelerating the discovery of novel therapeutics and enabling more personalized biomedical interventions. The future of binding affinity refinement lies in multiscale, integrated computational strategies that closely mirror physiological complexity.

References

• 1. Understanding Computational Molecular Dynamics ... [https://www.ai-futureschool.com/en/physics/understanding-computational-molecular-dynamics-techniques.php]
• 2. Enhancing binding affinity predictions through efficient ... [https://pubs.rsc.org/en/content/articlehtml/2025/sc/d5sc01030f]
• 3. Addressing pathological deficiencies in the simulation of ... [https://www.nature.com/articles/s41598-025-04482-7]
• 4. drMD: Molecular Dynamics for Experimentalists [https://www.sciencedirect.com/science/article/pii/S0022283624005485]
• 5. A Record-Breaking Dataset to Train AI Models has Launched [https://newscenter.lbl.gov/2025/05/14/computational-chemistry-unlocked-a-record-breaking-dataset-to-train-ai-models-has-launched/]
• 6. Exploring Meta's Open Molecules 2025 (OMol25) & Universal ... [https://www.rowansci.com/blog/exploring-open-molecules-2025]
• 7. Boltz-2: Towards Accurate and Efficient Binding Affinity ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12262699/]
• 8. Molecular Dynamics - Statistical Ensembles [https://www.uoxray.uoregon.edu/local/manuals/biosym/discovery/General/Dynamics/Ensembles.html]
• 9. Which ensemble to use for a molecular dynamics simulation? [https://mattermodeling.stackexchange.com/questions/4252/which-ensemble-to-use-for-a-molecular-dynamics-simulation]
• 10. Molecular Dynamics: Ensembles [https://www.compchems.com/molecular-dynamics-ensembles/]
• 11. Canonical (NVT) ensemble — Atomistic Simulation Tutorial [https://docs.matlantis.com/atomistic-simulation-tutorial/en/6_2_md-nvt.html]
• 12. NVT ensemble - VASP Wiki [https://www.vasp.at/wiki/NVT_ensemble]
• 13. NpT ensemble - VASP Wiki [https://www.vasp.at/wiki/NpT_ensemble]
• 14. NPT ensemble — Atomistic Simulation Tutorial [https://docs.matlantis.com/atomistic-simulation-tutorial/en/6_3_md-npt.html]
• 15. NVE ensemble - VASP Wiki [https://vasp.at/wiki/NVE_ensemble]
• 16. Molecular Dynamics simulation [https://docs.matlantis.com/atomistic-simulation-tutorial/en/6_1_md-nve.html]
• 17. A review on description dynamics and conformational ... [https://www.sciencedirect.com/science/article/abs/pii/S0010482524013301]
• 18. - GROMACS 2024.4 documentation Molecular Dynamics [https://manual.gromacs.org/2024.4/reference-manual/algorithms/molecular-dynamics.html]
• 19. Exploring ligand binding pathways on proteins using ... [https://www.nature.com/articles/s41467-021-23157-1]
• 20. Conformational and functional analysis of molecular dynamics ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC3118354/]
• 21. Molecular dynamics simulations of proteins: an in-depth ... [https://pubmed.ncbi.nlm.nih.gov/41003729/]
• 22. Can molecular dynamics simulations improve predictions ... [https://pubmed.ncbi.nlm.nih.gov/36681903/]
• 23. Recent Developments in Free Energy Calculations for Drug ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8387144/]
• 24. Relative Binding Free Energy Calculations Applied to Protein ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC5777225/]
• 25. FEP+ - Schrödinger [https://www.schrodinger.com/platform/products/fep/]
• 26. An optimized thermodynamics integration protocol for ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC10235760/]
• 27. Accurate predictions of protein mutational effects ... [https://www.nature.com/articles/s42004-025-01771-0]
• 28. Bennett acceptance ratio [https://en.wikipedia.org/wiki/Bennett_acceptance_ratio]
• 29. Binding Affinity Prediction: From Conventional to Machine ... [https://arxiv.org/html/2410.00709v2]
• 30. The MM/PBSA and MM/GBSA methods to estimate ligand- ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC4487606/]
• 31. Enhancing binding affinity predictions through efficient ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12118577/]
• 32. Recent Developments and Applications of the MMPBSA ... [https://www.frontiersin.org/journals/molecular-biosciences/articles/10.3389/fmolb.2017.00087/full]
• 33. Enhancing MM/P(G)BSA Methods: Integration of Formulaic ... [https://pubmed.ncbi.nlm.nih.gov/40197754/]
• 34. Assessing the performance of MM/PBSA and MM/GBSA ... [https://pubs.rsc.org/en/content/articlelanding/2024/cp/d3cp04366e]
• 35. Molecular Insights into Structural Dynamics and Binding ... [https://www.mdpi.com/1422-0067/25/24/13482]
• 36. Article Extending MMPBSA for membrane proteins [https://www.sciencedirect.com/science/article/pii/S0006349525003182]
• 37. Advances in G Protein-Coupled Receptor High-throughput ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC7483569/]
• 38. Recent Trends and Applications of Molecular Modeling in ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6073596/]
• 39. Large scale investigation of GPCR molecular dynamics ... [https://www.nature.com/articles/s41467-025-57034-y]
• 40. High-Throughput Fluorescence Polarization Assay to ... [https://pubmed.ncbi.nlm.nih.gov/30925845/]
• 41. TrGPCR: GPCR-Ligand Binding Affinity Prediction Based ... [https://pubmed.ncbi.nlm.nih.gov/37610904/]
• 42. Can molecular dynamics simulations improve the structural ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8186765/]
• 43. Enhancing binding affinity predictions through efficient ... [https://pubs.rsc.org/en/content/articlelanding/2025/sc/d5sc01030f]
• 44. Probing a Model of a GPCR/Ligand Complex in an Explicit ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC1367274/]
• 45. Resolving data bias improves generalization in binding ... [https://www.nature.com/articles/s42256-025-01124-5]
• 46. Accelerated aptamer selection via SELEX and molecular ... [https://pubs.rsc.org/en/content/articlehtml/2025/ra/d5ra06759f]
• 47. Molecular dynamics simulations reveal mechanistic ... [https://link.springer.com/article/10.1007/s10822-025-00633-0]
• 48. Recent Developments and Applications of the MMPBSA Method [https://pmc.ncbi.nlm.nih.gov/articles/PMC5768160/]
• 49. Molecular dynamics simulations of aptamer-binding reveal ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6876308/]
• 50. Multi-pronged molecular insights into flavonoid-mediated inhibition of... [https://pubs.rsc.org/en/content/articlehtml/2025/ra/d4ra09076d]
• 51. Exploring Protein–Ligand Binding-Affinity Prediction | Rowan [https://rowansci.com/blog/exploring-protein-ligand-binding-affinity-prediction]
• 52. Integrating machine learning and molecular dynamics ... [https://www.nature.com/articles/s41598-025-25116-y]
• 53. leveraging AI-agents for customized research workflows [https://pubs.rsc.org/en/content/articlehtml/2025/me/d5me00062a]
• 54. Enabling Scalable AI-Driven Molecular Dynamics ... [https://developer.nvidia.com/blog/enabling-scalable-ai-driven-molecular-dynamics-simulations/]
• 55. A Metadynamics-Based Protocol for the Determination of ... [https://www.mdpi.com/1422-0067/20/8/1970]
• 56. Comparison between metadynamics and umbrella sampling [https://pubmed.ncbi.nlm.nih.gov/26472364/]
• 57. Enhanced Sampling in the Age of Machine Learning [https://arxiv.org/abs/2509.04291]
• 58. Why Are We Focused on Simulation-based Methods for ... [https://www.bio-itworld.com/news/2025/11/14/why-are-we-focused-on-simulation-based-methods-for-predicting-binding-affinity]
• 59. EMFF-2025: a general neural network potential for ... [https://www.nature.com/articles/s41524-025-01809-w]
• 60. An improved reactive force field parameter optimization ... [https://www.sciencedirect.com/science/article/abs/pii/S0927025625001193]
• 61. Evaluating Data Partitioning Strategies for Accurate ... [https://www.sciencedirect.com/science/article/pii/S2001037025004222]
• 62. Automated optimization of force field parameters against ... [https://arxiv.org/html/2402.11169v2]
• 63. Evaluation of Physics‐Based Protein Design Methods for ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12202725/]
• 64. Multiscale Methods in Drug Design Bridge Chemical and ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6445369/]
• 65. Theory and Practice of Coarse-Grained Molecular Dynamics ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8465211/]
• 66. Navigating protein landscapes with a machine-learned ... [https://www.nature.com/articles/s41557-025-01874-0]
• 67. Evaluating ligand docking methods for drugging protein ... [https://jcheminf.biomedcentral.com/articles/10.1186/s13321-025-01067-4]
• 68. Molecular Dynamics Simulations in Drug Discovery and ... [https://www.mdpi.com/2227-9717/9/1/71]
• 69. Role of Protonation States in the Stability of Molecular ... [https://pubmed.ncbi.nlm.nih.gov/38430110/]
• 70. Constant pH Molecular Dynamics Reveals pH-modulated ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC5713896/]
• 71. Protocols Utilizing Constant pH Molecular Dynamics to ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC4306499/]
• 72. 2.5. CpHMD Simulations [https://bio-protocol.org/exchange/minidetail?id=7534506&type=30]
• 73. Crash testing machine learning force fields for molecules ... [https://pubs.rsc.org/en/content/articlehtml/2025/sc/d4sc06529h]
• 74. Force Field Optimization by End-to-end Differentiable ... [https://arxiv.org/html/2409.13844v1]
• 75. What are Machine-Learned Force Fields and How Does It ... [https://www.synopsys.com/glossary/what-are-machine-learned-force-fields.html]
• 76. DPmoire: a tool for constructing accurate machine learning ... [https://www.nature.com/articles/s41524-025-01740-0]
• 77. Accurate machine learning force fields via experimental ... [https://www.nature.com/articles/s41524-024-01251-4]
• 78. Machine Learning and Integrative Structural Dynamics Identify ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12389032/]
• 79. Leak Proof PDBBind: A Reorganized Dataset of Protein-Ligand Complexes for More Generalizable Binding Affinity Prediction [https://arxiv.org/html/2308.09639v2]
• 80. A workflow to create a high-quality protein–ligand binding ... [https://pubs.rsc.org/en/content/articlehtml/2025/dd/d4dd00357h]
• 81. PDBBind Optimization to Create a High-Quality Protein-Ligand Binding Dataset for Binding Affinity Prediction [https://arxiv.org/html/2411.01223v1]
• 82. Augmented BindingNet dataset for enhanced ligand ... [https://www.nature.com/articles/s44386-024-00003-0]
• 83. Predicting receptor-ligand pairing preferences in plant ... [https://www.sciencedirect.com/science/article/pii/S2001037025002454]
• 84. Benchmarking the Structure-Based Virtual Screening ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12363558/]
• 85. Combined machine learning models, docking analysis ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC11490052/]
• 86. 8.2 Correlation analysis and visualization [https://fiveable.me/data-visualization/unit-8/correlation-analysis-visualization/study-guide/DtkSkNxN52HflJB8]
• 87. PLAS-20k: Extended Dataset of Protein-Ligand Affinities ... [https://www.nature.com/articles/s41597-023-02872-y]
• 88. 15 Visualization [https://experimentology.io/015-viz.html]
• 89. Correlation Analysis Guide: Understand Data Relationships [https://innerview.co/blog/correlation-analysis-a-comprehensive-guide-to-understanding-data-relationships]
• 90. Scoring Functions for Protein-Ligand Binding Affinity ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC7613667/]
• 91. Improving Binding Affinity Prediction with Graph‐Based ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC11516055/]
• 92. Narrowing the gap between machine learning scoring ... [https://www.nature.com/articles/s42004-025-01428-y]
• 93. A machine learning approach to predicting protein-ligand ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC3524828/]
• 94. A beginner's approach to deep learning applied to VS and MD ... [https://jcheminf.biomedcentral.com/articles/10.1186/s13321-025-00985-7]
• 95. The sequence and structural integrity of the SARS-CoV-2 ... [https://www.nature.com/articles/s42003-025-08974-0]
• 96. Entry Inhibitors of SARS-CoV-2 Targeting the Transmembrane ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12300635/]
• 97. Structure Refinement of Membrane Proteins via Molecular ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6013386/]
• 98. Rising Drug Development Activities Fuel Growth In The ... [https://www.openpr.com/news/4274973/rising-drug-development-activities-fuel-growth-in-the-molecular]
• 99. Molecular Dynamics Simulations for Materials and ... [https://matlantis.com/en/resources/blog/md-intro/]
• 100. Leading AI-Driven Drug Discovery Platforms [https://www.sciencedirect.com/science/article/abs/pii/S0031699725075118]
• 101. DRUG-TARGET INTERACTION/AFFINITY PREDICTION [https://arxiv.org/html/2502.15346v1]
• 102. Running Molecular Dynamics on Alliance clusters with AMBER [https://computecanada.github.io/molmodsim-amber-md-lesson/26-MD_Software_Benchmarks/index.html]
• 103. Best CPU, GPU, RAM for Molecular Dynamics in 2024 ... [https://bizon-tech.com/blog/best-cpu-gpu-ram-for-molecular-dynamics?srsltid=AfmBOoq7AjpMWo6a7h1LPKvUXflaNjdctBgJarLxU7mRG3C51Ht6atYG]