Predicting Receptor Binding Specificity of Zoonotic Influenza A Viruses Using Molecular Dynamics Simulations

1. Introduction

Influenza A viruses (IAVs) circulate in a wide range of avian and mammalian hosts, with waterfowl and shorebirds serving as the primary natural reservoir [1]. Zoonotic transmission of IAVs from avian species to mammals, including swine, mustelids, and occasionally humans, poses a continuous pandemic threat [1, 2]. The host range of an IAV is largely determined by the binding specificity of its major surface glycoprotein, hemagglutinin (HA), to sialic acid (Sia) receptors on host epithelial cells [2]. Avian-adapted HAs preferentially bind to Sia linked to galactose via an α2-3 linkage (Siaα2-3Gal), whereas human-adapted HAs bind to Siaα2-6Gal [3]. This tropism barrier is a critical bottleneck for cross-species transmission [1, 3].

Molecular dynamics (MD) simulations provide a powerful computational framework for investigating the atomic-level interactions between HA and host glycan receptors [2]. By simulating the conformational dynamics of HA–receptor complexes in explicit solvent, MD can predict binding affinities, identify key residues that determine specificity, and assess the impact of mutations that may shift receptor preference [3]. This article discusses the theoretical basis of HA–receptor interactions, the workflow of MD simulations applied to IAV receptor binding, and the implications of these computational predictions for veterinary surveillance and pandemic risk assessment.

2. Biological and Structural Basis of HA–Receptor Interactions

2.1 Hemagglutinin Structure

Hemagglutinin is a homotrimeric glycoprotein embedded in the viral envelope [1]. Each monomer consists of two domains: a globular head domain carrying the receptor binding site (RBS) and a stem domain responsible for membrane fusion [1, 2]. The RBS is a shallow pocket formed by three structural elements: the 190-helix (residues 190–198), the 130-loop (residues 134–138), and the 220-loop (residues 221–228) [3]. Within the RBS, a set of conserved residues directly contacts the sialic acid moiety, including Tyr98, Trp153, His183, and Tyr195 (H3 numbering) [2]. A secondary binding site, termed the vestigial esterase domain, may also interact with extended glycan chains [3].

2.2 Sialic Acid Receptors

Sialic acids are nine-carbon monosaccharides that cap glycans on host cell surfaces [1]. The linkage between Sia and the penultimate galactose is a key determinant of host tropism: avian enteric epithelia predominantly express Siaα2-3Gal, while human respiratory epithelia express both Siaα2-3Gal and Siaα2-6Gal, with the latter predominant in the upper airways [2]. Swine tracheal epithelium contains both linkages, making swine potential mixing vessels for reassortment [3]. The fine structure of the glycan (e.g., fucosylation, sulfation, and the underlying sugar chain) further modulates binding [1].

2.3 Key Residues Determining Specificity

Mutational and structural studies have identified several amino acid positions in the RBS that govern receptor preference [2, 3]. For example, substitutions at residues 226 and 228 (H3 numbering) are strongly associated with α2-6 specificity in human seasonal H3N2 and H2N2 viruses [2]. Avian H5N1 and H7N9 HAs typically possess Gln226 and Gly228, which favor α2-3 binding; changes to Leu226 and Ser228 are linked to increased α2-6 binding [3]. Other residues, such as Glu190, Asp225, and Ser193, also contribute to linkage preference by forming additional hydrogen bonds with the glycan [1, 2].

Table 1 summarizes the major RBS residues and their linkage preference in selected IAV subtypes.

Table 1: Selected HA RBS Residues and Receptor Preference

3. Molecular Dynamics Simulations of HA–Receptor Complexes

3.1 Overview of the MD Workflow

Molecular dynamics simulations numerically integrate Newton's equations of motion for a system of atoms over time, typically on the nanosecond to microsecond scale [2]. For HA–receptor binding prediction, the workflow includes the following steps:

1. Selection of high-resolution crystal or cryo-EM structures of HA from relevant subtypes (e.g., H5, H7, H9) [3].
2. Preparation of glycan receptor models representing α2-3-linked sialyllactose (3-SL) and α2-6-linked sialyllactose (6-SL) [1].
3. Docking or manual placement of the glycan into the RBS, followed by solvation in a water box with neutralizing ions [2].
4. Energy minimization, equilibration, and production runs using a classical force field (e.g., CHARMM, AMBER, GROMACS) [3].
5. Analysis of trajectory data including root-mean-square deviation (RMSD), hydrogen bond occupancy, binding free energy via MM/GBSA or MM/PBSA, and per-residue energy decomposition [2].

The following Mermaid diagram illustrates this pipeline.

flowchart TD
    A[Select HA crystal structure], > B[Prepare glycan receptor models (3-SL, 6-SL)]
    B, > C[Dock glycan into RBS or align from known complex]
    C, > D[Solvate, add ions, parameterize with force field]
    D, > E[Energy minimization and equilibration (NVT, NPT)]
    E, > F[Production MD (ns to μs simulations)]
    F, > G[Trajectory analysis: RMSD, H-bonds, MM/GBSA, per-residue energy]
    G, > H{Predict dominant receptor preference}
    H, > I[α2-3 favored], avian-like
    H, > J[α2-6 favored], mammalian-like
    H, > K[Dual or intermediate], flexibility

3.2 Force Fields and Glycan Parameterization

The accuracy of MD simulations depends on the quality of the force field parameters used for proteins and glycans. Common force fields for protein simulations include CHARMM36 and AMBER ff14SB [3]. For carbohydrates, force fields such as GLYCAM06 or CHARMM-GUI glycan builder provide comprehensive parameter sets [2]. Simulation temperature should be maintained at 310 K to approximate mammalian body temperature, while avian body temperature (around 313–315 K) can be simulated for comparative studies [1].

3.3 Binding Free Energy Calculations

End-point free energy methods such as Molecular Mechanics Generalized Born Surface Area (MM/GBSA) and Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) are widely used to estimate relative binding affinities of HA for different glycans [2]. These methods compute the difference between the free energy of the complex and the sum of free energies of the free protein and free ligand in solution [3]. MM/GBSA can rank the preference of an HA for α2-3 versus α2-6 glycans with moderate computational cost [2]. More rigorous approaches, such as free energy perturbation (FEP) or thermodynamic integration, provide higher accuracy but are computationally intensive [3].

3.4 Analysis of Key Interactions

Hydrogen bond occupancy analysis reveals which residues contribute most to stabilizing the bound glycan. For avian-adapted HAs, a network of hydrogen bonds typically involves Gln226, and the α2-3 linkage forms a bent conformation that fits snugly into the RBS [2]. For human-adapted HAs, Leu226 and Ser228 promote a more extended conformation of α2-6 glycans, with additional interactions from residues Gly225 and Glu190 [3]. MD simulations can capture the dynamic formation and breakage of these contacts, offering insights into the role of water-mediated bridges [1].

4. Predicting Host Tropism and Pandemic Risk

4.1 Computational Signatures of Zoonotic Potential

A virus that acquires mutations enabling α2-6 binding while retaining replicative fitness in mammalian respiratory tract is considered to have heightened pandemic potential [1, 2]. MD simulations can identify such mutations before they appear in surveillance data. For example, the Q226L and G228S substitutions in H5 and H7 HAs have been simulated to show enhanced α2-6 binding [3]. Simulations of HA from avian H9N2 viruses have indicated that the T189A and Q227M mutations (H3 numbering) shift preference toward human-type receptors [2].

4.2 Integration with Sequence Surveillance

MD predictions complement phylogenetic and deep mutational scanning data by providing structural context [2, 3]. When novel HA sequences from field samples are obtained, homology modeling can generate 3D structures for MD assessment of receptor binding [3]. This approach is particularly valuable for emerging subtypes such as H10N8, H5N6, and H7N4, where experimental glycan array data may be limited [1].

4.3 Limitations and Considerations

While MD simulations offer atomic-level detail, they are subject to limitations [2]. Force field inaccuracies for glycans, insufficient sampling of conformational states, and the absence of the full viral membrane and host cell context can affect predictive accuracy [3]. Simulations of large systems (e.g., full-length HA in a membrane) are computationally demanding and often require coarse-grained methods [1]. Therefore, MD predictions should be validated with experimental binding assays such as glycan arrays and surface plasmon resonance [2].

5. Implications for Veterinary Medicine and One Health

Understanding HA receptor specificity is crucial for veterinary surveillance of IAVs in poultry, swine, and other livestock [1]. Avian influenza strains that acquire even partial mammalian receptor binding ability require enhanced monitoring and biosecurity measures [2]. MD simulations can be used to pre-screen viral isolates from outbreaks to assess their zoonotic risk, guiding decisions on culling, vaccination, and trade restrictions [3].

The computational workflow described here is conceptually similar to methods applied to other zoonotic viruses, such as the modeling of coronavirus spike protein interactions with host receptors, as detailed in related articles on this portal (e.g., Computational Insights into Host Receptor Binding Dynamics of Emerging Zoonotic Coronaviruses Using Molecular Dynamics Simulations and Molecular Modeling of Zoonotic Spillover: Predicting Receptor-Binding Dynamics in Bat-Borne Coronaviruses). The integration of MD with deep learning and large-scale mutational scanning, as discussed in Deep Learning for Predicting Receptor-Binding Domain Dynamics in Emerging Zoonotic Coronaviruses, further enhances predictive power.

6. Conclusion

Molecular dynamics simulations provide a powerful, hypothesis-driven approach to predict receptor binding specificity of zoonotic influenza A viruses. By modeling the atomic interactions between HA and host glycans, these simulations can identify key residues that govern host tropism and assess the impact of mutations that may facilitate cross-species transmission. Although computational predictions require experimental validation, they serve as a valuable tool for veterinary surveillance, outbreak risk assessment, and pandemic preparedness. Continued advances in force field accuracy, simulation length, and integration with machine learning will strengthen the role of MD simulations in computational virology.

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.

References

[1] Swayne DE, Suarez DL, Sims LD. Influenza A viruses. In: Swayne DE, editor. Diseases of Poultry. 14th ed. Ames, IA: Wiley-Blackwell; 2020. p. 210–256.

[2] Zambon MC, Ellis JS. Influenza virology and the role of hemagglutinin receptor binding. In: Webster RG, Monto AS, Braciale TJ, Lamb RA, editors. Textbook of Influenza. 2nd ed. Oxford: Wiley-Blackwell; 2013. p. 35–51.

[3] Taubenberger JK, Kash JC. Influenza hemagglutinin structure and receptor binding specificity. In: Compans RW, Oldstone MBA, editors. Influenza: Pathogenesis and Control – Volume I. Current Topics in Microbiology and Immunology. Vol. 385. Cham: Springer; 2014. p. 1–30.