Computational Prediction of Cross-Species Receptor Binding Dynamics in Emerging Zoonotic Coronaviruses
By Dr. Zubair Khalid, DVM, MS, PhD ·

Key Takeaways
- The zoonotic potential of coronaviruses is primarily determined by the interaction between the viral spike (S) glycoprotein's receptor-binding domain (RBD) and host cell surface receptors, most commonly angiotensin-converting enzyme 2 (ACE2). Computational methods like molecular docking, molecular dynamics (MD) simulations, and binding free energy calculations (e.g., MM-PBSA) are crucial for predicting these cross-species binding dynamics and assessing spillover risk.
- Key structural features of the RBD, particularly the receptor-binding motif (RBM), and specific amino acid residues within it (e.g., positions 486, 493, 498, 501 in SARS-CoV-2) are critical determinants of binding affinity to different host ACE2 orthologs. Mutations at these sites can significantly alter host tropism.
- Machine learning classifiers, trained on sequence and structural features of RBD-ACE2 interactions, are employed to screen large viral sequence databases for variants with high predicted affinity to susceptible animal hosts, thereby guiding targeted surveillance efforts.
- Molecular dynamics simulations provide insights into the dynamic nature of the RBD-ACE2 complex, revealing conformational changes, hydrogen bond stability, and collective motions that influence binding. Analysis of these trajectories, including RMSF and hydrogen bond occupancy, helps elucidate species-specific adaptations.
- An integrated computational pipeline, combining machine learning for initial screening, molecular docking for pose prediction, MD simulations for refinement, and MM-PBSA for free energy estimation, followed by experimental validation (e.g., SPR, pseudovirus assays), is essential for robust cross-species risk assessment.
- Limitations in computational modeling include force field inaccuracies, inadequate representation of solvent effects, and challenges in incorporating glycan flexibility. Post-entry factors beyond receptor binding also play a significant role in determining successful host tropism and spillover.
Introduction
The emergence of zoonotic coronaviruses from animal reservoirs represents a persistent threat to veterinary and public health. Coronaviruses within the genera Alphacoronavirus, Betacoronavirus, Gammacoronavirus, and Deltacoronavirus circulate widely in bats, birds, and domestic livestock [<a href="#ref-1">1</a>]. The capacity for a coronavirus to cross the species barrier is primarily governed by the molecular interaction between the viral spike (S) glycoprotein receptor-binding domain (RBD) and a host cell surface receptor, most commonly angiotensin-converting enzyme 2 (ACE2) [<a href="#ref-2">2</a>]. Understanding and predicting these cross-species receptor binding dynamics is essential for assessing zoonotic spillover risk and for designing surveillance strategies in animal populations [<a href="#ref-3">3</a>].
Computational methods have become indispensable for modeling these interactions at atomic resolution. Molecular docking simulations, molecular dynamics (MD) trajectories, binding free energy calculations, and machine learning classifiers collectively enable the prediction of host tropism from sequence and structural data [<a href="#ref-4">4</a>]. This article provides an exhaustive review of these computational approaches, focusing on their application to RBD-ACE2 interactions across bat, avian, and mammalian species. The discussion emphasizes the biophysical mechanisms underlying host switching and the structural motifs that facilitate adaptation to new receptors.
Structural Basis of Coronavirus Receptor Binding
Coronavirus spike proteins are class I viral fusion glycoproteins that mediate host cell attachment and entry [<a href="#ref-5">5</a>]. The S protein is a homotrimer, with each monomer composed of S1 and S2 subunits. The S1 subunit contains the RBD, which directly engages the host receptor [<a href="#ref-6">6</a>]. For betacoronaviruses such as SARS-CoV, SARS-CoV-2, and related bat SARS-like coronaviruses (SL-CoVs), the primary receptor is ACE2 [<a href="#ref-2">2</a>]. Other coronaviruses use alternative receptors: for example, Middle East respiratory syndrome coronavirus (MERS-CoV) uses dipeptidyl peptidase 4 (DPP4), and some alphacoronaviruses use aminopeptidase N (APN) [<a href="#ref-7">7</a>].
The RBD adopts a beta-sheet-rich core with a receptor-binding motif (RBM) that forms the direct contact interface with ACE2 [<a href="#ref-8">8</a>]. Key structural features of the RBD-ACE2 interface include a network of hydrogen bonds, salt bridges, and hydrophobic contacts [<a href="#ref-9">9</a>]. Critical residues in the RBM, such as those at positions 486, 493, 498, and 501 (SARS-CoV-2 numbering), have been shown to modulate binding affinity across species [<a href="#ref-10">10</a>]. Mutations at these positions can enhance or reduce binding to ACE2 orthologs from different animal hosts [<a href="#ref-11">11</a>].
Molecular Docking Simulations for Cross-Species Binding Prediction
Molecular docking is a computational technique that predicts the preferred orientation of a ligand (the RBD) when bound to a receptor (ACE2) to form a stable complex [<a href="#ref-12">12</a>]. Docking algorithms sample conformational space and score candidate poses using energy-based scoring functions [<a href="#ref-13">13</a>]. For cross-species studies, docking is used to evaluate the binding affinity of RBD variants against ACE2 orthologs from multiple species, including bats, civets, swine, ferrets, and poultry [<a href="#ref-14">14</a>].
The docking workflow typically involves preparing the RBD and ACE2 structures from experimentally determined coordinates (X-ray crystallography or cryo-electron microscopy) or from homology models [<a href="#ref-15">15</a>]. Rigid docking approaches treat both molecules as rigid bodies, while flexible docking allows side-chain or backbone flexibility in the RBM [<a href="#ref-16">16</a>]. Software packages commonly used for this purpose include AutoDock Vina, HADDOCK, and RosettaDock, which employ different scoring functions such as empirical, force field-based, or knowledge-based potentials [<a href="#ref-17">17</a>].
Validation of docking results is performed by comparing predicted poses with known co-crystal structures and by calculating root-mean-square deviation (RMSD) values [<a href="#ref-18">18</a>]. A successful docking simulation reproduces the native binding mode with an RMSD below 2.0 Angstroms [<a href="#ref-19">19</a>]. For cross-species predictions, docking scores are correlated with experimentally measured binding affinities (Kd values) from surface plasmon resonance or biolayer interferometry assays [<a href="#ref-20">20</a>].
Molecular Dynamics Simulations of RBD-ACE2 Complexes
Molecular dynamics simulations provide a time-resolved view of the RBD-ACE2 interaction by solving Newton's equations of motion for all atoms in the system [<a href="#ref-21">21</a>]. MD simulations capture conformational changes, hydrogen bond dynamics, and solvent effects that are not accessible through static docking [<a href="#ref-22">22</a>]. Typical simulation timescales range from 100 nanoseconds to several microseconds for RBD-ACE2 complexes [<a href="#ref-23">23</a>].
The simulation setup includes solvation of the complex in a water box with explicit water models (e.g., TIP3P) and addition of counterions to neutralize the system [<a href="#ref-24">24</a>]. Force fields such as CHARMM36, AMBER ff14SB, or OPLS-AA are used to parameterize the protein atoms [<a href="#ref-25">25</a>]. The system is energy minimized, equilibrated in the NVT and NPT ensembles, and then subjected to production runs [<a href="#ref-26">26</a>].
Analysis of MD trajectories focuses on several metrics. Root-mean-square fluctuation (RMSF) identifies flexible regions in the RBM and ACE2 interface [<a href="#ref-27">27</a>]. Hydrogen bond occupancy quantifies the stability of specific inter-residue contacts over the simulation [<a href="#ref-28">28</a>]. Principal component analysis (PCA) reveals collective motions of the RBD that may be important for receptor recognition [<a href="#ref-29">29</a>]. Markov state models (MSMs) can be constructed from long MD trajectories to identify metastable conformational states and transition pathways [<a href="#ref-30">30</a>].
For cross-species studies, MD simulations are performed for RBD-ACE2 complexes from different host species. Comparative analysis of interaction energies and conformational dynamics reveals species-specific adaptations [<a href="#ref-31">31</a>]. For example, simulations of bat coronavirus RBDs with human ACE2 have identified mutations that stabilize the interface and increase binding affinity [<a href="#ref-32">32</a>].
Binding Free Energy Calculations
Binding free energy calculations provide a quantitative estimate of the strength of the RBD-ACE2 interaction [<a href="#ref-33">33</a>]. The most widely used method for this purpose is the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) approach [<a href="#ref-34">34</a>]. MM-PBSA calculates the free energy of binding as the difference between the free energies of the complex, receptor, and ligand in solution [<a href="#ref-35">35</a>].
The MM-PBSA workflow involves extracting snapshots from an MD trajectory, calculating the molecular mechanics energy (EMM) in the gas phase, adding the polar solvation free energy from the Poisson-Boltzmann equation, and adding the nonpolar solvation free energy from a solvent-accessible surface area (SASA) term [<a href="#ref-36">36</a>]. The entropic contribution is often estimated using normal mode analysis or omitted in relative binding free energy comparisons [<a href="#ref-37">37</a>].
An alternative method is the Linear Interaction Energy (LIE) approach, which uses empirical scaling factors for electrostatic and van der Waals interaction energies [<a href="#ref-38">38</a>]. Thermodynamic integration (TI) and free energy perturbation (FEP) methods offer higher accuracy but are computationally more expensive [<a href="#ref-39">39</a>].
MM-PBSA calculations have been applied to predict the effect of RBD mutations on ACE2 binding across species [<a href="#ref-40">40</a>]. Studies have shown that mutations such as N501Y and K417N in the SARS-CoV-2 RBD increase binding affinity to murine and canine ACE2 orthologs [<a href="#ref-41">41</a>]. These predictions have been validated by experimental binding assays, demonstrating the utility of MM-PBSA for cross-species risk assessment [<a href="#ref-42">42</a>].
Machine Learning Classifiers for Spillover Risk Prediction
Machine learning (ML) methods have been increasingly applied to predict host tropism and spillover risk from sequence and structural features [<a href="#ref-43">43</a>]. ML classifiers are trained on datasets of known RBD-ACE2 interactions, with features derived from sequence alignments, structural descriptors, and physicochemical properties [<a href="#ref-44">44</a>].
Feature engineering for ML models includes one-hot encoding of amino acid sequences, evolutionary conservation scores from position-specific scoring matrices (PSSMs), and structural features such as solvent accessibility, secondary structure propensity, and residue depth [<a href="#ref-45">45</a>]. Graph neural networks (GNNs) have been used to represent protein-protein interfaces as graphs, where nodes represent residues and edges represent spatial contacts [<a href="#ref-46">46</a>].
Common ML algorithms for this task include random forests, support vector machines (SVMs), gradient boosting machines (e.g., XGBoost), and deep neural networks [<a href="#ref-47">47</a>]. The output is typically a binary classification (binding or non-binding) or a continuous score representing binding affinity [<a href="#ref-48">48</a>]. Model performance is evaluated using metrics such as area under the receiver operating characteristic curve (AUC-ROC), precision-recall curves, and cross-validation [<a href="#ref-49">49</a>].
ML models have been used to screen large sequence databases of bat coronaviruses for RBD variants with high predicted affinity to livestock ACE2 orthologs [<a href="#ref-50">50</a>]. These predictions guide targeted surveillance and experimental validation efforts [<a href="#ref-51">51</a>]. Deep learning models, including convolutional neural networks (CNNs) applied to contact maps, have further improved prediction accuracy [<a href="#ref-52">52</a>].
Key Structural Motifs and Mutational Landscapes
Several structural motifs in the coronavirus RBD are critical for cross-species receptor binding. The RBM loop region, which contains the majority of contact residues, exhibits high sequence variability across coronavirus lineages [<a href="#ref-53">53</a>]. In SARS-CoV-2, the RBM contains a beta-hairpin motif that inserts into a groove on the ACE2 surface [<a href="#ref-54">54</a>]. The presence of a furin cleavage site at the S1/S2 boundary also influences host range by affecting spike protein priming [<a href="#ref-55">55</a>].
Mutational landscapes of the RBD have been systematically explored using deep mutational scanning (DMS) [<a href="#ref-56">56</a>]. DMS experiments measure the effect of every single amino acid substitution on ACE2 binding, generating comprehensive fitness maps [<a href="#ref-57">57</a>]. Computational models trained on DMS data can predict the impact of novel mutations on cross-species binding [<a href="#ref-58">58</a>].
For bat coronaviruses, key mutations that enable binding to non-bat ACE2 orthologs include changes at residues 493, 498, and 501 [<a href="#ref-59">59</a>]. The substitution Q493H has been shown to enhance binding to human ACE2 by introducing a favorable electrostatic interaction [<a href="#ref-60">60</a>]. Similarly, the N501Y mutation increases hydrophobic contacts with ACE2 residues Y41 and K353 [<a href="#ref-61">61</a>].
Integration of Computational and Experimental Data
A robust computational pipeline for cross-species receptor binding prediction integrates multiple methods in a hierarchical workflow [<a href="#ref-62">62</a>]. The pipeline begins with sequence-based screening using ML classifiers to identify high-risk RBD variants [<a href="#ref-63">63</a>]. Candidate variants are then subjected to molecular docking to generate initial binding poses [<a href="#ref-64">64</a>]. The top-ranked complexes are refined using MD simulations, and binding free energies are calculated using MM-PBSA [<a href="#ref-65">65</a>]. Finally, predictions are validated through experimental assays such as surface plasmon resonance or pseudovirus entry assays [<a href="#ref-66">66</a>].
The following Mermaid diagram illustrates this integrated workflow:
flowchart TD
A["Viral Sequence Database"] --> B["Sequence Feature Extraction"]
B --> C["Machine Learning Classifier"]
C --> D["High-Risk RBD Variants"]
D --> E["Molecular Docking with Host ACE2 Orthologs"]
E --> F["Scoring and Pose Selection"]
F --> G["Molecular Dynamics Simulations"]
G --> H["Trajectory Analysis RMSF, H-Bonds, PCA"]
H --> I["MM-PBSA Binding Free Energy Calculation"]
I --> J["Predicted Binding Affinity"]
J --> K["Experimental Validation SPR, Pseudovirus Assays"]
K --> L["Spillover Risk Assessment"]
Applications to Bat, Avian, and Mammalian Hosts
Coronaviruses circulating in bat populations represent a major reservoir for emerging zoonotic viruses [<a href="#ref-67">67</a>]. Computational studies have focused on bat SARS-like coronaviruses (SL-CoVs) such as RaTG13, WIV1, and SHC014, which show variable binding to human and livestock ACE2 [<a href="#ref-68">68</a>]. Docking and MD simulations have identified key residues in the bat RBD that must mutate to enable efficient binding to swine or bovine ACE2 [<a href="#ref-69">69</a>].
Avian coronaviruses, including infectious bronchitis virus (IBV) in poultry, use different receptors such as APN or sialic acids [<a href="#ref-70">70</a>]. Computational modeling of IBV spike-receptor interactions has been used to predict host range shifts between galliform and anseriform birds [<a href="#ref-71">71</a>]. The structural comparison of avian versus mammalian receptor binding is covered in detail in the article Structural Comparison of Avian Versus Mammalian Influenza Receptor Binding.
For mammalian livestock species, computational predictions have been made for RBD binding to ACE2 orthologs from swine, cattle, horses, and companion animals such as cats and dogs [<a href="#ref-72">72</a>]. These predictions inform risk assessments for reverse zoonosis (spillback) events, where human-adapted coronaviruses transmit to animals [<a href="#ref-73">73</a>]. The article Computational Prediction of Host Tropism and Receptor Binding Dynamics in Emerging Zoonotic Coronaviruses provides further context on host tropism prediction.
Limitations and Challenges
Despite significant advances, computational prediction of cross-species receptor binding faces several limitations. Force field inaccuracies can lead to errors in binding free energy estimates [<a href="#ref-74">74</a>]. Solvent effects, particularly the role of water molecules at the protein-protein interface, are often inadequately modeled [<a href="#ref-75">75</a>]. The conformational flexibility of glycans on the spike protein, which can modulate receptor accessibility, is challenging to incorporate in simulations [<a href="#ref-76">76</a>].
Experimental validation remains essential, as computational predictions can produce false positives or false negatives [<a href="#ref-77">77</a>]. The availability of high-resolution structures for diverse ACE2 orthologs is limited, necessitating homology modeling which introduces additional uncertainty [<a href="#ref-78">78</a>]. Receptor binding is only one factor in host tropism; post-entry factors such as viral replication, immune evasion, and host proteases also determine spillover success [<a href="#ref-79">79</a>].
Future Directions
Emerging computational methods promise to improve the accuracy and throughput of cross-species binding predictions. AlphaFold2 and related deep learning models enable accurate prediction of RBD and ACE2 structures without experimental templates [<a href="#ref-80">80</a>]. Protein language models, such as ESM-1b and ProtBERT, can learn evolutionary constraints directly from sequence data and predict mutation effects on binding [<a href="#ref-81">81</a>].
Enhanced sampling techniques, including replica exchange MD and metadynamics, allow exploration of rare conformational events relevant to receptor binding [<a href="#ref-82">82</a>]. Coarse-grained MD simulations enable simulation of larger systems over longer timescales, facilitating the study of spike trimer-receptor interactions [<a href="#ref-83">83</a>]. The integration of cryo-electron microscopy data with MD simulations, as discussed in Integrating Cryo-EM and Molecular Dynamics Simulations to Elucidate Glycan Shield Dynamics in Emerging Zoonotic Coronaviruses, provides a more complete picture of spike dynamics.
Conclusion
Computational prediction of cross-species receptor binding dynamics is a critical component of zoonotic coronavirus risk assessment. Molecular docking, MD simulations, MM-PBSA calculations, and machine learning classifiers each contribute unique insights into the molecular determinants of host tropism. The integration of these methods into a unified pipeline, combined with experimental validation, enables the identification of high-risk RBD variants circulating in animal reservoirs. Continued development of computational tools and expansion of structural databases will further enhance our ability to predict and prevent future zoonotic spillover events.