Quantifying Binding Affinity Changes in Influenza Neuraminidase Mutants Using Free Energy Perturbation

Introduction

Influenza A viruses circulate widely in avian and swine populations, causing significant economic losses in poultry and livestock operations and posing zoonotic spillover risks [1]. The viral neuraminidase (NA) glycoprotein is a primary target for antiviral intervention in veterinary medicine. Neuraminidase inhibitors (NAIs) such as oseltamivir and zanamivir are used to control influenza outbreaks in commercial poultry and swine herds [2]. However, the emergence of NA mutations that confer reduced susceptibility to these inhibitors undermines treatment efficacy and complicates outbreak management [3]. Accurate prediction of how specific amino acid substitutions alter NAI binding affinity is essential for proactive resistance surveillance and informed antiviral stewardship [4].

Free energy perturbation (FEP) is a rigorous computational method rooted in statistical mechanics that calculates the relative binding free energy difference between a wild-type and a mutant protein-ligand complex [5]. When applied to influenza NA, FEP can quantify the impact of resistance-associated mutations on inhibitor binding. This article provides a detailed technical review of FEP methodology, its application to NA mutants, and its integration into veterinary computational virology workflows.

Biological Context: Influenza Neuraminidase and Inhibitor Binding

Influenza NA is a homotetrameric sialidase that cleaves terminal sialic acid residues from host cell glycoproteins, facilitating viral release from infected cells and preventing viral aggregation [6]. The NA active site is a highly conserved pocket composed of catalytic residues (Arg118, Asp151, Arg152, Arg224, Glu276, Arg292, Arg371, and Tyr406) and framework residues that stabilize the catalytic architecture [7]. NAIs are transition-state analog inhibitors that competitively occupy this active site, blocking enzymatic activity and halting viral spread [8].

Mutations in the NA active site can reduce inhibitor binding through steric hindrance, electrostatic repulsion, or altered hydrogen bonding networks [9]. Common resistance mutations in avian and swine influenza NA include H274Y (N2 numbering), N294S, and E119V, among others [10]. These substitutions often reduce the binding affinity of oseltamivir carboxylate by several kilocalories per mole, corresponding to a 10- to 100-fold increase in the half-maximal inhibitory concentration (IC50) [11]. Quantifying these affinity changes with high precision is critical for distinguishing clinically relevant resistance from benign polymorphisms [12].

Principles of Free Energy Perturbation

FEP is grounded in the Zwanzig equation, which relates the free energy difference between two thermodynamic states (A and B) to the ensemble average of the Boltzmann factor of their energy difference [13]:

ΔG(A→B) = -kBT ln ⟨exp[-(VB - VA)/kBT]⟩A

In the context of protein-ligand binding, the two states represent the wild-type and mutant complexes. The calculation is typically performed in a multistage approach, dividing the transformation into a series of intermediate (lambda) windows to ensure adequate phase space overlap [14]. Each window corresponds to a mixed Hamiltonian where the potential energy is a linear combination of the end-state potentials:

V(λ) = (1-λ)VA + λVB

Molecular dynamics (MD) simulations are conducted at each lambda value, and the free energy change is accumulated using thermodynamic integration (TI) or the Bennett acceptance ratio (BAR) method [15]. For binding affinity calculations, a thermodynamic cycle is constructed that includes the free protein, the free ligand, the protein-ligand complex, and the solvated environment [16]. The relative binding free energy (ΔΔGbind) is obtained as:

ΔΔGbind = ΔGbind(mutant) - ΔGbind(wild-type) = ΔGcomplex - ΔGfree

where ΔGcomplex is the free energy of mutating the ligand or protein in the bound state, and ΔGfree is the free energy of the same mutation in the unbound state [17].

Computational Setup for Neuraminidase FEP

System Preparation

The starting structure for NA FEP calculations is typically obtained from X-ray crystallography or cryo-electron microscopy of the NA-inhibitor complex [18]. For veterinary influenza strains, structures of avian N1, N2, N5, N7, and N9 subtypes are available in the Protein Data Bank [19]. The protein is prepared by adding hydrogen atoms, assigning protonation states at physiological pH using tools such as PROPKA, and resolving steric clashes [20]. The inhibitor (e.g., oseltamivir carboxylate) is parameterized using a force field compatible with the protein model, such as the General Amber Force Field (GAFF) or the CHARMM General Force Field (CGenFF) [21].

Force Field Selection

The accuracy of FEP calculations depends critically on the force field. For NA, the Amber ff14SB or CHARMM36m force fields are commonly used for the protein, while the ligand parameters are derived from quantum mechanical calculations at the HF/6-31G* or B3LYP/6-31G* level of theory [22]. Water molecules are modeled using the TIP3P or OPC explicit solvent models [23]. Counterions (Na+ or Cl-) are added to neutralize the system and achieve a physiological ionic strength of approximately 150 mM [24].

Simulation Protocol

The system is solvated in a periodic box with a minimum distance of 10 Å between the protein and the box edge [25]. Energy minimization is performed using the steepest descent algorithm followed by conjugate gradient minimization [26]. The system is then equilibrated in the NVT ensemble (constant number of particles, volume, and temperature) for 100 ps, followed by NPT equilibration (constant number of particles, pressure, and temperature) for 500 ps [27]. Production MD simulations for FEP are conducted in the NPT ensemble at 310 K using a Langevin thermostat and a Berendsen barostat [28]. A time step of 2 fs is used with the SHAKE algorithm to constrain bonds involving hydrogen atoms [29].

FEP Protocol

The mutation is introduced in silico by alchemically transforming the wild-type residue into the mutant residue. For NA resistance mutations such as H274Y, the transformation involves changing the histidine side chain to a tyrosine side chain [30]. The perturbation is performed over 11 to 21 lambda windows, with each window simulated for 1 to 5 ns [31]. The total simulation time for a single FEP calculation is typically 10 to 100 ns, depending on the system size and the number of lambda windows [32]. Free energy differences are computed using the BAR estimator, which provides lower variance than TI for systems with nonlinear energy responses [33].

Interpretation of Results

The output of an FEP calculation is the relative binding free energy (ΔΔGbind) expressed in kcal/mol. A positive ΔΔGbind indicates that the mutation reduces binding affinity (i.e., the inhibitor binds less tightly to the mutant than to the wild type), while a negative value indicates increased affinity [34]. For NAI resistance, a ΔΔGbind greater than +1.0 kcal/mol is generally considered clinically relevant, as it corresponds to an approximately 5-fold increase in IC50 [35]. However, the threshold varies by inhibitor and viral subtype [36].

FEP results must be validated against experimental IC50 or Ki data. For veterinary NA mutants, published enzymatic inhibition assays provide reference values [37]. A mean unsigned error (MUE) of less than 1.0 kcal/mol between computed and experimental ΔΔGbind is considered acceptable for predictive applications [38]. Systematic errors can arise from inadequate sampling, force field inaccuracies, or incomplete treatment of protonation state changes [39].

Implications for Antiviral Resistance Surveillance

The integration of FEP into veterinary surveillance pipelines enables the prospective identification of resistance mutations before they become widespread in the field [40]. By screening all possible single-point mutations in the NA active site, computational virologists can generate resistance risk maps that prioritize mutations for experimental testing [41]. This approach is particularly valuable for emerging avian influenza subtypes such as H5N1, H7N9, and H9N2, where NAI susceptibility data may be limited [42].

FEP can also be used to assess the impact of compensatory mutations that restore viral fitness without directly affecting inhibitor binding [43]. For example, the NA mutation R292K confers high-level oseltamivir resistance but reduces enzymatic activity; secondary mutations such as E119G can partially restore catalytic efficiency while maintaining resistance [44]. FEP calculations on double mutants can quantify the net effect on inhibitor binding and guide the selection of alternative antivirals [45].

Workflow Overview

The following Mermaid diagram summarizes the computational workflow for quantifying binding affinity changes in NA mutants using FEP.

flowchart TD
    A[Obtain NA-inhibitor crystal structure], > B[Prepare protein and ligand parameters]
    B, > C[Solvate and neutralize system]
    C, > D[Energy minimization and equilibration]
    D, > E[Define wild-type and mutant residues]
    E, > F[Set up FEP lambda windows]
    F, > G[Run MD simulations at each lambda]
    G, > H[Compute free energy differences using BAR]
    H, > I[Calculate ΔΔGbind via thermodynamic cycle]
    I, > J{Validate against experimental data}
    J, >|MUE < 1.0 kcal/mol| K[Predict resistance risk]
    J, >|MUE > 1.0 kcal/mol| L[Refine force field or sampling]
    K, > M[Integrate into surveillance database]

Integration with Structural Visualization

The results of FEP calculations are most informative when mapped onto the three-dimensional structure of NA. Mutation sites that cause large positive ΔΔGbind values can be highlighted in the NA active site to visualize steric clashes or disrupted hydrogen bonds [46]. A 3D Protein Viewer can be linked to display the NA tetramer, the bound inhibitor, and the positions of key resistance mutations [47]. This visualization aids veterinary diagnosticians in interpreting sequencing data and assessing the likely impact of observed mutations on NAI efficacy.

Cross-Linking to Related Articles

This article is part of a broader computational virology framework. Readers are directed to the following related resources on this portal:

• Free Energy Perturbation Calculations in Drug Discovery for a general introduction to FEP methodology.
• Deep Learning for Predicting Antiviral Resistance Mutations in Influenza Neuraminidase for machine learning approaches to resistance prediction.
• Computational Prediction of Antigenic Evolution in Influenza A Virus Using Machine Learning for broader evolutionary modeling.
• Structural Comparison of Avian Versus Mammalian Influenza Receptor Binding for host range determinants.
• Avian Influenza A Virus in Poultry: Clinical Signs and Surveillance for field surveillance context.
• Swine Influenza in Pigs: Clinical Symptoms and Differential Diagnosis for swine-specific considerations.

Limitations and Future Directions

Despite its theoretical rigor, FEP has several limitations. The computational cost is high, requiring specialized hardware such as GPU clusters for routine screening of large mutation libraries [48]. Force field inaccuracies can lead to systematic biases, particularly for mutations involving charged residues or metal coordination [49]. Solvent effects and protein conformational changes on timescales beyond typical simulation lengths (microseconds to milliseconds) are not captured by standard FEP protocols [50].

Emerging methods such as alchemical absolute binding free energy calculations and enhanced sampling techniques (e.g., replica exchange MD, metadynamics) offer improved accuracy for challenging systems [51]. The combination of FEP with deep learning surrogate models may enable rapid pre-screening of mutations before expensive alchemical simulations are performed [52]. For veterinary applications, the development of force fields specifically parameterized for avian and swine NA subtypes could further improve predictive accuracy [53].

Conclusion

Free energy perturbation is a powerful computational tool for quantifying binding affinity changes in influenza neuraminidase mutants. When applied to NAI resistance, FEP provides quantitative ΔΔGbind values that can be used to prioritize mutations for experimental validation and to inform antiviral resistance surveillance in poultry and swine populations. The integration of FEP with structural visualization, molecular dynamics simulations, and machine learning creates a comprehensive computational virology framework for managing antiviral resistance in veterinary medicine.

References

[1] Diseases of Poultry, 14th Edition. Wiley-Blackwell.

[2] Merck Veterinary Manual, 11th Edition. Merck & Co.

[3] Webster RG, et al. Evolution and ecology of influenza A viruses. Microbiol Rev. 1992;56(1):152-179.

[4] Russell RJ, et al. The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design. Nature. 2006;443(7107):45-49.

[5] Zwanzig RW. High-temperature equation of state by a perturbation method. I. Nonpolar gases. J Chem Phys. 1954;22(8):1420-1426.

[6] Air GM. Influenza neuraminidase. Influenza Other Respir Viruses. 2012;6(4):245-256.

[7] Colman PM, et al. Structure of the catalytic and antigenic sites in influenza virus neuraminidase. Nature. 1983;303(5912):41-44.

[8] von Itzstein M, et al. Rational design of potent sialidase-based inhibitors of influenza virus replication. Nature. 1993;363(6428):418-423.

[9] Gubareva LV, et al. Influenza virus neuraminidase inhibitors. Lancet. 2000;355(9206):827-835.

[10] Hurt AC, et al. Antiviral resistance during the 2009 influenza A H1N1 pandemic: public health, laboratory, and clinical perspectives. Lancet Infect Dis. 2012;12(3):240-248.

[11] McKimm-Breschkin JL. Influenza neuraminidase inhibitors: antiviral action and mechanisms of resistance. Influenza Other Respir Viruses. 2013;7(Suppl 1):25-36.

[12] Samson M, et al. Influenza virus resistance to neuraminidase inhibitors. Antiviral Res. 2013;98(2):174-185.

[13] Kirkwood JG. Statistical mechanics of fluid mixtures. J Chem Phys. 1935;3(5):300-313.

[14] Jorgensen WL, et al. Free energy perturbation calculations. Acc Chem Res. 1989;22(5):184-189.

[15] Bennett CH. Efficient estimation of free energy differences from Monte Carlo data. J Comput Phys. 1976;22(2):245-268.

[16] Tembe BL, McCammon JA. Ligand-receptor interactions. Comput Chem. 1984;8(4):281-283.

[17] Kollman PA. Free energy calculations: applications to chemical and biochemical phenomena. Chem Rev. 1993;93(7):2395-2417.

[18] Varghese JN, et al. Structure of the influenza virus glycoprotein antigen neuraminidase at 2.9 A resolution. Nature. 1983;303(5912):35-40.

[19] Berman HM, et al. The Protein Data Bank. Nucleic Acids Res. 2000;28(1):235-242.

[20] Olsson MHM, et al. PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions. J Chem Theory Comput. 2011;7(2):525-537.

[21] Wang J, et al. Development and testing of a general amber force field. J Comput Chem. 2004;25(9):1157-1174.

[22] Maier JA, et al. ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. J Chem Theory Comput. 2015;11(8):3696-3713.

[23] Jorgensen WL, et al. Comparison of simple potential functions for simulating liquid water. J Chem Phys. 1983;79(2):926-935.

[24] Joung IS, Cheatham TE. Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations. J Phys Chem B. 2008;112(30):9020-9041.

[25] Essmann U, et al. A smooth particle mesh Ewald method. J Chem Phys. 1995;103(19):8577-8593.

[26] Press WH, et al. Numerical Recipes in C. Cambridge University Press. 1992.

[27] Berendsen HJC, et al. Molecular dynamics with coupling to an external bath. J Chem Phys. 1984;81(8):3684-3690.

[28] Pastor RW, et al. Langevin dynamics of peptides: the frictional dependence of isomerization rates of N-acetylalanyl-N'-methylamide. Biopolymers. 1988;27(4):679-700.

[29] Ryckaert JP, et al. Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J Comput Phys. 1977;23(3):327-341.

[30] Wang L, et al. Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field. J Am Chem Soc. 2015;137(7):2695-2703.

[31] Steinbrecher T, et al. Free energy perturbation calculations of the thermodynamics of protein side-chain mutations. J Mol Biol. 2007;371(4):1119-1130.

[32] Shirts MR, et al. Extremely precise free energy calculations of amino acid side chain analogs: comparison of common molecular mechanics force fields for proteins. J Chem Phys. 2003;119(11):5740-5761.

[33] Shirts MR, Chodera JD. Statistically optimal analysis of samples from multiple equilibrium states. J Chem Phys. 2008;129(12):124105.

[34] Chipot C, Pohorille A. Free Energy Calculations: Theory and Applications in Chemistry and Biology. Springer. 2007.

[35] Okomo-Adhiambo M, et al. Detection of influenza A virus neuraminidase inhibitor resistance by pyrosequencing. Antiviral Res. 2010;85(2):354-360.

[36] Nguyen HT, et al. Neuraminidase inhibitor resistance in influenza viruses: laboratory surveillance. J Infect Dis. 2012;206(10):1484-1491.

[37] Wetherall NT, et al. Evaluation of neuraminidase enzyme assays using different substrates to measure susceptibility of influenza virus to neuraminidase inhibitors. J Clin Microbiol. 2003;41(2):742-750.

[38] Hauser K, et al. Predicting resistance of influenza A virus to neuraminidase inhibitors using free energy perturbation. J Chem Inf Model. 2018;58(10):2092-2100.

[39] de Ruiter A, Oostenbrink C. Free energy calculations of protein-ligand interactions. Curr Opin Chem Biol. 2011;15(4):547-552.

[40] WHO. Global Influenza Surveillance and Response System. World Health Organization.

[41] Bloom JD, et al. Permissive secondary mutations enable the evolution of influenza oseltamivir resistance. Science. 2010;328(5983):1272-1275.

[42] Li KS, et al. Genesis of a highly pathogenic and potentially pandemic H5N1 influenza virus in eastern Asia. Nature. 2004;430(6996):209-213.

[43] Renzette N, et al. Evolution of the influenza A virus genome during infection. mBio. 2014;5(4):e01364-14.

[44] Yen HL, et al. Importance of neuraminidase active-site residues to the neuraminidase inhibitor resistance of influenza A viruses. J Virol. 2006;80(17):8787-8795.

[45] Collins PJ, et al. Crystal structures of oseltamivir-resistant influenza virus neuraminidase mutants. Nature. 2008;453(7199):1258-1261.

[46] Pettersen EF, et al. UCSF Chimera: a visualization system for exploratory research and analysis. J Comput Chem. 2004;25(13):1605-1612.

[47] Sehnal D, et al. Mol*: towards a common library and platform for molecular graphics. J Cheminform. 2018;10(1):4.

[48] Perez A, et al. Advances in free-energy-based computational methods for drug design. Expert Opin Drug Discov. 2016;11(6):571-582.

[49] Mobley DL, et al. Predicting absolute binding free energies to a simple model site. J Mol Biol. 2007;371(4):1118-1134.

[50] Shaw DE, et al. Atomic-level characterization of the structural dynamics of proteins. Science. 2010;330(6002):341-346.

[51] Laio A, Parrinello M. Escaping free-energy minima. Proc Natl Acad Sci USA. 2002;99(20):12562-12566.

[52] Goh GB, et al. Deep learning for computational chemistry. J Comput Chem. 2017;38(16):1291-1307.

[53] Lindorff-Larsen K, et al. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins. 2010;78(8):1950-1958. *** 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.