Molecular Dynamics Simulations of Bat Coronavirus Spike Protein-Receptor Interactions: Implications for Zoonotic Risk Assessment

Introduction

Bat coronaviruses (CoVs) represent a substantial reservoir of genetic diversity with demonstrated capacity for cross-species transmission into domestic animals and, in some cases, humans [1, 2]. The spike glycoprotein (S protein) mediates host cell entry by binding to specific receptors, most commonly angiotensin-converting enzyme 2 (ACE2) for alphacoronaviruses and betacoronaviruses of the subgenus Sarbecovirus [1, 3]. The receptor-binding domain (RBD) of the S protein undergoes conformational rearrangements that modulate receptor engagement and immune evasion [1, 3]. Understanding the biophysical determinants of this interaction at atomic resolution is critical for assessing zoonotic risk [2].

Molecular dynamics (MD) simulations provide a computational framework to probe the time-dependent behavior of the spike protein-receptor complex under solvated, physiological conditions [4, 5]. These simulations capture atomic-level fluctuations, transient binding intermediates, and energetic landscapes that are inaccessible to static structural techniques such as X-ray crystallography or cryo-electron microscopy alone [1, 6]. By integrating MD with free energy perturbation (FEP) and thermodynamic integration (TI) methods, researchers can quantify binding affinities and identify key hotspot residues that govern cross-species compatibility [3, 7]. This article reviews the methodological principles of MD simulations applied to bat coronavirus spike-receptor systems, presents representative case studies from recent bat CoVs, and outlines how these computational approaches inform veterinary zoonotic risk assessments.

MD Simulation Methodology for Spike-Receptor Complexes

System Preparation and Force Fields

A typical MD simulation begins with a high-resolution three-dimensional structure of the spike RBD bound to its cognate receptor, obtained from the Protein Data Bank (PDB) or generated via homology modeling or AlphaFold2 [1, 8]. The complex is embedded in a periodic water box with explicit solvent models (e.g., TIP3P or TIP4P) and neutralized with counterions [4, 9]. The choice of force field is critical: all-atom force fields such as CHARMM36, AMBER ff14SB, or OPLS-AA are commonly employed for protein simulations, with CHARMM36 offering specific parameterization for glycan moieties that are abundant on the spike protein [4, 6, 5]. For membrane-embedded full-length spikes, the protein is inserted into a lipid bilayer (e.g., POPC or a mixed membrane mimetic) using tools such as CHARMM-GUI or MEMBPLUGIN [6, 7]. The system size typically ranges from 200,000 to over 1 million atoms for full spike trimers [6, 7].

Equilibration and Production Runs

Following energy minimization, the system is equilibrated in multiple stages: first with positional restraints on protein heavy atoms, then gradual release over several hundred picoseconds, and finally an unrestrained equilibration under constant temperature (300 K) and pressure (1 bar) using a barostat such as Parrinello-Rahman and a thermostat such as Nosé-Hoover [9, 5]. Production runs for spike-receptor complexes often span 100 ns to several microseconds, with multiple independent replicas to ensure statistical convergence [4, 10]. Specialized hardware (e.g., GPU-accelerated molecular dynamics) and software packages such as GROMACS, NAMD, or AMBER are standard [9, 5]. For rare events such as RBD opening or receptor dissociation, enhanced sampling methods (metadynamics, umbrella sampling, or Markov state models) are essential [10, 3].

Free Energy Calculations

The binding free energy (ΔG_bind) between the bat CoV RBD and a host ACE2 ortholog can be computed using several protocols. The molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) method provides a computationally efficient estimate with moderate accuracy [3, 7]. More rigorous approaches include thermodynamic integration (TI) and free energy perturbation (FEP), which evaluate alchemical transformations of residues at the interface [3, 7]. These methods yield estimates of relative binding affinity changes upon mutation, enabling the identification of residues that are permissive or restrictive for cross-species interactions [3]. The table below summarizes common computational approaches.

Table 1. Overview of free energy methods used in analyzing spike-receptor interactions.

Key Biophysical Parameters in Bat CoV Spike Receptor Recognition

The spike RBD adopts two major conformations in most betacoronaviruses: a "standing-up" state that exposes the receptor-binding motif (RBM) and a "lying-down" state that conceals it [1, 3]. MD simulations have revealed that this conformational equilibrium is sensitive to pH, temperature, and the presence of glycan shields [6, 7]. For bat coronaviruses, the RBM often contains insertions or deletions relative to human-adapted strains, and these structural variations can be analyzed through root-mean-square fluctuation (RMSF) and principal component analysis (PCA) of simulation trajectories [4, 3].

Key residue-residue contacts at the interface include polar, hydrophobic, and electrostatic interactions. For example, a critical lysine residue at position 31 of human ACE2 forms a salt bridge with a glutamate on the RBD of many sarbecoviruses, but this interaction may be absent in bat ACE2 due to a histidine or arginine substitution [1, 3]. MD simulations can quantify the effect of such substitutions on binding free energy [3]. Additionally, hydrogen-bond occupancy analysis across trajectories pinpoints stable versus transient contacts [4, 5].

Case Studies in Bat Coronavirus Spike MD Simulations

SARS-Related Bat Coronaviruses (SARSr-CoV)

Bat CoVs closely related to SARS-CoV-1, such as WIV1 and SHC014, have been studied using MD to assess their potential for human ACE2 utilization [1, 3]. Simulations comparing the RBD of WIV1 with that of SARS-CoV-1 showed that a single mutation at residue 479 (serine to asparagine) substantially enhanced binding to human ACE2, a finding corroborated by MM/PBSA calculations [3]. These studies highlight how MD can identify single-amino-acid determinants of host range expansion [3].

Bat MERS-Related Coronaviruses (MERSr-CoV)

MERS-CoV uses dipeptidyl peptidase 4 (DPP4) as its receptor, and bat relatives of MERS-CoV, such as those found in Neoromicia and Pipistrellus bats, have been structural characterized [1, 2]. MD simulations of the bat MERSr-CoV RBD complexed with human DPP4 indicated that a loop region in the RBM exhibited higher flexibility and reduced binding affinity relative to the camel-adapted virus, suggesting a barrier to direct human spillover [1, 2]. Umbrella sampling provided a potential of mean force that showed a higher dissociation barrier for the bat variant [10, 3].

Novel Bat Coronaviruses (e.g., RaTG13, RmYN02)

The bat coronavirus RaTG13 shares approximately 96% genome identity with SARS-CoV-2 but exhibits differences in the RBD that affect ACE2 binding [1, 3]. MD simulations of RaTG13 RBD with human ACE2 revealed a lower number of inter-residue contacts and a higher RMSD compared to SARS-CoV-2 RBD, correlating with experimental binding data [3]. Free energy decomposition attributed the difference primarily to a lack of key aromatic stacking interactions (e.g., F486 in SARS-CoV-2) [3]. Such computational analyses are now routinely integrated with in vitro binding assays for zoonotic risk scoring [2].

Integration with Zoonotic Risk Assessment

Zoonotic spillover risk is a function of viral prevalence in reservoir hosts, ecological exposure, and molecular compatibility with the recipient host [2]. MD simulations address the molecular compatibility component by providing quantitative biophysical metrics:

• Binding free energy (ΔG_bind) between bat CoV RBD and livestock or companion animal ACE2 orthologs (e.g., canine, feline, swine, bovine).
• Residue interaction networks that identify permissive or restrictive residues.
• Conformational dynamics that influence RBD accessibility and antibody neutralization sensitivity [6, 7].

These metrics can be incorporated into machine learning pipelines that rank emerging viruses by spillover probability [2]. A standardized workflow is depicted in Figure 1.

flowchart TD
    A[Bat Coronavirus Sequence], > B[Homology Modeling / AlphaFold2]
    B, > C[Spike RBD-Receptor Complex Structure]
    C, > D[MD Simulation Setup (Force Field, Solvation, Equilibration)]
    D, > E[Production MD & Enhanced Sampling]
    E, > F[Free Energy Calculation (MM/PBSA, FEP)]
    F, > G[Binding Affinity (ΔG) & Contact Analysis]
    G, > H[Comparison Across Host ACE2 Orthologs]
    H, > I[Zoonotic Risk Score]
    I, > J[Pandemic Preparedness Prioritization]

Figure 1. Computational workflow for evaluating zoonotic risk of bat coronaviruses using molecular dynamics simulations.

Cross-links to related resources on this portal include: Structural Dynamics of Avian Influenza Hemagglutinin: Molecular Modeling and Receptor Binding Predictions for Pandemic Risk Assessment for comparative methodology, and Zoonotic Spillover Pathways and Receptor Binding Evolution in Bat Reservoirs for ecological context. For technical details on simulation protocols, see GROMACS Molecular Dynamics: Setting Up, Simulating, and Analyzing Protein-Water Systems and Molecular Dynamics Simulations of Proteins and Force Fields.

Limitations and Future Directions

MD simulations of bat coronavirus spike-receptor interactions are constrained by force field accuracy, sampling limitations, and the need for experimental validation. Current force fields may not fully capture electronic polarization effects at protein-protein interfaces, though polarizable force fields (e.g., AMOEBA) are under development [4, 5]. Enhanced sampling methods such as replica exchange MD and Markov state models can extend effective simulation timescales but require careful validation of state definitions [10]. Integration with deep learning architectures, such as AlphaFold3 for complex prediction and BindCraft for binder design, offers a promising avenue for high-throughput risk screening [8].

For veterinary applications, extending MD studies to include non-ACE2 receptors used by other coronavirus genera (e.g., APN for alphacoronaviruses, DPP4 for merbecoviruses) is essential to capture the full host range of bat CoVs [1, 2]. Additionally, simulating the spike protein in the context of the full viral membrane environment, including lipid composition effects on fusion, remains a computational challenge [6, 7].

Conclusion

Molecular dynamics simulations provide a powerful biophysical lens through which bat coronavirus spike protein-receptor interactions can be dissected at atomic resolution. By quantifying binding free energies, identifying critical residue contacts, and elucidating conformational dynamics, MD contributes directly to the assessment of cross-species transmission potential. When integrated with ecological surveillance and machine learning frameworks, these computational methods enhance our ability to prioritize emerging bat coronaviruses for veterinary and public health preparedness. Continued methodological advances in force field development, enhanced sampling, and artificial intelligence will further refine the predictive capacity of MD-based zoonotic risk assessment.

References

[1] Structural and Evolutionary Dynamics of Coronavirus Spike Protein: Integrating Cryo-EM, Molecular Dynamics, and Phylogenetic Surveillance. Available at: /knowledge/bioinformatics/structural-evolutionary-dynamics-coronavirus-spike-protein

[2] Zoonotic Spillover Pathways and Receptor Binding Evolution in Bat Reservoirs. Available at: /knowledge/bioinformatics/zoonotic-spillover-pathways-and-receptor-binding-evolution-in-bat-reservoirs

[3] Computational Prediction of Spike Protein Mutations and ACE2 Binding Dynamics in Emerging Coronaviruses. Available at: /knowledge/bioinformatics/computational-prediction-spike-protein-mutations-ace2-binding-dynamics-emerging-coronaviruses

[4] Molecular Dynamics Simulations of Proteins and Force Fields. Available at: /knowledge/bioinformatics/molecular-dynamics-simulations-of-proteins-and-force-fields

[5] Molecular Dynamics Simulations in Biochemistry. Available at: /knowledge/bioinformatics/molecular-dynamics-simulations-in-biochemistry

[6] Molecular Dynamics Simulations of Membrane-Bound Viral Glycoproteins. Available at: /knowledge/bioinformatics/molecular-dynamics-simulations-of-membrane-bound-viral-glycoproteins

[7] Molecular Dynamics Simulations of Viral Envelope Proteins: Insights into Host Recognition and Drug Design. Available at: /knowledge/bioinformatics/molecular-dynamics-simulations-viral-envelope-proteins

[8] Structural Prediction of Viral Envelope Glycoproteins Using AlphaFold2: Implications for Host Receptor Binding and Vaccine Design. Available at: /knowledge/bioinformatics/structural-prediction-viral-envelope-glycoproteins-alphafold2 *** 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.

[9] GROMACS Molecular Dynamics: Setting Up, Simulating, and Analyzing Protein-Water Systems. Available at: /knowledge/bioinformatics/gromacs-molecular-dynamics-simulation-protocols

[10] Markov State Models in Molecular Dynamics Simulations. Available at: /knowledge/bioinformatics/markov-state-models-in-molecular-dynamics-simulations