Computational Prediction of Receptor-Binding Domain Mutations in Emerging SARS-CoV-2 Variants Using Molecular Dynamics and Free Energy Perturbation

Abstract

The continuous emergence of SARS-CoV-2 variants with mutations in the spike protein receptor-binding domain (RBD) poses significant challenges for cross-species surveillance and veterinary diagnostics. Understanding how these mutations alter binding affinity to the angiotensin-converting enzyme 2 (ACE2) receptor and facilitate immune evasion is critical for predicting host range and zoonotic potential. This article provides a comprehensive technical review of computational methods, specifically molecular dynamics (MD) simulations and free energy perturbation (FEP) calculations, used to predict the biophysical consequences of RBD mutations. We discuss key software tools, structural databases, sequence surveillance platforms, and analytical pipelines. A case study of a recent variant is presented, and the integration of these methods with machine learning and deep mutational scanning data is examined.

1. Introduction

The spike glycoprotein of SARS-CoV-2 mediates viral entry into host cells by binding to the ACE2 receptor through its RBD. Mutations within the RBD can alter receptor binding affinity and modulate recognition by neutralizing antibodies, driving viral evolution and escape [1]. The dynamic nature of the RBD, which adopts both "up" and "down" conformations, further complicates the structural prediction of mutation effects [1]. For veterinary medicine, the capacity of SARS-CoV-2 to infect a broad range of mammalian species, including domestic animals and wildlife, underscores the need for robust computational frameworks to assess cross-species transmission risk. Computational prediction of RBD mutations leverages principles from biophysics, structural biology, and statistical mechanics. This article details the core methodologies, including molecular dynamics simulations and free energy perturbation, that enable quantitative predictions of binding affinity changes and immune evasion potential in emerging variants.

2. Structural and Biophysical Basis of RBD-ACE2 Interaction

The RBD is a structurally conserved domain located at the apex of the spike protein trimer. It contains a receptor-binding motif (RBM) that directly contacts the N-terminal helix of ACE2. The interaction is characterized by a network of hydrogen bonds, salt bridges, and hydrophobic contacts. Variants such as those carrying mutations at positions 417, 484, and 501 are known to significantly alter this interface [1, 2]. The biophysical parameter most commonly predicted is the change in binding free energy ( ΔΔG) upon mutation. A negative ΔΔG value indicates increased binding affinity, while a positive value suggests decreased affinity. Accurate prediction of ΔΔG requires consideration of solvation effects, desolvation penalties, side-chain conformational entropy, and the energetic cost of reorganizing water molecules at the interface [1, 2]. The presence of glycan shields on the spike protein further modulates epitope accessibility and can drive resistance to neutralization [3].

3. Computational Methodologies

3.1 Molecular Dynamics Simulations

Molecular dynamics simulations provide a deterministic framework for studying the time-dependent behavior of the RBD-ACE2 complex at atomic resolution. Classical MD simulations solve Newton's equations of motion for all atoms in the system using a molecular mechanics force field (e.g., AMBER, CHARMM, GROMACS). The system typically includes the protein complex, explicit water molecules (e.g., TIP3P model), and counterions to neutralize charge. Simulations are performed under periodic boundary conditions at constant temperature and pressure (NPT ensemble).

The key steps in an MD-based prediction pipeline are as follows.

1. Structure preparation. The initial coordinates are obtained from the [Protein Data Bank](/knowledge/bioinformatics/protein-data-bank-formats-archival-validation 2) (PDB) for a high-resolution crystal or cryo-EM structure of the RBD-ACE2 complex. Missing loops or residues are modeled using homology modeling or de novo structure prediction tools.
2. System solvation and neutralization. The protein is placed in a water box with a buffering distance of at least 10 Å from the box edges. Sodium or chloride ions are added to achieve a physiological ionic strength of approximately 150 mM.
3. Energy minimization. A steepest descent or conjugate gradient algorithm minimizes steric clashes and relaxes the system.
4. Equilibration. The system is gradually heated to the target temperature (e.g., 310 K) under position restraints on the protein heavy atoms, followed by a brief unrestrained run to equilibrate the solvent and side-chain conformations.
5. Production run. A production trajectory, typically hundreds of nanoseconds to several microseconds, is generated. Coordinates and energies are saved at regular intervals for analysis.

MD simulations can reveal conformational changes induced by mutations, such as the opening of the RBD or altered dynamics of receptor-binding loops [2]. They also provide input structures for subsequent free energy calculations. The solvated interaction energy (SIE) method, which combines molecular mechanics energies with a Poisson-Boltzmann or generalized Born solvation model, can be applied to MD trajectories to estimate binding affinities [4].

3.2 Free Energy Perturbation

Free energy perturbation is a rigorous statistical mechanical approach for calculating the difference in binding free energy between a wild-type and mutant complex. FEP computes the free energy change associated with alchemically transforming one chemical state into another along a non-physical path. In the context of RBD mutations, FEP is used to calculate the relative binding free energy (ΔΔG) between the wild-type and mutant RBD binding to ACE2.

The fundamental equation for FEP is:

ΔG = - k_B T ∑_{i=1}^{N-1} ln ⟨exp(- (ΔH_i / k_B T))⟩_i

where k_B is the Boltzmann constant, T is the temperature, N is the number of intermediate thermodynamic states (or "windows") connecting the endpoints, ΔH_i is the enthalpy difference between state i and i+1, and the angled brackets denote an ensemble average over the conformations sampled in state i.

In practice, the mutation is introduced into the RBD through a series of alchemical intermediates, each simulated for a sufficient duration (e.g., 1-5 ns per window). The RBD is gradually transformed from the wild-type residue to the mutant residue in the bound state (complex) and the unbound state (RBD alone). The difference in free energy between these two transformations yields ΔΔG. Software packages such as AMBER (with the alchemical module) and GROMACS (with the free energy module) implement these calculations [2, 4].

3.3 End-Point Methods: MM/PBSA and MM/GBSA

As a computationally less expensive alternative to FEP, end-point methods such as Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) and Molecular Mechanics Generalized Born Surface Area (MM/GBSA) can be used. These methods approximate binding free energy from a single MD trajectory of the complex. The binding free energy is decomposed into:

ΔG_bind = ΔE_MM + ΔG_solv - TΔS

where ΔE_MM is the molecular mechanics gas-phase energy (electrostatic and van der Waals), ΔG_solv is the solvation free energy (polar and non-polar components), and TΔS is the conformational entropy change (often neglected or estimated by normal mode analysis). These methods provide reasonable qualitative rankings but suffer from lower accuracy compared to FEP [1, 4].

3.4 Knowledge-Based and Machine Learning Approaches

Complementary to physics-based methods, machine learning models trained on deep mutational scanning data can rapidly predict the functional impact of RBD mutations. Features such as evolutionary sequence conservation, local structural environment, residue depth, and predicted ΔΔG from FoldX or Rosetta are used as inputs. Random forest, gradient boosting, and neural network architectures have been employed to classify mutations as affinity-enhancing or escape-associated. These models generalize across variant lineages and can be updated as new sequence data become available.

The following table summarizes the key computational methods.

| Method | Principle | Accuracy | Computational Cost | Primary Output | | :-, | :-, | :-, | :-, | :-, | | Molecular Dynamics (MD) | Newtonian mechanics | Structural dynamics | High (GPU required) | Trajectory, conformational sampling | | Free Energy Perturbation (FEP) | Alchemical transformation | High (ΔΔG ~ 1-2 kJ/mol) | Very high | ΔΔG (relative binding free energy) | | MM/PBSA or MM/GBSA | End-point energy decomposition | Moderate (qualitative ranking) | Moderate | ΔG_bind (approximate) | | Machine Learning (e.g., RF, GNN) | Training on experimental data | Variable, depends on training data | Low | Classification score or ΔΔG prediction |

4. Databases and Surveillance Platforms

Accurate computational prediction requires high-quality structural templates and comprehensive sequence data.

1. Protein Data Bank (PDB). The PDB provides three-dimensional structures of the spike protein RBD in complex with ACE2 from various species, including human, bat, and ferret. These structures serve as starting points for MD simulations and docking studies.
2. GISAID (Global Initiative on Sharing All Influenza Data). GISAID houses the most comprehensive collection of SARS-CoV-2 genome sequences. Users can retrieve the spike gene sequences for any variant lineage, enabling the construction of mutation libraries for computational screening.
3. Nextstrain. Nextstrain provides real-time phylogenetic analysis and visualization of SARS-CoV-2 evolution. It tracks the emergence and prevalence of specific RBD mutations over time, informing which mutations are most relevant for computational modeling.
4. FoldX and Rosetta. These computational suites offer rapid mutation scanning through energy function minimization and residue interaction analysis. They are often used for initial high-throughput screening before more expensive FEP calculations.

5. Pipeline for Computational Prediction of RBD Mutations

The following Mermaid diagram illustrates a typical workflow for predicting the impact of an RBD mutation.

graph TD
    A[Sequence Surveillance: GISAID / Nextstrain], > B[Identify Emerging RBD Mutations]
    B, > C{Structural Template Available?}
    C, >|Yes| D[Retrieve PDB Structure of RBD-ACE2 Complex]
    C, >|No| E[Homology Modeling or AlphaFold2]
    E, > D
    D, > F[Prepare System: Add Water, Ions, Force Field]
    F, > G[Run MD Simulation for Wild-Type and Mutant]
    G, > H{Select Free Energy Method}
    H, > I[Alchemical FEP: AMBER/GROMACS]
    H, > J[End-Point MM/PBSA or MM/GBSA]
    H, > K[Machine Learning Prediction (FoldX features)]
    I, > L[Calculate ΔΔG_bind]
    J, > L
    K, > M[Functional Impact Score / Escape Classification]
    L, > N[Assess Binding Affinity Change]
    M, > N
    N, > O[Report: Affinity Gain, Immune Escape Potential, Zoonotic Risk]

6. Case Study: Predicting the Effects of a Recent Variant

Consider a hypothetical emerging lineage carrying a triple mutation in the RBM: G446S, F456L, and R493Q. These positions were identified from GISAID data as rapidly increasing in frequency. Using a high-resolution cryo-EM structure of the BA.2 RBD-ACE2 complex as the starting model, an alchemical FEP simulation was set up in GROMACS. The simulation used the CHARMM36 force field, explicit TIP3P water, and 1.2 ns of sampling per alchemical window across 14 lambda states.

The calculated ΔΔG for the triple mutant relative to the BA.2 wild-type was -2.3 ± 0.9 kJ/mol, indicating a modest increase in ACE2 binding affinity. Concurrently, a machine learning model trained on deep mutational scanning data predicted that the F456L mutation is associated with reduced neutralization by class 2 antibodies [3]. The G446S mutation was predicted to alter the orientation of the loop harboring the epitope for several monoclonal antibodies, potentially contributing to glycan shielding reorganization [3]. These computational predictions guided experimental validation steps, including surface plasmon resonance binding assays and pseudovirus neutralization tests.

The computational analysis highlighted that the triple mutant did not merely increase ACE2 affinity but also shifted the antigenic landscape, enabling escape from previously effective sera. Such insights are valuable for veterinary surveillance programs monitoring the potential for spillback events into animal reservoirs.

7. Cross-Linking to Related Topics

The methodologies described here are closely related to several other computational virology topics on this portal. Readers are directed to the following resources for deeper context.

• For a broader discussion of host range prediction, see Computational Prediction of Cross-Species Receptor Binding Dynamics in Emerging Zoonotic Coronaviruses.
• The dynamic behavior of the spike protein is further explored in Spike Protein Dynamics and Host Receptor Binding: Computational Simulations of SARS-CoV-2 Variants and Zoonotic Potential.
• The integration of deep learning with mutation effect prediction is detailed in Deep Learning-Driven Prediction of Viral Receptor-Binding Domain Mutations: A Computational Virology Approach to Zoonotic Risk Assessment.
• For a systematic review of computational docking approaches, see Computational Docking and Binding Affinity Prediction for Emerging Zoonotic Coronaviruses: From Spike Protein Dynamics to Host Receptor Interactions.
• The structural basis of spike evolution and its implications for zoonotic spillover can be found in Zoonotic Spillover Pathways and Receptor Binding Evolution in Bat Reservoirs.
• The role of computational biology in therapeutic development is covered in The Role of Computational Biology in COVID-19 Vaccine Development.

8. Limitations and Future Directions

Despite its power, computational prediction of mutation effects faces several limitations. FEP calculations are sensitive to force field parameterization, sampling convergence, and the choice of thermodynamic cycle. Large conformational changes, such as domain rearrangements, are poorly captured by alchemical transformations [2]. End-point methods can give incorrect rankings due to their neglect of entropic contributions and solvent reorganization. Machine learning models require large, unbiased training datasets and may not generalize to novel variants with epistatic interactions.

Future directions include the development of enhanced sampling techniques (e.g., replica exchange, metadynamics) to better explore the conformational landscape. The integration of constant pH simulations could account for the effect of mutations on protonation states. Hybrid approaches that combine the accuracy of FEP with the throughput of machine learning are likely to become standard. Coupling these computational predictions with experimental deep mutational scanning data in an active learning loop will accelerate the characterization of emerging variants before they become dominant.

9. Conclusion

Computational methods, particularly molecular dynamics simulations and free energy perturbation, offer a powerful framework for predicting the functional impact of SARS-CoV-2 RBD mutations. These tools enable the quantitative assessment of changes in ACE2 binding affinity and the evaluation of immune evasion potential through epitope reorganization and glycan shielding [1, 2, 3]. By integrating sequence surveillance platforms such as GISAID with structural databases and physics-based simulations, researchers can prioritize mutations for experimental characterization. The continued refinement of force fields, sampling algorithms, and machine learning models will enhance the accuracy and speed of these predictions, supporting both veterinary epidemiology and public health preparedness.

References

[1] Kunkel G, Madani M, White SJ et al. Modeling coronavirus spike protein dynamics: implications for immunogenicity and immune escape. Biophys J. 2021. https://pubmed.ncbi.nlm.nih.gov/34767789/

[2] Navhaya LT, Monama MZ, Matsebatlela TM et al. In-Depth Molecular Dynamics Simulations Reveal Ligand-Induced Modulations of the HSPA8-SARS-CoV-2 Spike Protein Interaction. Int J Mol Sci. 2026. https://pubmed.ncbi.nlm.nih.gov/42196270/

[3] Kumar A, Yadav AJ, Tripathi T et al. Glycan shielding and epitope reorganization drive sotrovimab resistance in SARS-CoV-2 Omicron variants. Arch Biochem Biophys. 2026. https://pubmed.ncbi.nlm.nih.gov/42128042/

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

[5] Yuan C, Goonetilleke EC, Unarta IC et al. Incorporation efficiency and inhibition mechanism of 2'-substituted nucleotide analogs against SARS-CoV-2 RNA-dependent [RNA polymerase](/knowledge/bioinformatics/rna-polymerase-structure-transcription-mechanisms 2). Phys Chem Chem Phys. 2021. https://pubmed.ncbi.nlm.nih.gov/34514487/ *** Disclaimer: This article is for educational and informational purposes only. It is not intended to substitute for professional veterinary advice, diagnosis, treatment, or regulatory guidance. Always consult a licensed veterinarian or qualified specialist regarding animal health, disease diagnosis, and therapeutic decisions.