Spike Protein Dynamics and Host Range Prediction Using Molecular Dynamics and Machine Learning
By Dr. Zubair Khalid, DVM, MS, PhD ·

Key Takeaways
- Molecular Dynamics (MD) simulations provide atomistic resolution of spike (S) protein conformational changes, revealing how mutations like D614G and N501Y alter flexibility and enhance receptor binding affinity by influencing thermodynamic stability and allosteric communication.
- Machine Learning (ML) classifiers, trained on structural and sequence features derived from MD simulations and genomic data, enable quantitative prediction of host range and cross-species transmission potential by identifying key determinants of species tropism.
- Glycan shielding and specific N-glycan modifications, such as those observed in PEDV, significantly influence RBD opening and immune evasion, impacting host recognition and viral infectivity, which can be modeled using integrated computational approaches.
- The integration of MD-derived conformational ensembles with ML classifiers facilitates a systematic evaluation of zoonotic spillover risk by analyzing critical structural features like receptor-binding interface composition, loop dynamics, allosteric networks, and sialic acid binding capacity.
- This computational framework is directly applicable to veterinary surveillance, allowing for the prediction of host tropism in animal coronaviruses and other viral pathogens, thereby guiding proactive risk assessment and pandemic preparedness efforts.
Introduction
The spike (S) glycoprotein of enveloped viruses mediates host cell attachment and membrane fusion, constituting the primary determinant of species tropism and zoonotic potential [<a href="#ref-1">1</a>, <a href="#ref-2">2</a>]. In coronaviruses, the S protein exists as a trimeric class I fusion protein whose large ectodomain comprises a receptor-binding subunit (S1) and a fusion subunit (S2) [<a href="#ref-3">3</a>, <a href="#ref-4">4</a>]. Conformational transitions between closed (prefusion) and open (receptor-accessible) states govern infectivity and are influenced by mutations that alter thermodynamic stability, glycan shielding, and allosteric communication [<a href="#ref-5">5</a>, <a href="#ref-6">6</a>, <a href="#ref-7">7</a>]. Molecular dynamics (MD) simulations provide atomistic resolution of these processes, while machine learning (ML) classifiers trained on structural and sequence features enable quantitative host range prediction [<a href="#ref-8">8</a>, <a href="#ref-9">9</a>]. This article reviews the integrated computational framework for analyzing spike protein dynamics and forecasting cross-species transmission in veterinary contexts, with emphasis on coronaviruses affecting swine, poultry, and companion animals.
Molecular Dynamics Simulations of Spike Protein Conformational Dynamics
MD simulations have been instrumental in characterizing the conformational landscape of spike proteins. Full-length, fully glycosylated models embedded in lipid bilayers reveal that the stalk region contains two independent joints that allow the head domain to sample orientations compatible with receptor engagement [<a href="#ref-10">10</a>]. The receptor-binding domain (RBD) within the S1 subunit undergoes hinge-like movements that expose or conceal the receptor-binding motif [<a href="#ref-4">4</a>, <a href="#ref-11">11</a>]. Coarse-grained normal-mode analyses of 17,081 SARS-CoV-2 spike variants demonstrate that mutations such as D614G increase open-state occupancy by enhancing closed-state flexibility and reducing open-state flexibility [<a href="#ref-7">7</a>]. All-atom simulations further show that the N501Y mutation in the RBD introduces hydrophobic contacts with Y41 and K353 of human ACE2, raising binding affinity by approximately 0.81 kcal/mol [<a href="#ref-12">12</a>].
The interplay between furin cleavage site mutations and the D614G substitution modulates S protein dynamics by altering the conformational equilibrium of the S1/S2 interface [<a href="#ref-5">5</a>]. Candidate adaptive polymorphisms (CAPs) identified through evolutionary models exhibit allosteric coupling with ACE2-binding residues, with Delta strain variants showing epistatic interactions that rigidify the RBD, whereas Omicron variants accumulate additive mutations that preserve flexibility [<a href="#ref-6">6</a>]. MD simulations of variant spikes from Pakistani strains (I210T, V213G, S371F, S373P, T478K, F486V, Y505H, D796Y) reveal that Y505H significantly enhances binding to both ACE2 and neuropilin-1 (NRP1), with mutant complexes stabilizing earlier than wild-type in 40 ns trajectories [<a href="#ref-13">13</a>].
Glycan moieties modulate RBD opening and immune evasion. Simulations of the PEDV spike (alphacoronavirus) integrated with cryo-electron tomography and mass spectrometry demonstrate that a key N-glycan at the D0 domain controls conformational switching necessary for receptor binding [<a href="#ref-14">14</a>]. In SARS-CoV-2, glycans stabilize specific open and closed states; antibody epitope exposure is more accurately predicted by glycan-antibody clash analysis than by solvent accessibility alone [<a href="#ref-10">10</a>]. The N-terminal domain (NTD) contains a sialic acid-binding site that exhibits short-lived interactions, contributing to avidity-driven attachment [<a href="#ref-15">15</a>].
External perturbations such as electric fields and temperature also influence spike dynamics. Prefusion RBDs are highly susceptible to structural damage by moderate electric fields (10^4 V/m), whereas the postfusion S2 conformation remains intact at four orders of magnitude higher field strength [<a href="#ref-16">16</a>]. Febrile temperatures (311-312 K) enhance the binding affinity of monoclonal antibodies to the RBD by increasing intraprotein dynamics and conformational sampling of complementary epitopes [<a href="#ref-17">17</a>].
Machine Learning for Binding Affinity and Host Range Prediction
ML models exploit sequence and structural features to predict spike, receptor interactions and host tropism. Deep learning-based QSAR models trained on approved drug libraries (4,388 compounds) identified protein, protein interaction modulators that were subsequently docked against the Delta spike RBD; fexofenadine, an antihistamine, showed stable binding in MD simulations and MM/GBSA calculations [<a href="#ref-8">8</a>]. Virtual screening of 2,456 approved drugs against the S1, ACE2 interface using hierarchical docking (HTVS, SP, XP) and 100 ns MD simulations revealed four stable binders (riboflavin, fenoterol, cangrelor, vidarabine) that interact with physiologically important interfacial residues.
Graph-theoretic frameworks provide a multiscale representation of spike structure. A three-level model incorporating 19 top-level vertices, intermediate proximity graphs (6 Å cutoff), and molecular descriptors identified N501Y and L452R as mutations with the most pronounced effects on RBD conformation and stability, corroborated by I-TASSER modeling and MD [<a href="#ref-9">9</a>]. Immunoinformatics pipelines combining B-cell, CD4+, and CD8+ T-cell epitope prediction with MD-validation have been used to design multi-epitope vaccine constructs; one such construct containing 4 CD4+ and 4 CD8+ epitopes linked to adjuvant proteins showed strong binding to TLR3 and TLR4 in simulations [<a href="#ref-18">18</a>].
Predictive models of antibody dynamics following infection use hierarchical gamma models to distinguish between continuous decay and decay-to-plateau kinetics. In healthcare workers, spike (S) antibodies persisted with a predicted half-life exceeding 465 days, correlating tightly with surrogate neutralization (R^2 = 0.72), whereas nucleoprotein antibodies decayed with a half-life of 60 days [<a href="#ref-19">19</a>]. Such frameworks can be adapted to veterinary species to estimate duration of protective immunity after natural infection or vaccination.
Integration of MD and ML for Zoonotic Spillover Risk Assessment
Combining MD-derived conformational ensembles with ML classifiers enables systematic evaluation of host range determinants. Key structural features that govern species tropism include:
- Receptor-binding interface composition: stable hydrogen bonds and salt bridges that differ between SARS-CoV-2 and SARS-CoV [<a href="#ref-20">20</a>].
- Loop dynamics: a flexible loop (Loop 3) in the unbound RBD can adopt substates that occlude the ACE2-binding interface, representing cryptic targets for therapeutic design [<a href="#ref-11">11</a>].
- Allosteric networks: CAP sites in the spike control the dynamics of binding residues in the open state, suggesting mechanisms for adaptation to new hosts [<a href="#ref-6">6</a>].
- Sialic acid binding capacity: NTD insertions in SARS-CoV-2, absent in SARS-CoV, confer attachment to sialylated glycans on host cells [<a href="#ref-15">15</a>].
The following workflow (Figure 1) integrates these components into a prediction pipeline:
flowchart TD
A["Viral genome sequence databases: GISAID, NCBI"] --> B["Homology modeling / AlphaFold2 structure prediction"]
B --> C["All-atom or coarse-grained MD simulations"]
C --> D["Conformational ensemble analysis: RMSD, RMSF, PCA, FEL"]
D --> E["Feature extraction: binding free energy, interface contacts, allosteric coupling"]
E --> F["ML classifier training: random forest, deep neural network, graph network"]
F --> G["Host range prediction: probability of ACE2 binding, receptor ortholog compatibility"]
G --> H["Risk assessment: spillover potential, surveillance prioritization"]
Figure 1. Integrated computational workflow for host range prediction using MD simulations and machine learning.
Table 1 summarizes selected mutations and their effects as resolved by MD studies.
| Mutation | Variant Lineage | Effect on Dynamics | Binding Affinity Change | Reference |
|---|---|---|---|---|
| D614G | Early pandemic | Increases open-state occupancy; allosterically modulates furin site | Neutral or slight increase | [<a href="#ref-5">5</a>, <a href="#ref-7">7</a>] |
| N501Y | Alpha, Beta, Gamma | Introduces hydrophobic contacts with ACE2 Y41, K353 | +0.81 kcal/mol | [<a href="#ref-12">12</a>] |
| T478K | Delta | Enhances electrostatic interactions with ACE2 | Moderate increase | [<a href="#ref-13">13</a>] |
| Y505H | Pakistani isolates | Disrupts protein function; stabilizes RBD, ACE2 complex | Large increase | [<a href="#ref-13">13</a>] |
| E484Q+L452R | Kappa, Delta | Reduces loop flexibility; increases hydrophobic packing | Moderate increase | [<a href="#ref-21">21</a>] |
Table 1. Impact of spike protein mutations on conformational dynamics and receptor binding.
Implications for Veterinary Surveillance and Pandemic Preparedness
The computational tools described here are directly applicable to veterinary virology. For example, the PEDV spike structure solved by cryo-ET revealed domain motions and glycosylation patterns that affect host recognition in swine [<a href="#ref-14">14</a>]. Similarly, MD simulations of avian influenza hemagglutinin (HA) have been used to predict receptor-binding specificity for avian versus mammalian sialic acid linkages see related article: [Structural Dynamics of Avian Influenza Hemagglutinin: Molecular Modeling and Receptor Binding Predictions for Pandemic Risk Assessment]. For bat coronaviruses, computational predictions of ACE2 binding affinity using MD and ML can flag strains with high zoonotic potential see related article: [Predicting Zoonotic Spillover: Computational Modeling of Receptor-Binding Dynamics in Emerging Bat Coronaviruses].
Sequence databases such as GISAID and NCBI GenBank provide the raw material for variant surveillance. Routine reanalysis of spike protein dynamics with newly emerging mutations can be automated using Markov state models or deep learning surrogate models to avoid costly all-atom simulations see related article: [Markov State Models in Molecular Dynamics Simulations]. Predictive models of antibody escape, such as those based on deep mutational scanning, can be integrated with MD to forecast vaccine failure in livestock and companion animals see related article: [Deep Mutational Scanning and Machine Learning for Predicting SARS-CoV-2 Spike Protein Evolution and Antibody Escape].
Conclusion
Molecular dynamics simulations and machine learning classifiers constitute a powerful synergistic framework for predicting host range from spike protein structure and dynamics. By resolving the atomic-level mechanisms of receptor binding, conformational transitions, and mutational effects, these methods enable proactive risk assessment of emerging viral pathogens in veterinary populations. Continued integration with cryo-EM, glycomics, and large-scale sequence surveillance will further refine spillover predictions and guide countermeasure development.
Related Clinical & Scientific Guides
- A Practical Guide to Detecting Antimicrobial Resistance Genes in Shotgun Metagenomic Data
- Computational Immunology: Modeling the Immune System
- How to Set Hard Filters for Germline Variant Calling: A Practical Guide to GATK Best Practices