Dynamic residue interaction network analysis of the oseltamivir binding site of N1 neuraminidase and its H274Y mutation site conferring drug resistance in influenza A virus

Initial structures

Structure preparation and MD simulations were performed using the Amber 20 package (Case et al., 2020). The crystal structures of wild-type (WT) avian influenza virus A/H5N1 NA and its H274Y mutant in complex with OTV were obtained from the Protein Data Bank (PDB code: 2HU4 and 3CL0) (Russell et al., 2006; Collins et al., 2008), using a single monomer for modeling. H5N1 NA contains one calcium ion that is necessary for structural stability (Smith et al., 2006). As no calcium ions were found in the crystal structure of H5N1 NA registered as 2HU4, the coordinates of the calcium ion were obtained from the corresponding structure registered as 3CL0. The protonation state of His in the modeled complex at pH 7 was determined using the PDB2PQR server (Dolinsky et al., 2004), and the other ionized residues, that is, Arg, Lys, Asp, and Glu, were treated as charged entities. The LEaP module in Amber 20 was used to add missing hydrogen atoms to the proteins and OTV. H5N1 NA contains eight disulfide bonds. For each disulfide bond, the “bond” command was executed in the tLEaP program to create a covalent bond between the SG atoms of the proximate cysteine residues. The FF14SB variant of the AMBER force field was used to describe the proteins (Maier et al., 2015). The generalized AMBER force field was applied to OTV (Wang et al., 2004). The partial atomic charges in OTV were determined according to the restrained electrostatic potential fitting procedure (Bayly et al., 1993), based on the quantum chemistry calculations at the HF/6-31G(d) level with the Gaussian 16 program (Frisch et al., 2016). The complexes of the WT and H274Y mutant NA with OTV were solvated in a truncated octahedral box of TIP3P water molecules with a distance of at least 10 Å around them. The total charges of the complexes were neutralized by adding sodium counter ions. The total number of atoms in the complexes of the WT and H274Y mutant NA were 30,252 and 29,635, respectively.

Molecular dynamics (MD) simulations

Each system was energy minimized using the steepest descent method for 500 steps, followed by the conjugate gradient method for 4,500 steps, with a harmonic restraint weight of 10 kcal mol⁻¹ Å⁻² on the complexes, except for hydrogen atoms. After energy minimization, each system was gradually heated to 300 K over a period of 200 ps in the NVT ensemble. All MD simulations were performed using the PMEMD module in the AMBER 20 package. Temperature was regulated using the weak-coupling algorithm (Berendsen, Postma & Funsteren, 1984). All bond lengths including hydrogen atoms were constrained using the SHAKE algorithm (Ryckaert, Ciccotti & Berendsen, 1977), allowing an MD time-step of 2 fs to be used. Periodic boundary conditions were adopted. A cutoff of 8 Å was set for non-bonded interactions. Long-range electrostatics were treated using the particle-mesh Ewald method (Darden, York & Pedersen, 1993). After heating, MD simulations were performed for 80 ns in the NpT ensemble at a pressure of 1.0 atm and a temperature of 300 K. The pressure was maintained using a Berendsen barostat. The production phase to be analyzed was the last 40 ns of MD simulations, which was determined based on the root mean square displacement (RMSD) for the backbone atoms in the proteins with respect to the initial structure along the simulation time. The time series of RMSD for the backbone atoms in the WT and H274Y mutant NA are shown in Fig. 1. The changes in RMSD were almost constant after 40 ns, indicating that the MD simulations properly converged in this time region.

Binding free energy calculations

Binding free energies were estimated for OTV complexed with WT and H274Y mutant NA using the Molecular Mechanics Poisson Boltzmann Surface Area (MM-PBSA) method in the Amber 20 package, MMPBSA.py (Miller et al., 2012). The electrostatic contribution to the solvation free energy was estimated using the adaptive Poisson Boltzmann (PB) solver (Baker et al., 2001). The dielectric constants in the protein and water were set to 4 and 80, respectively. A relatively large dielectric constant is desirable in NA, considering that its binding site contains many charged residues (Hou et al., 2011). The ionic strength was set at 150 mM. The ratio between the longest dimension of the rectangular finite-difference grid and that of the solute was set to 4. The linear PB equation was solved using a maximum of 1000 iterations. MM-PBSA calculations were performed over 2,000 frames extracted from the production phase in the last 40 ns of MD simulations.

The entropy due to the vibrational degrees of freedom was calculated by normal mode analysis of 100 configurations using the mmpbsa_py_nabnmode program in the Amber 20 package. Each configuration was energy minimized with a generalized Born solvent model, using a maximum of 10,000 steps with a target root-mean-square gradient of 10⁻³ kcal mol⁻¹ Å⁻¹.

Dynamic residue interaction network (dRIN) analysis

Residue interactions were analyzed for 2,000 three-dimensional structures of the WT and H274Y mutant NA-OTV complexes generated in the 40-ns production phase of the MD simulations. In this study, ligand-residue interactions were included in the residue interactions. The RING software was used to categorize residue-residue and residue-ligand interactions into specific types, such as hydrogen bonds, van der Waals (vdW) interactions, disulfide bonds, salt bridges, π-cation interactions, and π–π stacking interactions (Piovesan, Minervini & Tosatto, 2016). As the RING software missed hydrogen bonds between residues and the ligand in the NA-OTV complexes, the cpptraj program in the Amber 20 package was used to detect hydrogen bonds. Although the RING software identified several different types of interactions per residue pair, only one interaction per interaction type was considered. We further investigated the interactions between residues that were not included in the attractive interactions identified by the RING software but were in close contact with each other. In this study, two residues were identified as being in close contact, if $d_{ij}-\left({\sigma }_i+{\sigma }_j\right)$ < 0.4 Å, where d_ij denotes the interatomic distance between the i- and j-th atoms in the two residues and σ_i is the vdW radius of the i-th atom. This definition of close contact has been adopted by several programs for the visualization and analysis of molecular structures, such as the UCSF Chimera software (Pettersen et al., 2004). By repeatedly applying the above procedure, the dRINs were constructed based on the sets of residue interactions and close contacts obtained from the MD simulations for the WT and H274Y mutant NA-OTV complexes. The occupancy of residue interactions was calculated as the percentage of frames with interactions between residues in the MD simulation. The dRINs were visualized using Cytoscape software (Shannon et al., 2003) with residues and their interactions represented by nodes and edges, respectively.

Binding structures and energies

Figure 2 represents snapshot images obtained from the MD simulations for the WT and H274Y mutant NA-OTV complexes, showing the OTV binding site and the region adjacent to residue 274. Table 1 summarizes the computed binding free energies (ΔG) of OTV for the WT and H274Y mutant NA obtained from the MM-PBSA calculations, along with the enthalpy (ΔH) and entropy (TΔS). The detailed results for the energetic components divided by each intermolecular interaction are shown in Table S1. The binding free energies of OTV were computed to be −11.54 and −4.34 kcal mol⁻¹ for the WT and H274Y mutant NA, respectively. The 7.20 kcal mol⁻¹ increase in the binding free energy of OTV due to the H274Y mutation could significantly reduce the efficiency of this inhibitor against NA. This is supported by the experimental fact that the H274Y mutant NA has a 300- to 1700-fold decrease in relative susceptibility to OTV compared with the WT NA in H5N1 viruses (Ives et al., 2002; Wang, Tai & Mendel, 2002; Kiso et al., 2004; De Jong et al., 2005; Mishin, Hayden & Gubareva, 2005; Yen et al., 2007). Several computational studies have been conducted on the change in the binding free energy of OTV due to the H274Y mutation (Malaisree et al., 2009; Wang & Zheng, 2009; Park & Jo, 2009; Nguyen, Mai & Li, 2011; Li et al., 2012; Vergara-Jaque et al., 2012; Woods et al., 2012; Ripoll et al., 2012; Woods et al., 2013; Yusuf et al., 2016). The current results are qualitatively consistent with previous experimental and computational studies, indicating that the MD simulations, which form the basis for the subsequent analyses, were reliable.

The binding free energy computed using the MM-PBSA approach can be sensitive to numerous factors, such as the value of the dielectric constant (ε) chosen for the protein. Generally, ε = 1.0 is used for proteins; however, in this study, ε = 4.0 was used for NA based on previous studies, considering that the binding site of NA contains many charged residues (Hou et al., 2011). The results of binding free energies obtained using the MM-PBSA method with ε = 1.0 to confirm the validity of this factor are shown in Table S2. The results calculated with ε = 1.0 showed a larger difference in binding free energies (10.48 kcal mol⁻¹) than those calculated with ε = 4.0 (7.20 kcal mol⁻¹), further reducing the affinity of NA for OTV after the H274Y mutation. The difference in the binding free energy of the WT and H274Y mutant NA to OTV was calculated to be approximately 4.4 kcal mol⁻¹, based on the experimentally measured IC₅₀ value (Yen et al., 2007). This clearly indicates that it is better to use ε = 4.0 for NA in the MM-PBSA calculation, as mentioned in the previous study (Hou et al., 2011).

To further analyze the structural changes associated with the H274Y mutation, the pocket cavity volumes of the WT and mutant NA were examined using the POVME 3.0 software (Wagner et al., 2017). The results analyzed for the last 40 ns of the MD simulations showed that the average pocket cavity volumes were computed to be 394 and 853 ų, for the WT and H274Y mutant NA, respectively. The pocket cavity volume of the H274Y mutant NA was increased more than twofold compared to that of the WT NA. In the following sections, such structural changes in the drug-binding site of NA have been analyzed in more detail in terms of residue interactions.

Dynamic residue interaction network (dRIN) analysis

Figure 3 shows the dRIN graphs for the WT and H274Y mutant NA complexed with OTV. In a dRIN graph, a node represents an amino acid residue or a ligand, and an edge connecting two nodes represents a residue-residue or a residue-ligand interaction. In most current approaches, RIN is modeled using a single structure (Piovesan, Minervini & Tosatto, 2016; Shcherbinin & Veselovsky, 2019). dRIN extends the original version as it models using multiple structures, providing a statistical perspective on residue interactions. The thickness of the edge represents the occupancy at which an interaction occurs between residues. In Fig. 3, only residue interactions with total occupancy > 10% are shown.

Figure 4A shows the occupancies at which some interactions occur between residues. Occupancies for each type of residue interaction are summarized in Tables 2 and 3. Figure 4B shows the change in the occupancies of residues interacting with each of the WT and H274Y mutant NA.

Figure 5A shows the decomposed protein-ligand complex binding free energies of the WT and H274Y mutant NA on a per-residue basis according to the MM-PBSA calculations. Detailed information on the analysis of the protein-ligand interaction, divided by component per residue, is summarized in Tables 4 and 5. Figure 5B shows the change in the decomposed binding free energies between the WT and H274Y mutant NA.

dRIN of WT NA

Figure 3A shows that the OTV binding site in WT NA has five charged residues, E119, D151, R152, R292, and R371, which form strong hydrogen-bond interactions with OTV. Of these, the hydrogen-bond interaction with three charged residues, E119, R292, and R371 were very robust, accounting for 100% of the interactions formed; the hydrogen-bond interaction with the remaining two of the five charged residues, D151 and R152, have slightly lower occupancies, around 90% (Fig. 4A and Table 2). As shown in Fig. 5A, the major contributors to the binding free energy of OTV for WT NA were E119, D151, R292, and R371. Table 4 shows that the main component of the contribution from these residues to the binding free energy of OTV for WT NA was the electrostatic interaction. The residues S246 and E276 were in close contact with OTV at occupancies of around 50%, and no specific type of attractive interaction was identified for OTV at these residues (Table 2). The main interaction between Y406 and OTV is also close contact (35%), but additionally form a small amount of vdW interaction (9%) (Table 2). Figure 5A shows that these three closely contacting residues, S246, E276, and Y406 contributed little to the binding free energy of OTV for WT NA.

In the adjacent site of residue 274 shown in Fig. 3A, the aromatic residues, Y252 and H296 formed π–π stacking interactions, and W295 formed vdW interaction with H274. Although the edges corresponding to these residue-residue pairs in Fig. 3A only represented major interactions, Table 2 shows that the contributions of both the vdW and π–π stacking interactions were significant among these aromatic residues. The π–π stacking interaction between Y252 and H274 was very robust, accounting for 100% of the interactions formed, whereas the formation of vdW interactions was also observed with an occupancy of about 50% (Table 2). The formation of both vdW and π–π stacking interactions was observed between H296 and H274, but the occupancy of π–π stacking interaction was 7% larger than that of vdW interaction (Table 2). The vdW interaction was mainly observed between W295 and H274 with an occupancy of about 90%, but the π–π stacking interaction was also formed between them with an occupancy of about 50%, in addition to a small amount of hydrogen-bond interaction (12%) (Table 2). The residues N294 and Y316 formed hydrogen-bond interactions with H274 (Fig. 3A). The hydrogen-bond interaction between N294 and H274 was very robust, with nearly 100% of the interaction formed, and the vdW interactions were also frequently observed between them with an occupancy of 85% (Table 2). The occupancy of hydrogen-bond interaction between Y316 and H274 was about 50% (Table 2). The residue P301 was mainly in close contact with H274 (72%), and also formed a weak vdW interaction (17%) (Table 2). The formation of both vdW interaction and close contact were observed between D293 and H274, but the occupancy of vdW interaction was about 10% larger than that of close contact (Table 2). The residue S246 formed interactions with H274, but in much smaller fractions (Table 2).

As shown in Figure 3A, the dRIN graph indicates that residue 274 was not observed to interact directly with OTV. Therefore, the decrease in binding affinity between OTV and NA due to the H274Y mutation can be the result of indirect effects caused by changes in the residue interaction network.

dRIN of H274Y mutant NA

As shown in Fig. 3B, the characteristics of dRIN in the H274Y mutant NA changed significantly as compared to the dRIN in the WT NA. Figure 4B indicates that the occupancy of ligand-residue interactions between the OTV and its surrounding residues, E119, D151, R152, and S246 decreased after H274Y mutation, whereas the occupancy of the OTV-Y406 interaction increased. In particular, it was observed that the occupancies of interactions with OTV were greatly reduced by more than 80% for the residues D151 and R152, and could be expected to contribute significantly to the decreased binding affinity of OTV with H274Y mutant NA. As shown in Fig. 5B, the H274Y mutation resulted in significantly increased components of decomposed binding free energy at E119, D151, and R152, destabilizing the interaction between OTV and H274Y mutant NA. D151 and R152 are residues located in the 150-loop region in NA. Previous studies have reported that the H274Y mutation in NA can change the structure of the 150-loop region from a closed to an open form (Kar & Knecht, 2012). Thus, the H274Y mutation reduced the binding affinity of OTV to NA by cleaving the ligand-residue interaction.

The dRIN structure of the H274Y mutation site remained almost unchanged, even after the His-to-Tyr mutation, except for a few residue-residue pairs and the interface bordering the OTV (Fig. 3B). In particular, two histidine moieties overlapped with each other through vdW and π–π stacking interactions in the pair of H296 and H274 in WT NA, whereas no interaction was observed between the two aromatic rings in H274Y mutant NA (Fig. 3 and Tables 2 and 3). Notably, high proportions of both hydrogen-bond and vdW interactions were observed for the pair of N294 and H274 in WT NA, as shown in Table 2. However, the H274Y mutation resulted in decreased hydrogen-bond interaction occupancy compared to that in the WT NA, thereby making vdW interaction more prevalent (Tables 2 and 3). Similarly, for the interaction between D293 and residue 274, the overall occupancy of the interactions was about 80% in both the WT and H274Y mutant NA, but the preferred interaction in WT NA was the vdW interaction, and that in the H274Y mutant NA was close contact (Tables 2 and 3).

As shown in Fig. 3B, there were three residues, S246, E276, and R292, at the interface between the OTV binding site of NA and its H274Y mutation site, that were observed to interact directly with both OTV and residue 274 in the H274Y mutant NA. His-to-Tyr mutation of residue 274 caused significant changes in the dRIN at the interface region. As shown in Fig. 4B, the occupancies of interactions between residue 274 and the residues S246, E276, and R292, which form the interface, increased after the H274Y mutation, whereas the occupancies of interactions with OTV remained almost unchanged. Notably, no interaction was observed between E276 and H274 in WT NA; however, in the H274Y mutant NA, a pair of residues E276 and Y274 formed a solid hydrogen-bond interaction with an occupancy of about 100% (Fig. 3 and Tables 2 and 3). On the contrary, there were close contacts between E276 and OTV, with an occupancy of approximately 50% in WT NA, and this occupancy remained unchanged after the H274Y mutation (Fig. 3 and Tables 2 and 3). In the pair of residues 246 and 274, weak interactions were observed between S246 and H274, with an overall occupancy of approximately 23% in WT NA, and H274Y mutation increased close contacts to approximately 80% and an overall occupancy to approximately 90% (Tables 2 and 3). In contrast, the occupancy of close contacts between S246 and OTV was 54% in WT NA; this occupancy decreased after H274Y mutation, but the difference was small, about 17% (Tables 2 and 3). The pair of residues 292 and 274 in WT NA had no interaction between R292 and H274, but that in H274Y mutant NA had close contact with an occupancy of about 40% (Tables 2 and 3). On the contrary, R292 formed a robust hydrogen-bond interaction with OTV with an occupancy of 100% in the WT NA, and this occupancy remained unchanged in the H274Y mutant NA (Tables 2 and 3).