Despite the success of these retrospective analyses, prospective prediction of H5N1 mutations remains a much more difficult task. Influenza hemagglutinin is heavily glycosylated ¹¹, and the viral glycans can affect both affinity and specificity for host glycans ^12-14. Even using simplified ligands, large-scale expression and experimental screening of hemagglutinin glycoprotein mutants for specificity changes remain challenging due to biosafety issues and the difficulty of doing large-scale mutagenesis in cell-culture systems that will produce the relevant hemagglutinin glycosylation patterns. We have therefore designed an approach to large-scale computational screening of hemagglutinin mutants that will allow more directed experimental validation.

A number of experimental and computational methods have been developed to examine receptor-binding-domain mutants, but ligand-binding mutants from clinical isolates of influenza virus encompass both receptor-binding-domain and allosteric sites. Experimental data for allosteric sites are particularly sparse due to the challenges of high-throughput mutagenesis and screening of influenza hemagglutinin. We have designed a molecular-dynamics approach to score potential mutants with robust predictive power for both receptor-binding-domain and allosteric mutations. We perform thousands of simulations of 17 hemagglutinin mutants totaling >1 millisecond in length and employ a Bayesian model to rank mutations that disrupt hemagglutinin-ligand complex stability. Based on our analysis, we predict a significantly increased k_off for 7 of these mutants. This means of analyzing molecular-dynamics data to make experimentally verifiable predictions offers a potentially general method to identify ligand-binding mutants, particularly allosteric ones. Our analysis provides a robust means to evaluate mutants prior to experimental mutagenesis and testing; these results also constitute an important step towards understanding the determinants of ligand binding by H5N1 influenza.

Dissociation rates were chosen as a means to evaluate predicted ligand-binding mutants because the association and dissociation rates (k_on and k_off) of ligands from wild-type hemagglutinins is relatively slow—data on monovalent k_on and k_off are not available, but X-31 hemagglutinin rosettes bind fetuin with multivalent rates of k_on = 2× 10³ M⁻¹ s⁻¹ and k_off = 2× 10⁻⁴ s⁻¹.¹⁵ Experimental dissociation rates reported for hemagglutinin vary by up to 10,000-fold, however, based in part on the surface conjugation^15,16. Depending on whether a ligand-binding mutation alters the transition-state free energy or only the free energy of the bound state, it would alter both k_on and k_off or k_on alone. We can sample and estimate fast processes more accurately via molecular dynamics than we can slow ones, so acceleration of k_off is a more accessible parameter than deceleration of k_on. Computational methods to predict free energies of binding under active development,¹⁷ but predicting binding of charged, flexible ligands by a flexible protein is extremely challenging for current methods. Methods based on molecular mechanics-Generalized Born calculations have recently been applied to predict hemagglutinin glycan binding.¹⁸ These show promise for predicting receptor-binding-domain mutations, while our molecular-dynamics-based calculations are designed also to detect allosteric mutants in a robust fashion.

The estimated dissociation rates incorporate both positive (dissociation observed) and negative simulation data of varying length; see the supplement for k_off probability density functions of each mutant and starting conformation. We predict the probability that each mutation accelerates k_off, destabilizing the bound complex (Fig. S2). We then perform a bootstrap analysis to identify mutations that significantly speed k_off compared to the wild-type VN1194 complex.

We predict that 7 mutations significantly (p<10⁻⁵) perturb binding of α2,3-sialyllactose by hemagglutinin: E231V, G134V, S136A, N224R, N224V, L48P, and Y82V (Fig. 2). More dissociation events were observed from starting conformation C than A or B, suggesting that the three protein-ligand complexes derived from the crystal structure differ in baseline stability. The effects of each mutation were relatively consistent across starting conformations, however—all but one of these mutants (Y82V) significantly perturb binding in at least two of the three conformations tested. This conformation dependence emphasizes the need for extensive molecular dynamics sampling; we have performed additional simulations of six mutants using 30 additional starting states (Fig. S3).

The three mutations most strongly predicted to destabilize ligand binding by hemagglutinin are S136A, N224V, and L48P. As shown in Fig. 1, S136 lies within the binding pocket, N224 is a loop residue at the edge of the ligand-binding pocket, and L48 is distant from the binding pocket and does not directly contact the bound ligand. The mutations we predict to affect the stability of ligand binding thus include both binding-pocket and allosteric sites. S136A and N224V lie adjacent to experimentally identified ligand-specificity mutation sites 138, 225, 226, 227, and 228. The E231V mutant, which scored fifth in our analysis, lies within the monomer-monomer interface; it is possible that this mutation acts to alter interactions between monomers.

As negative controls, we have tested two mutants, S127P and N197K, where the EC50 in a plate-binding assay using live whole virus is within 10-fold of wild-type VN1194¹⁰ and we thus expect the monomeric k_off to be minimally perturbed. Our analysis does not predict acceleration of either k_off (Fig. S4). In addition, one of our three top-scoring mutants, S136A, has been tested in H3N2 influenza. S136 is conserved between VN1194 and human H3N2 strains; the S136A mutant of A/Aichi/68 has been studied experimentally and has 30% of wild-type activity in an erythrocyte binding assay.²⁷

Here, we present a new means to predict ligand-binding affinity mutations in influenza hemagglutinin. Our method combines molecular dynamics analysis with sequence data and employs information-theoretic methods to score mutation sites based on their relatedness to ligand binding in a robust manner. We combine this scoring with a sequence-based analysis of residue co-variation to produce a focused set of mutation sites. In this report, we predict mutants with altered binding affinity, but the underlying method is designed to be applicable to predicting altered binding specificity as well.

We also demonstrate a new statistical methodology for computational evaluation of ligand-binding mutants by estimating changes to k_off. Dissociation is a more accessible process than association, and it may also be more relevant, as studies of binding of small flexible ligands by MHC molecules²⁸ and RNaseS²⁹ show dissociation rates to be sensitive to mutation while association rates remain relatively constant. This occurs when the transition state is dominated by nonspecific interactions. The bound state free energy is then more mutation-sensitive than the transition state, so k_off is affected much more than k_on. Though plausible, it is unknown if hemagglutinin-ligand binding has such a transition state.

We model the dissociation reaction as approximately two-state kinetically, and a Bayesian framework allows rate estimation based on both positive and negative observations of dissociation in a heterogeneous dataset of many molecular dynamics simulations. This approach is particularly helpful in comparing two rates, as one can utilize the entire probability distribution rather than only the maximum-likelihood estimate. To account for small-sample-size effects, we encapsulate the Bayesian analysis in a bootstrap error analysis to give robust estimates of statistical significance. We have tested this analytic procedure by retrospective validation against mutants that bind similarly to wild-type hemagglutinin in experimental assays and comparing our top-scoring mutant to experimental data on the analogous mutation in H3N2 influenza.

1_si_001